Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.vector;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.LogManager;
012import org.apache.logging.log4j.Logger;
013
014import edu.jas.kern.PrettyPrint;
015import edu.jas.structure.AlgebraElem;
016import edu.jas.structure.NotInvertibleException;
017import edu.jas.structure.RingElem;
018
019
020/**
021 * GenMatrix implements a generic matrix algebra over RingElem entries. Matrix
022 * has n columns and m rows over C.
023 * @author Heinz Kredel
024 */
025
026public class GenMatrix<C extends RingElem<C>> implements AlgebraElem<GenMatrix<C>, C> {
027
028
029    private static final Logger logger = LogManager.getLogger(GenMatrix.class);
030
031
032    public final GenMatrixRing<C> ring;
033
034
035    public final ArrayList<ArrayList<C>> matrix;
036
037
038    private int hashValue = 0;
039
040
041    /**
042     * Constructor for zero GenMatrix.
043     * @param r matrix ring
044     */
045    public GenMatrix(GenMatrixRing<C> r) {
046        this(r, r.getZERO().matrix);
047    }
048
049
050    /**
051     * Constructor for GenMatrix.
052     * @param r matrix ring
053     * @param m matrix
054     */
055    public GenMatrix(GenMatrixRing<C> r, List<List<C>> m) {
056        ring = r;
057        matrix = new ArrayList<ArrayList<C>>(r.rows);
058        for (List<C> row : m) {
059            ArrayList<C> nr = new ArrayList<C>(row);
060            matrix.add(nr);
061        }
062        logger.info("{} x {} matrix constructed", ring.rows, ring.cols);
063    }
064
065
066    /**
067     * Constructor for GenMatrix.
068     * @param r matrix ring
069     * @param m matrix
070     */
071    public GenMatrix(GenMatrixRing<C> r, ArrayList<ArrayList<C>> m) {
072        if (r == null || m == null) {
073            throw new IllegalArgumentException("Empty r or m not allowed, r = " + r + ", m = " + m);
074        }
075        ring = r;
076        matrix = new ArrayList<ArrayList<C>>(m);
077        logger.info("{} x {} matrix constructed", ring.rows, ring.cols);
078    }
079
080
081    /**
082     * Constructor for GenMatrix.
083     * @param r matrix ring
084     * @param m matrix
085     */
086    public GenMatrix(GenMatrixRing<C> r, C[][] m) {
087        ring = r;
088        matrix = new ArrayList<ArrayList<C>>(r.rows);
089        for (C[] row : m) {
090            ArrayList<C> nr = new ArrayList<C>(r.cols);
091            for (int i = 0; i < row.length; i++) {
092                nr.add(row[i]);
093            }
094            matrix.add(nr);
095        }
096        logger.info("{} x {} matrix constructed", ring.rows, ring.cols);
097    }
098
099
100    /**
101     * Get element at row i, column j.
102     * @param i row index.
103     * @param j column index.
104     * @return this(i,j).
105     */
106    public C get(int i, int j) {
107        return matrix.get(i).get(j);
108    }
109
110
111    /**
112     * Set element at row i, column j. Mutates this matrix.
113     * @param i row index.
114     * @param j column index.
115     * @param el element to set.
116     */
117    public void setMutate(int i, int j, C el) {
118        ArrayList<C> ri = matrix.get(i);
119        ri.set(j, el);
120        hashValue = 0; // invalidate
121    }
122
123
124    /**
125     * Set element at row i, column j.
126     * @param i row index.
127     * @param j column index.
128     * @param el element to set.
129     * @return new matrix m, with m(i,j) == el.
130     */
131    public GenMatrix<C> set(int i, int j, C el) {
132        GenMatrix<C> mat = this.copy();
133        mat.setMutate(i, j, el);
134        return mat;
135    }
136
137
138    /**
139     * Get column i.
140     * @param i column index.
141     * @return this(*,i) as vector.
142     */
143    public GenVector<C> getColumn(int i) {
144        List<C> cl = new ArrayList<C>(ring.rows);
145        for (int k = 0; k < ring.rows; k++) {
146            cl.add(matrix.get(k).get(i));
147        }
148        GenVectorModul<C> vfac = new GenVectorModul<C>(ring.coFac, ring.rows);
149        GenVector<C> col = new GenVector<C>(vfac, cl);
150        return col;
151    }
152
153
154    /**
155     * Get row i.
156     * @param i row index.
157     * @return this(i,*) as vector.
158     */
159    public GenVector<C> getRow(int i) {
160        List<C> cl = new ArrayList<C>(ring.cols);
161        cl.addAll(matrix.get(i));
162        GenVectorModul<C> vfac = new GenVectorModul<C>(ring.coFac, ring.cols);
163        GenVector<C> row = new GenVector<C>(vfac, cl);
164        return row;
165    }
166
167
168    /**
169     * Get diagonal.
170     * @return diagonal(this) as vector.
171     */
172    public GenVector<C> getDiagonal() {
173        List<C> cl = new ArrayList<C>(ring.rows);
174        for (int i = 0; i < ring.rows; i++) {
175            cl.add(matrix.get(i).get(i));
176        }
177        GenVectorModul<C> vfac = new GenVectorModul<C>(ring.coFac, ring.rows);
178        GenVector<C> dia = new GenVector<C>(vfac, cl);
179        return dia;
180    }
181
182
183    /**
184     * Trace.
185     * @return sum of diagonal elements.
186     */
187    public C trace() {
188        C cl = ring.coFac.getZERO();
189        for (int i = 0; i < ring.rows; i++) {
190            cl = cl.sum(matrix.get(i).get(i));
191        }
192        return cl;
193    }
194
195
196    /**
197     * Get upper triangular U matrix.
198     * @return U as matrix with equal length rows.
199     */
200    public GenMatrix<C> getUpper() {
201        final C zero = ring.coFac.getZERO();
202        final C one = ring.coFac.getONE();
203        List<List<C>> cl = new ArrayList<List<C>>(ring.rows);
204        for (int k = 0; k < ring.rows; k++) {
205            List<C> ul = matrix.get(k);
206            List<C> rl = new ArrayList<C>(ring.cols);
207            for (int i = 0; i < ring.cols; i++) {
208                if (i < k) {
209                    rl.add(zero);
210                } else if (i >= k) {
211                    rl.add(ul.get(i));
212                // } else {
213                //     if (ul.get(i).isZERO()) {
214                //         rl.add(zero);
215                //     } else {
216                //         rl.add(one); // mat(k,k).inverse()
217                //     }
218                }
219            }
220            cl.add(rl);
221        }
222        GenMatrix<C> U = new GenMatrix<C>(ring, cl);
223        return U;
224    }
225
226
227    /**
228     * Get upper triangular U matrix with diagonal 1.
229     * @return U as matrix with equal length rows and diagonal 1.
230     */
231    public GenMatrix<C> getUpperScaled() {
232        final C zero = ring.coFac.getZERO();
233        final C one = ring.coFac.getONE();
234        List<List<C>> cl = new ArrayList<List<C>>(ring.rows);
235        for (int k = 0; k < ring.rows; k++) {
236            List<C> ul = matrix.get(k);
237            C kk = ul.get(k);
238            if (kk.isZERO()) {
239                kk = one;
240            } else {
241                kk = kk.inverse();
242            }
243            List<C> rl = new ArrayList<C>(ring.cols);
244            for (int i = 0; i < ring.cols; i++) {
245                if (i < k) {
246                    rl.add(zero);
247                } else { // if (i >= k)
248                    rl.add( ul.get(i).multiply(kk) );
249                }
250            }
251            cl.add(rl);
252        }
253        GenMatrix<C> U = new GenMatrix<C>(ring, cl);
254        return U;
255    }
256
257
258    /**
259     * Get lower triangular L matrix.
260     * @return L as matrix with equal length rows.
261     */
262    public GenMatrix<C> getLower() {
263        final C zero = ring.coFac.getZERO();
264        List<List<C>> cl = new ArrayList<List<C>>(ring.rows);
265        for (int k = 0; k < ring.rows; k++) {
266            List<C> ul = matrix.get(k);
267            List<C> rl = new ArrayList<C>(ring.cols);
268            for (int i = 0; i < ring.cols; i++) {
269                if (i <= k) {
270                    rl.add(ul.get(i)); // can be zero
271                } else if (i > k) {
272                    rl.add(zero);
273                }
274            }
275            cl.add(rl);
276        }
277        GenMatrix<C> L = new GenMatrix<C>(ring, cl);
278        return L;
279    }
280
281
282    /**
283     * Get the String representation as RingElem.
284     * @see java.lang.Object#toString()
285     */
286    @Override
287    public String toString() {
288        StringBuffer s = new StringBuffer();
289        boolean firstRow = true;
290        s.append("[\n");
291        for (List<C> val : matrix) {
292            if (firstRow) {
293                firstRow = false;
294            } else {
295                s.append(",\n");
296            }
297            boolean first = true;
298            s.append("[ ");
299            for (C c : val) {
300                if (first) {
301                    first = false;
302                } else {
303                    s.append(", ");
304                }
305                s.append(c.toString());
306            }
307            s.append(" ]");
308        }
309        s.append(" ] ");
310        if (!PrettyPrint.isTrue()) {
311            s.append(":: " + ring.toString());
312            s.append("\n");
313        }
314        return s.toString();
315    }
316
317
318    /**
319     * Get a scripting compatible string representation.
320     * @return script compatible representation for this Element.
321     * @see edu.jas.structure.Element#toScript()
322     */
323    @Override
324    public String toScript() {
325        // Python case
326        StringBuffer s = new StringBuffer();
327        boolean firstRow = true;
328        s.append("( ");
329        for (List<C> val : matrix) {
330            if (firstRow) {
331                firstRow = false;
332            } else {
333                s.append(", ");
334            }
335            boolean first = true;
336            s.append("( ");
337            for (C c : val) {
338                if (first) {
339                    first = false;
340                } else {
341                    s.append(", ");
342                }
343                s.append(c.toScript());
344            }
345            s.append(" )");
346        }
347        s.append(" ) ");
348        return s.toString();
349    }
350
351
352    /**
353     * Get a scripting compatible string representation of the factory.
354     * @return script compatible representation for this ElemFactory.
355     * @see edu.jas.structure.Element#toScriptFactory()
356     */
357    @Override
358    public String toScriptFactory() {
359        // Python case
360        return factory().toScript();
361    }
362
363
364    /**
365     * Get the corresponding element factory.
366     * @return factory for this Element.
367     * @see edu.jas.structure.Element#factory()
368     */
369    public GenMatrixRing<C> factory() {
370        return ring;
371    }
372
373
374    /**
375     * Copy method.
376     * @see edu.jas.structure.Element#copy()
377     */
378    @Override
379    public GenMatrix<C> copy() {
380        //return ring.copy(this);
381        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
382        ArrayList<C> v;
383        for (ArrayList<C> val : matrix) {
384            v = new ArrayList<C>(val);
385            m.add(v);
386        }
387        return new GenMatrix<C>(ring, m);
388    }
389
390
391    /**
392     * Stack method.
393     * @param st stacked matrix ring.
394     * @param b other matrix.
395     * @return stacked matrix, this on top of other.
396     */
397    public GenMatrix<C> stack(GenMatrixRing<C> st, GenMatrix<C> b) {
398        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(st.rows);
399        ArrayList<C> v;
400        for (ArrayList<C> val : matrix) {
401            v = new ArrayList<C>(val);
402            m.add(v);
403        }
404        for (ArrayList<C> val : b.matrix) {
405            v = new ArrayList<C>(val);
406            m.add(v);
407        }
408        return new GenMatrix<C>(st, m);
409    }
410
411
412    /**
413     * Concat method.
414     * @param cc concated matrix ring.
415     * @param b other matrix.
416     * @return concated matrix, this before of other.
417     */
418    public GenMatrix<C> concat(GenMatrixRing<C> cc, GenMatrix<C> b) {
419        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(cc.rows);
420        ArrayList<ArrayList<C>> bm = b.matrix;
421        ArrayList<C> v, o;
422        int i = 0;
423        for (ArrayList<C> val : matrix) {
424            v = new ArrayList<C>(val);
425            o = bm.get(i++);
426            v.addAll(o);
427            m.add(v);
428        }
429        return new GenMatrix<C>(cc, m);
430    }
431
432
433    /**
434     * Test if this is equal to a zero matrix.
435     */
436    public boolean isZERO() {
437        for (List<C> row : matrix) {
438            for (C elem : row) {
439                if (!elem.isZERO()) {
440                    return false;
441                }
442            }
443        }
444        return true;
445    }
446
447
448    /**
449     * Test if this is one.
450     * @return true if this is 1, else false.
451     */
452    public boolean isONE() {
453        int i = 0;
454        for (List<C> row : matrix) {
455            int j = 0;
456            for (C elem : row) {
457                if (i == j) {
458                    if (!elem.isONE()) {
459                        //System.out.println("elem.isONE = " + elem);
460                        return false;
461                    }
462                } else if (!elem.isZERO()) {
463                    //System.out.println("elem.isZERO = " + elem);
464                    return false;
465                }
466                j++;
467            }
468            i++;
469        }
470        return true;
471    }
472
473
474    /**
475     * Test if this is a non-zero diagonal matrix.
476     * @return true if this is non-zero diagonal, else false.
477     */
478    public boolean isDiagonal() {
479        int z = 0;
480        int i = 0;
481        for (List<C> row : matrix) {
482            int j = 0;
483            for (C elem : row) {
484                if (i == j) {
485                    if (!elem.isZERO()) {
486                        z++;
487                    }
488                } else if (!elem.isZERO()) {
489                    //System.out.println("elem.isZERO = " + elem);
490                    return false;
491                }
492                j++;
493            }
494            i++;
495        }
496        if (z == 0) { // matrix is zero
497            return false;
498        }
499        return true;
500    }
501
502
503    /**
504     * Comparison with any other object.
505     * @see java.lang.Object#equals(java.lang.Object)
506     */
507    @Override
508    @SuppressWarnings("unchecked")
509    public boolean equals(Object other) {
510        if (!(other instanceof GenMatrix)) {
511            return false;
512        }
513        GenMatrix om = (GenMatrix) other;
514        if (!ring.equals(om.ring)) {
515            return false;
516        }
517        if (!matrix.equals(om.matrix)) {
518            return false;
519        }
520        return true;
521    }
522
523
524    /**
525     * Hash code for this GenMatrix.
526     * @see java.lang.Object#hashCode()
527     */
528    @Override
529    public int hashCode() {
530        if (hashValue == 0) {
531            hashValue = 37 * matrix.hashCode() + ring.hashCode();
532            if (hashValue == 0) {
533                hashValue = 1;
534            }
535        }
536        return hashValue;
537    }
538
539
540    /**
541     * compareTo, lexicogaphical comparison.
542     * @param b other
543     * @return 1 if (this &lt; b), 0 if (this == b) or -1 if (this &gt; b).
544     */
545    @Override
546    public int compareTo(GenMatrix<C> b) {
547        if (!ring.equals(b.ring)) {
548            return -1;
549        }
550        ArrayList<ArrayList<C>> om = b.matrix;
551        int i = 0;
552        for (ArrayList<C> val : matrix) {
553            ArrayList<C> ov = om.get(i++);
554            int j = 0;
555            for (C c : val) {
556                int s = c.compareTo(ov.get(j++));
557                if (s != 0) {
558                    return s;
559                }
560            }
561        }
562        return 0;
563    }
564
565
566    /**
567     * Test if this is a unit. I.e. there exists x with this.multiply(x).isONE()
568     * == true. Tests if matrix is not singular.
569     * Was previously a test if all diagonal elements are units and all other
570     *  elements are zero.
571     * @return true if this is a unit, else false.
572     */
573    public boolean isUnit() {
574        LinAlg<C> la = new LinAlg<C>();
575        GenMatrix<C> mat = this.copy();
576        List<Integer> P = la.decompositionLU(mat);
577        if (P == null || P.isEmpty()) {
578            return false;
579        }
580        return true;
581    }
582
583
584    /**
585     * sign of matrix.
586     * @return 1 if (this &lt; 0), 0 if (this == 0) or -1 if (this &gt; 0).
587     */
588    public int signum() {
589        return compareTo(ring.getZERO());
590    }
591
592
593    /**
594     * Sum of matrices.
595     * @param b other matrix.
596     * @return this+b
597     */
598    public GenMatrix<C> sum(GenMatrix<C> b) {
599        ArrayList<ArrayList<C>> om = b.matrix;
600        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
601        int i = 0;
602        for (ArrayList<C> val : matrix) {
603            ArrayList<C> ov = om.get(i++);
604            ArrayList<C> v = new ArrayList<C>(ring.cols);
605            int j = 0;
606            for (C c : val) {
607                C e = c.sum(ov.get(j++));
608                v.add(e);
609            }
610            m.add(v);
611        }
612        return new GenMatrix<C>(ring, m);
613    }
614
615
616    /**
617     * Difference of matrices.
618     * @param b other matrix.
619     * @return this-b
620     */
621    public GenMatrix<C> subtract(GenMatrix<C> b) {
622        ArrayList<ArrayList<C>> om = b.matrix;
623        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
624        int i = 0;
625        for (ArrayList<C> val : matrix) {
626            ArrayList<C> ov = om.get(i++);
627            ArrayList<C> v = new ArrayList<C>(ring.cols);
628            int j = 0;
629            for (C c : val) {
630                C e = c.subtract(ov.get(j++));
631                v.add(e);
632            }
633            m.add(v);
634        }
635        return new GenMatrix<C>(ring, m);
636    }
637
638
639    /**
640     * Negative of this matrix.
641     * @return -this
642     */
643    public GenMatrix<C> negate() {
644        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
645        //int i = 0;
646        for (ArrayList<C> val : matrix) {
647            ArrayList<C> v = new ArrayList<C>(ring.cols);
648            for (C c : val) {
649                C e = c.negate();
650                v.add(e);
651            }
652            m.add(v);
653        }
654        return new GenMatrix<C>(ring, m);
655    }
656
657
658    /**
659     * Absolute value of this matrix.
660     * @return abs(this)
661     */
662    public GenMatrix<C> abs() {
663        if (signum() < 0) {
664            return negate();
665        }
666        return this;
667    }
668
669
670    /**
671     * Product of this matrix with scalar.
672     * @param s scalar.
673     * @return this*s
674     */
675    public GenMatrix<C> multiply(C s) {
676        return scalarMultiply(s);
677    }
678
679
680    /**
681     * Product of this matrix with scalar.
682     * @param s scalar
683     * @return this*s
684     */
685    public GenMatrix<C> scalarMultiply(C s) {
686        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
687        //int i = 0;
688        for (ArrayList<C> val : matrix) {
689            ArrayList<C> v = new ArrayList<C>(ring.cols);
690            for (C c : val) {
691                C e = c.multiply(s);
692                v.add(e);
693            }
694            m.add(v);
695        }
696        return new GenMatrix<C>(ring, m);
697    }
698
699
700    /**
701     * Left product of this matrix with scalar.
702     * @param s scalar
703     * @return s*this
704     */
705    public GenMatrix<C> leftScalarMultiply(C s) {
706        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
707        //int i = 0;
708        for (ArrayList<C> val : matrix) {
709            ArrayList<C> v = new ArrayList<C>(ring.cols);
710            for (C c : val) {
711                C e = s.multiply(c);
712                v.add(e);
713            }
714            m.add(v);
715        }
716        return new GenMatrix<C>(ring, m);
717    }
718
719
720    /**
721     * Linear compination of this matrix with scalar multiple of other matrix.
722     * @param s scalar
723     * @param t scalar
724     * @param b other matrix.
725     * @return this*s+b*t
726     */
727    public GenMatrix<C> linearCombination(C s, GenMatrix<C> b, C t) {
728        ArrayList<ArrayList<C>> om = b.matrix;
729        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
730        int i = 0;
731        for (ArrayList<C> val : matrix) {
732            ArrayList<C> ov = om.get(i++);
733            ArrayList<C> v = new ArrayList<C>(ring.cols);
734            int j = 0;
735            for (C c : val) {
736                C c1 = c.multiply(s);
737                C c2 = ov.get(j++).multiply(t);
738                C e = c1.sum(c2);
739                v.add(e);
740            }
741            m.add(v);
742        }
743        return new GenMatrix<C>(ring, m);
744    }
745
746
747    /**
748     * Linear combination of this matrix with scalar multiple of other matrix.
749     * @param t scalar
750     * @param b other matrix.
751     * @return this+b*t
752     */
753    public GenMatrix<C> linearCombination(GenMatrix<C> b, C t) {
754        ArrayList<ArrayList<C>> om = b.matrix;
755        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
756        int i = 0;
757        for (ArrayList<C> val : matrix) {
758            ArrayList<C> ov = om.get(i++);
759            ArrayList<C> v = new ArrayList<C>(ring.cols);
760            int j = 0;
761            for (C c : val) {
762                C c2 = ov.get(j++).multiply(t);
763                C e = c.sum(c2);
764                v.add(e);
765            }
766            m.add(v);
767        }
768        return new GenMatrix<C>(ring, m);
769    }
770
771
772    /**
773     * Left linear combination of this matrix with scalar multiple of other
774     * matrix.
775     * @param t scalar
776     * @param b other matrix.
777     * @return this+t*b
778     */
779    public GenMatrix<C> linearCombination(C t, GenMatrix<C> b) {
780        ArrayList<ArrayList<C>> om = b.matrix;
781        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
782        int i = 0;
783        for (ArrayList<C> val : matrix) {
784            ArrayList<C> ov = om.get(i++);
785            ArrayList<C> v = new ArrayList<C>(ring.cols);
786            int j = 0;
787            for (C c : val) {
788                C c2 = t.multiply(ov.get(j++));
789                C e = c.sum(c2);
790                v.add(e);
791            }
792            m.add(v);
793        }
794        return new GenMatrix<C>(ring, m);
795    }
796
797
798    /**
799     * left linear compination of this matrix with scalar multiple of other
800     * matrix.
801     * @param s scalar
802     * @param t scalar
803     * @param b other matrix.
804     * @return s*this+t*b
805     */
806    public GenMatrix<C> leftLinearCombination(C s, C t, GenMatrix<C> b) {
807        ArrayList<ArrayList<C>> om = b.matrix;
808        ArrayList<ArrayList<C>> m = new ArrayList<ArrayList<C>>(ring.rows);
809        int i = 0;
810        for (ArrayList<C> val : matrix) {
811            ArrayList<C> ov = om.get(i++);
812            ArrayList<C> v = new ArrayList<C>(ring.cols);
813            int j = 0;
814            for (C c : val) {
815                C c1 = s.multiply(c);
816                C c2 = t.multiply(ov.get(j++));
817                C e = c1.sum(c2);
818                v.add(e);
819            }
820            m.add(v);
821        }
822        return new GenMatrix<C>(ring, m);
823    }
824
825
826    /**
827     * Transposed matrix.
828     * @param tr transposed matrix ring.
829     * @return transpose(this)
830     */
831    public GenMatrix<C> transpose(GenMatrixRing<C> tr) {
832        GenMatrix<C> t = tr.getZERO().copy();
833        ArrayList<ArrayList<C>> m = t.matrix;
834        int i = 0;
835        for (ArrayList<C> val : matrix) {
836            int j = 0;
837            for (C c : val) {
838                (m.get(j)).set(i, c); //A[j,i] = A[i,j]
839                j++;
840            }
841            i++;
842        }
843        // return new GenMatrix<C>(tr,m);
844        return t;
845    }
846
847
848    /**
849     * Transposed matrix.
850     * @return transpose(this)
851     */
852    public GenMatrix<C> transpose() {
853        GenMatrixRing<C> tr = ring.transpose();
854        return transpose(tr);
855    }
856
857
858    /**
859     * Multiply this with S.
860     * @param S other matrix.
861     * @return this * S.
862     */
863    public GenMatrix<C> multiply(GenMatrix<C> S) {
864        int na = ring.blocksize;
865        int nb = ring.blocksize;
866        //System.out.println("#blocks = " + (matrix.size()/na) + ", na = " + na 
867        //    + " SeqMultBlockTrans");
868        ArrayList<ArrayList<C>> m = matrix;
869        //ArrayList<ArrayList<C>> s = S.matrix;
870
871        GenMatrixRing<C> tr = S.ring.transpose();
872        GenMatrix<C> T = S.transpose(tr);
873        ArrayList<ArrayList<C>> t = T.matrix;
874        //System.out.println("T = " + T); 
875
876        GenMatrixRing<C> pr = ring.product(S.ring);
877        GenMatrix<C> P = pr.getZERO().copy();
878        ArrayList<ArrayList<C>> p = P.matrix;
879        //System.out.println("P = " + P); 
880
881        for (int ii = 0; ii < m.size(); ii += na) {
882            for (int jj = 0; jj < t.size(); jj += nb) {
883
884                for (int i = ii; i < Math.min((ii + na), m.size()); i++) {
885                    ArrayList<C> Ai = m.get(i); //A[i];
886                    for (int j = jj; j < Math.min((jj + nb), t.size()); j++) {
887                        ArrayList<C> Bj = t.get(j); //B[j];
888                        C c = ring.coFac.getZERO();
889                        for (int k = 0; k < Bj.size(); k++) {
890                            c = c.sum(Ai.get(k).multiply(Bj.get(k)));
891                            //  c += Ai[k] * Bj[k];
892                        }
893                        (p.get(i)).set(j, c); // C[i][j] = c;
894                    }
895                }
896
897            }
898        }
899        return new GenMatrix<C>(pr, p);
900    }
901
902
903    /**
904     * Multiply this with S. Simple unblocked algorithm.
905     * @param S other matrix.
906     * @return this * S.
907     */
908    public GenMatrix<C> multiplySimple(GenMatrix<C> S) {
909        ArrayList<ArrayList<C>> m = matrix;
910        ArrayList<ArrayList<C>> B = S.matrix;
911
912        GenMatrixRing<C> pr = ring.product(S.ring);
913        GenMatrix<C> P = pr.getZERO().copy();
914        ArrayList<ArrayList<C>> p = P.matrix;
915
916        for (int i = 0; i < pr.rows; i++) {
917            ArrayList<C> Ai = m.get(i); //A[i];
918            for (int j = 0; j < pr.cols; j++) {
919                C c = ring.coFac.getZERO();
920                for (int k = 0; k < S.ring.rows; k++) {
921                    c = c.sum(Ai.get(k).multiply(B.get(k).get(j)));
922                    //  c += A[i][k] * B[k][j];
923                }
924                (p.get(i)).set(j, c); // C[i][j] = c;
925            }
926        }
927        return new GenMatrix<C>(pr, p);
928    }
929
930
931    /**
932     * Divide this by S.
933     * @param S other matrix.
934     * @return this * S^{-1}.
935     */
936    public GenMatrix<C> divide(GenMatrix<C> S) {
937        //throw new UnsupportedOperationException("divide not yet implemented");
938        return multiply(S.inverse());
939    }
940
941
942    /**
943     * Divide left this by S.
944     * @param S other matrix.
945     * @return S^{-1} * this.
946     */
947    public GenMatrix<C> divideLeft(GenMatrix<C> S) {
948        return S.inverse().multiply(this);
949    }
950
951
952    /**
953     * Remainder after division of this by S.
954     * @param S other matrix.
955     * @return this - (this / S) * S.
956     */
957    public GenMatrix<C> remainder(GenMatrix<C> S) {
958        throw new UnsupportedOperationException("remainder not implemented");
959    }
960
961
962    /**
963     * Quotient and remainder by division of this by S.
964     * @param S a GenMatrix
965     * @return [this/S, this - (this/S)*S].
966     */
967    public GenMatrix<C>[] quotientRemainder(GenMatrix<C> S) {
968        throw new UnsupportedOperationException("quotientRemainder not implemented, input = " + S);
969    }
970
971
972    /**
973     * Inverse of this.
974     * @return x with this * x = 1, if it exists.
975     * @see edu.jas.vector.LinAlg#inverseLU(edu.jas.vector.GenMatrix,java.util.List)
976     */
977    public GenMatrix<C> inverse() {
978        //throw new UnsupportedOperationException("inverse implemented in LinAlg.inverseLU()");
979        LinAlg<C> la = new LinAlg<C>();
980        GenMatrix<C> mat = this.copy();
981        List<Integer> P = la.decompositionLU(mat);
982        if (P == null || P.isEmpty()) {
983            throw new NotInvertibleException("matrix not invertible");
984        }
985        mat = la.inverseLU(mat,P);
986        return mat;
987    }
988
989
990    /**
991     * Greatest common divisor.
992     * @param b other element.
993     * @return gcd(this,b).
994     */
995    public GenMatrix<C> gcd(GenMatrix<C> b) {
996        throw new UnsupportedOperationException("gcd not implemented");
997    }
998
999
1000    /**
1001     * Extended greatest common divisor.
1002     * @param b other element.
1003     * @return [ gcd(this,b), c1, c2 ] with c1*this + c2*b = gcd(this,b).
1004     */
1005    public GenMatrix<C>[] egcd(GenMatrix<C> b) {
1006        throw new UnsupportedOperationException("egcd not implemented");
1007    }
1008
1009}