001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021/* Derived from the MathWorks and NIST implementation, which was released into
022 * the public domain.
023 *
024 * http://math.nist.gov/javanumerics/jama/
025 */
026package org.biojava.nbio.structure.jama;
027
028import java.io.BufferedReader;
029import java.io.PrintWriter;
030import java.io.StreamTokenizer;
031import java.io.StringWriter;
032import java.text.DecimalFormat;
033import java.text.DecimalFormatSymbols;
034import java.text.NumberFormat;
035import java.util.Locale;
036
037
038/**
039 *      Jama = Java Matrix class.
040 * <P>
041 *      The Java Matrix Class provides the fundamental operations of numerical
042 *      linear algebra.  Various constructors create Matrices from two dimensional
043 *      arrays of double precision floating point numbers.  Various "gets" and
044 *      "sets" provide access to submatrices and matrix elements.  Several methods
045 *      implement basic matrix arithmetic, including matrix addition and
046 *      multiplication, matrix norms, and element-by-element array operations.
047 *      Methods for reading and printing matrices are also included.  All the
048 *      operations in this version of the Matrix Class involve real matrices.
049 *      Complex matrices may be handled in a future version.
050 * <P>
051 *      Five fundamental matrix decompositions, which consist of pairs or triples
052 *      of matrices, permutation vectors, and the like, produce results in five
053 *      decomposition classes.  These decompositions are accessed by the Matrix
054 *      class to compute solutions of simultaneous linear equations, determinants,
055 *      inverses and other matrix functions.  The five decompositions are:
056 * <P><UL>
057 *      <LI>Cholesky Decomposition of symmetric, positive definite matrices.
058 *      <LI>LU Decomposition of rectangular matrices.
059 *      <LI>QR Decomposition of rectangular matrices.
060 *      <LI>Singular Value Decomposition of rectangular matrices.
061 *      <LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
062 * </UL>
063 * <DL>
064 * <DT><B>Example of use:</B></DT>
065 * <P>
066 * <DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||.
067 * <P><PRE>
068 *              double[][] vals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
069 *              Matrix A = new Matrix(vals);
070 *              Matrix b = Matrix.random(3,1);
071 *              Matrix x = A.solve(b);
072 *              Matrix r = A.times(x).minus(b);
073 *              double rnorm = r.normInf();
074 * </PRE></DD>
075 * </DL>
076 *
077 * @author The MathWorks, Inc. and the National Institute of Standards and Technology.
078 * @version 5 August 1998
079 */
080
081public class Matrix implements Cloneable, java.io.Serializable {
082
083         static final long serialVersionUID = 8492558293015348719l;
084
085
086/* ------------------------
087        Class variables
088 * ------------------------ */
089
090        /** Array for internal storage of elements.
091        @serial internal array storage.
092        */
093        private double[][] A;
094
095        /** Row and column dimensions.
096        @serial row dimension.
097        @serial column dimension.
098        */
099        private int m, n;
100
101/* ------------------------
102        Constructors
103 * ------------------------ */
104
105        /** Construct an m-by-n matrix of zeros.
106        @param m    Number of rows.
107        @param n    Number of colums.
108        */
109
110        public Matrix (int m, int n) {
111                this.m = m;
112                this.n = n;
113                A = new double[m][n];
114        }
115
116        /** Construct an m-by-n constant matrix.
117        @param m    Number of rows.
118        @param n    Number of colums.
119        @param s    Fill the matrix with this scalar value.
120        */
121
122        public Matrix (int m, int n, double s) {
123                this.m = m;
124                this.n = n;
125                A = new double[m][n];
126                for (int i = 0; i < m; i++) {
127                        for (int j = 0; j < n; j++) {
128                                A[i][j] = s;
129                        }
130                }
131        }
132
133        /** Construct a matrix from a 2-D array.
134        @param A    Two-dimensional array of doubles.
135        @exception  IllegalArgumentException All rows must have the same length
136        @see        #constructWithCopy
137        */
138
139        public Matrix (double[][] A) {
140                m = A.length;
141                n = A[0].length;
142                for (int i = 0; i < m; i++) {
143                        if (A[i].length != n) {
144                                throw new IllegalArgumentException("All rows must have the same length.");
145                        }
146                }
147                this.A = A;
148        }
149
150        /** Construct a matrix quickly without checking arguments.
151        @param A    Two-dimensional array of doubles.
152        @param m    Number of rows.
153        @param n    Number of colums.
154        */
155
156        public Matrix (double[][] A, int m, int n) {
157                this.A = A;
158                this.m = m;
159                this.n = n;
160        }
161
162        /** Construct a matrix from a one-dimensional packed array
163        @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
164        @param m    Number of rows.
165        @exception  IllegalArgumentException Array length must be a multiple of m.
166        */
167
168        public Matrix (double[] vals, int m) {
169                this.m = m;
170                n = (m != 0 ? vals.length/m : 0);
171                if (m*n != vals.length) {
172                        throw new IllegalArgumentException("Array length must be a multiple of m.");
173                }
174                A = new double[m][n];
175                for (int i = 0; i < m; i++) {
176                        for (int j = 0; j < n; j++) {
177                                A[i][j] = vals[i+j*m];
178                        }
179                }
180        }
181
182/* ------------------------
183        Public Methods
184 * ------------------------ */
185
186        /** Construct a matrix from a copy of a 2-D array.
187        @param A    Two-dimensional array of doubles.
188        @exception  IllegalArgumentException All rows must have the same length
189        @return a Matrix
190        */
191
192        public static Matrix constructWithCopy(double[][] A) {
193                int m = A.length;
194                int n = A[0].length;
195                Matrix X = new Matrix(m,n);
196                double[][] C = X.getArray();
197                for (int i = 0; i < m; i++) {
198                        if (A[i].length != n) {
199                                throw new IllegalArgumentException
200                                        ("All rows must have the same length.");
201                        }
202                        for (int j = 0; j < n; j++) {
203                                C[i][j] = A[i][j];
204                        }
205                }
206                return X;
207        }
208
209        /** Make a deep copy of a matrix
210         * @return a identical copy of the Matrix
211        */
212
213        public Matrix copy () {
214                Matrix X = new Matrix(m,n);
215                double[][] C = X.getArray();
216                for (int i = 0; i < m; i++) {
217                        for (int j = 0; j < n; j++) {
218                                C[i][j] = A[i][j];
219                        }
220                }
221                return X;
222        }
223
224        /** Clone the Matrix object.
225        */
226
227        @Override
228public Object clone () {
229                return this.copy();
230        }
231
232        /** Access the internal two-dimensional array.
233        @return     Pointer to the two-dimensional array of matrix elements.
234        */
235
236        public double[][] getArray () {
237                return A;
238        }
239
240        /** Copy the internal two-dimensional array.
241        @return     Two-dimensional array copy of matrix elements.
242        */
243
244        public double[][] getArrayCopy () {
245                double[][] C = new double[m][n];
246                for (int i = 0; i < m; i++) {
247                        for (int j = 0; j < n; j++) {
248                                C[i][j] = A[i][j];
249                        }
250                }
251                return C;
252        }
253
254        /** Make a one-dimensional column packed copy of the internal array.
255        @return     Matrix elements packed in a one-dimensional array by columns.
256        */
257
258        public double[] getColumnPackedCopy () {
259                double[] vals = new double[m*n];
260                for (int i = 0; i < m; i++) {
261                        for (int j = 0; j < n; j++) {
262                                vals[i+j*m] = A[i][j];
263                        }
264                }
265                return vals;
266        }
267
268        /** Make a one-dimensional row packed copy of the internal array.
269        @return     Matrix elements packed in a one-dimensional array by rows.
270        */
271
272        public double[] getRowPackedCopy () {
273                double[] vals = new double[m*n];
274                for (int i = 0; i < m; i++) {
275                        for (int j = 0; j < n; j++) {
276                                vals[i*n+j] = A[i][j];
277                        }
278                }
279                return vals;
280        }
281
282        /** Get row dimension.
283        @return     m, the number of rows.
284        */
285
286        public int getRowDimension () {
287                return m;
288        }
289
290        /** Get column dimension.
291        @return     n, the number of columns.
292        */
293
294        public int getColumnDimension () {
295                return n;
296        }
297
298        /** Get a single element.
299        @param i    Row index.
300        @param j    Column index.
301        @return     A(i,j)
302        @exception  ArrayIndexOutOfBoundsException
303        */
304
305        public double get (int i, int j) {
306                return A[i][j];
307        }
308
309        /** Get a submatrix.
310        @param i0   Initial row index
311        @param i1   Final row index
312        @param j0   Initial column index
313        @param j1   Final column index
314        @return     A(i0:i1,j0:j1)
315        @exception  ArrayIndexOutOfBoundsException Submatrix indices
316        */
317
318        public Matrix getMatrix (int i0, int i1, int j0, int j1) {
319                Matrix X = new Matrix(i1-i0+1,j1-j0+1);
320                double[][] B = X.getArray();
321                try {
322                        for (int i = i0; i <= i1; i++) {
323                                for (int j = j0; j <= j1; j++) {
324                                        B[i-i0][j-j0] = A[i][j];
325                                }
326                        }
327                } catch(ArrayIndexOutOfBoundsException e) {
328                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
329                }
330                return X;
331        }
332
333        /** Get a submatrix.
334        @param r    Array of row indices.
335        @param c    Array of column indices.
336        @return     A(r(:),c(:))
337        @exception  ArrayIndexOutOfBoundsException Submatrix indices
338        */
339
340        public Matrix getMatrix (int[] r, int[] c) {
341                Matrix X = new Matrix(r.length,c.length);
342                double[][] B = X.getArray();
343                try {
344                        for (int i = 0; i < r.length; i++) {
345                                for (int j = 0; j < c.length; j++) {
346                                        B[i][j] = A[r[i]][c[j]];
347                                }
348                        }
349                } catch(ArrayIndexOutOfBoundsException e) {
350                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
351                }
352                return X;
353        }
354
355        /** Get a submatrix.
356        @param i0   Initial row index
357        @param i1   Final row index
358        @param c    Array of column indices.
359        @return     A(i0:i1,c(:))
360        @exception  ArrayIndexOutOfBoundsException Submatrix indices
361        */
362
363        public Matrix getMatrix (int i0, int i1, int[] c) {
364                Matrix X = new Matrix(i1-i0+1,c.length);
365                double[][] B = X.getArray();
366                try {
367                        for (int i = i0; i <= i1; i++) {
368                                for (int j = 0; j < c.length; j++) {
369                                        B[i-i0][j] = A[i][c[j]];
370                                }
371                        }
372                } catch(ArrayIndexOutOfBoundsException e) {
373                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
374                }
375                return X;
376        }
377
378        /** Get a submatrix.
379        @param r    Array of row indices.
380        @param j0   Initial column index
381        @param j1   Final column index
382        @return     A(r(:),j0:j1)
383        @exception  ArrayIndexOutOfBoundsException Submatrix indices
384        */
385
386        public Matrix getMatrix (int[] r, int j0, int j1) {
387                Matrix X = new Matrix(r.length,j1-j0+1);
388                double[][] B = X.getArray();
389                try {
390                        for (int i = 0; i < r.length; i++) {
391                                for (int j = j0; j <= j1; j++) {
392                                        B[i][j-j0] = A[r[i]][j];
393                                }
394                        }
395                } catch(ArrayIndexOutOfBoundsException e) {
396                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
397                }
398                return X;
399        }
400
401        /** Set a single element.
402        @param i    Row index.
403        @param j    Column index.
404        @param s    A(i,j).
405        @exception  ArrayIndexOutOfBoundsException
406        */
407
408        public void set (int i, int j, double s) {
409                A[i][j] = s;
410        }
411
412        /** Set a submatrix.
413        @param i0   Initial row index
414        @param i1   Final row index
415        @param j0   Initial column index
416        @param j1   Final column index
417        @param X    A(i0:i1,j0:j1)
418        @exception  ArrayIndexOutOfBoundsException Submatrix indices
419        */
420
421        public void setMatrix (int i0, int i1, int j0, int j1, Matrix X) {
422                try {
423                        for (int i = i0; i <= i1; i++) {
424                                for (int j = j0; j <= j1; j++) {
425                                        A[i][j] = X.get(i-i0,j-j0);
426                                }
427                        }
428                } catch(ArrayIndexOutOfBoundsException e) {
429                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
430                }
431        }
432
433        /** Set a submatrix.
434        @param r    Array of row indices.
435        @param c    Array of column indices.
436        @param X    A(r(:),c(:))
437        @exception  ArrayIndexOutOfBoundsException Submatrix indices
438        */
439
440        public void setMatrix (int[] r, int[] c, Matrix X) {
441                try {
442                        for (int i = 0; i < r.length; i++) {
443                                for (int j = 0; j < c.length; j++) {
444                                        A[r[i]][c[j]] = X.get(i,j);
445                                }
446                        }
447                } catch(ArrayIndexOutOfBoundsException e) {
448                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
449                }
450        }
451
452        /** Set a submatrix.
453        @param r    Array of row indices.
454        @param j0   Initial column index
455        @param j1   Final column index
456        @param X    A(r(:),j0:j1)
457        @exception  ArrayIndexOutOfBoundsException Submatrix indices
458        */
459
460        public void setMatrix (int[] r, int j0, int j1, Matrix X) {
461                try {
462                        for (int i = 0; i < r.length; i++) {
463                                for (int j = j0; j <= j1; j++) {
464                                        A[r[i]][j] = X.get(i,j-j0);
465                                }
466                        }
467                } catch(ArrayIndexOutOfBoundsException e) {
468                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
469                }
470        }
471
472        /** Set a submatrix.
473        @param i0   Initial row index
474        @param i1   Final row index
475        @param c    Array of column indices.
476        @param X    A(i0:i1,c(:))
477        @exception  ArrayIndexOutOfBoundsException Submatrix indices
478        */
479
480        public void setMatrix (int i0, int i1, int[] c, Matrix X) {
481                try {
482                        for (int i = i0; i <= i1; i++) {
483                                for (int j = 0; j < c.length; j++) {
484                                        A[i][c[j]] = X.get(i-i0,j);
485                                }
486                        }
487                } catch(ArrayIndexOutOfBoundsException e) {
488                        throw new ArrayIndexOutOfBoundsException("Submatrix indices");
489                }
490        }
491
492        /** Matrix transpose.
493        @return    A'
494        */
495
496        public Matrix transpose () {
497                Matrix X = new Matrix(n,m);
498                double[][] C = X.getArray();
499                for (int i = 0; i < m; i++) {
500                        for (int j = 0; j < n; j++) {
501                                C[j][i] = A[i][j];
502                        }
503                }
504                return X;
505        }
506
507        /** One norm
508        @return    maximum column sum.
509        */
510
511        public double norm1 () {
512                double f = 0;
513                for (int j = 0; j < n; j++) {
514                        double s = 0;
515                        for (int i = 0; i < m; i++) {
516                                s += Math.abs(A[i][j]);
517                        }
518                        f = Math.max(f,s);
519                }
520                return f;
521        }
522
523        /** Two norm
524        @return    maximum singular value.
525        */
526
527        public double norm2 () {
528                return (new SingularValueDecomposition(this).norm2());
529        }
530
531        /** Infinity norm
532        @return    maximum row sum.
533        */
534
535        public double normInf () {
536                double f = 0;
537                for (int i = 0; i < m; i++) {
538                        double s = 0;
539                        for (int j = 0; j < n; j++) {
540                                s += Math.abs(A[i][j]);
541                        }
542                        f = Math.max(f,s);
543                }
544                return f;
545        }
546
547        /** Frobenius norm
548        @return    sqrt of sum of squares of all elements.
549        */
550
551        public double normF () {
552                double f = 0;
553                for (int i = 0; i < m; i++) {
554                        for (int j = 0; j < n; j++) {
555                                f = Maths.hypot(f,A[i][j]);
556                        }
557                }
558                return f;
559        }
560
561        /**  Unary minus
562        @return    -A
563        */
564
565        public Matrix uminus () {
566                Matrix X = new Matrix(m,n);
567                double[][] C = X.getArray();
568                for (int i = 0; i < m; i++) {
569                        for (int j = 0; j < n; j++) {
570                                C[i][j] = -A[i][j];
571                        }
572                }
573                return X;
574        }
575
576        /** C = A + B
577        @param B    another matrix
578        @return     A + B
579        */
580
581        public Matrix plus (Matrix B) {
582                checkMatrixDimensions(B);
583                Matrix X = new Matrix(m,n);
584                double[][] C = X.getArray();
585                for (int i = 0; i < m; i++) {
586                        for (int j = 0; j < n; j++) {
587                                C[i][j] = A[i][j] + B.A[i][j];
588                        }
589                }
590                return X;
591        }
592
593        /** A = A + B
594        @param B    another matrix
595        @return     A + B
596        */
597
598        public Matrix plusEquals (Matrix B) {
599                checkMatrixDimensions(B);
600                for (int i = 0; i < m; i++) {
601                        for (int j = 0; j < n; j++) {
602                                A[i][j] = A[i][j] + B.A[i][j];
603                        }
604                }
605                return this;
606        }
607
608        /** C = A - B
609        @param B    another matrix
610        @return     A - B
611        */
612
613        public Matrix minus (Matrix B) {
614                checkMatrixDimensions(B);
615                Matrix X = new Matrix(m,n);
616                double[][] C = X.getArray();
617                for (int i = 0; i < m; i++) {
618                        for (int j = 0; j < n; j++) {
619                                C[i][j] = A[i][j] - B.A[i][j];
620                        }
621                }
622                return X;
623        }
624
625        /** A = A - B
626        @param B    another matrix
627        @return     A - B
628        */
629
630        public Matrix minusEquals (Matrix B) {
631                checkMatrixDimensions(B);
632                for (int i = 0; i < m; i++) {
633                        for (int j = 0; j < n; j++) {
634                                A[i][j] = A[i][j] - B.A[i][j];
635                        }
636                }
637                return this;
638        }
639
640        /** Element-by-element multiplication, C = A.*B
641        @param B    another matrix
642        @return     A.*B
643        */
644
645        public Matrix arrayTimes (Matrix B) {
646                checkMatrixDimensions(B);
647                Matrix X = new Matrix(m,n);
648                double[][] C = X.getArray();
649                for (int i = 0; i < m; i++) {
650                        for (int j = 0; j < n; j++) {
651                                C[i][j] = A[i][j] * B.A[i][j];
652                        }
653                }
654                return X;
655        }
656
657        /** Element-by-element multiplication in place, A = A.*B
658        @param B    another matrix
659        @return     A.*B
660        */
661
662        public Matrix arrayTimesEquals (Matrix B) {
663                checkMatrixDimensions(B);
664                for (int i = 0; i < m; i++) {
665                        for (int j = 0; j < n; j++) {
666                                A[i][j] = A[i][j] * B.A[i][j];
667                        }
668                }
669                return this;
670        }
671
672        /** Element-by-element right division, C = A./B
673        @param B    another matrix
674        @return     A./B
675        */
676
677        public Matrix arrayRightDivide (Matrix B) {
678                checkMatrixDimensions(B);
679                Matrix X = new Matrix(m,n);
680                double[][] C = X.getArray();
681                for (int i = 0; i < m; i++) {
682                        for (int j = 0; j < n; j++) {
683                                C[i][j] = A[i][j] / B.A[i][j];
684                        }
685                }
686                return X;
687        }
688
689        /** Element-by-element right division in place, A = A./B
690        @param B    another matrix
691        @return     A./B
692        */
693
694        public Matrix arrayRightDivideEquals (Matrix B) {
695                checkMatrixDimensions(B);
696                for (int i = 0; i < m; i++) {
697                        for (int j = 0; j < n; j++) {
698                                A[i][j] = A[i][j] / B.A[i][j];
699                        }
700                }
701                return this;
702        }
703
704        /** Element-by-element left division, C = A.\B
705        @param B    another matrix
706        @return     A.\B
707        */
708
709        public Matrix arrayLeftDivide (Matrix B) {
710                checkMatrixDimensions(B);
711                Matrix X = new Matrix(m,n);
712                double[][] C = X.getArray();
713                for (int i = 0; i < m; i++) {
714                        for (int j = 0; j < n; j++) {
715                                C[i][j] = B.A[i][j] / A[i][j];
716                        }
717                }
718                return X;
719        }
720
721        /** Element-by-element left division in place, A = A.\B
722        @param B    another matrix
723        @return     A.\B
724        */
725
726        public Matrix arrayLeftDivideEquals (Matrix B) {
727                checkMatrixDimensions(B);
728                for (int i = 0; i < m; i++) {
729                        for (int j = 0; j < n; j++) {
730                                A[i][j] = B.A[i][j] / A[i][j];
731                        }
732                }
733                return this;
734        }
735
736        /** Multiply a matrix by a scalar, C = s*A
737        @param s    scalar
738        @return     s*A
739        */
740
741        public Matrix times (double s) {
742                Matrix X = new Matrix(m,n);
743                double[][] C = X.getArray();
744                for (int i = 0; i < m; i++) {
745                        for (int j = 0; j < n; j++) {
746                                C[i][j] = s*A[i][j];
747                        }
748                }
749                return X;
750        }
751
752        /** Multiply a matrix by a scalar in place, A = s*A
753        @param s    scalar
754        @return     replace A by s*A
755        */
756
757        public Matrix timesEquals (double s) {
758                for (int i = 0; i < m; i++) {
759                        for (int j = 0; j < n; j++) {
760                                A[i][j] = s*A[i][j];
761                        }
762                }
763                return this;
764        }
765
766        /** Linear algebraic matrix multiplication, A * B
767        @param B    another matrix
768        @return     Matrix product, A * B
769        @exception  IllegalArgumentException Matrix inner dimensions must agree.
770        */
771
772        public Matrix times (Matrix B) {
773                if (B.m != n) {
774                        throw new IllegalArgumentException("Matrix inner dimensions must agree.");
775                }
776                Matrix X = new Matrix(m,B.n);
777                double[][] C = X.getArray();
778                double[] Bcolj = new double[n];
779                for (int j = 0; j < B.n; j++) {
780                        for (int k = 0; k < n; k++) {
781                                Bcolj[k] = B.A[k][j];
782                        }
783                        for (int i = 0; i < m; i++) {
784                                double[] Arowi = A[i];
785                                double s = 0;
786                                for (int k = 0; k < n; k++) {
787                                        s += Arowi[k]*Bcolj[k];
788                                }
789                                C[i][j] = s;
790                        }
791                }
792                return X;
793        }
794
795        /** LU Decomposition
796        @return     LUDecomposition
797        @see LUDecomposition
798        */
799
800        public LUDecomposition lu () {
801                return new LUDecomposition(this);
802        }
803
804        /** QR Decomposition
805        @return     QRDecomposition
806        @see QRDecomposition
807        */
808
809        public QRDecomposition qr () {
810                return new QRDecomposition(this);
811        }
812
813        /** Cholesky Decomposition
814        @return     CholeskyDecomposition
815        @see CholeskyDecomposition
816        */
817
818        public CholeskyDecomposition chol () {
819                return new CholeskyDecomposition(this);
820        }
821
822        /** Singular Value Decomposition
823        @return     SingularValueDecomposition
824        @see SingularValueDecomposition
825        */
826
827        public SingularValueDecomposition svd () {
828                return new SingularValueDecomposition(this);
829        }
830
831        /** Eigenvalue Decomposition
832        @return     EigenvalueDecomposition
833        @see EigenvalueDecomposition
834        */
835
836        public EigenvalueDecomposition eig () {
837                return new EigenvalueDecomposition(this);
838        }
839
840        /** Solve A*X = B
841        @param B    right hand side
842        @return     solution if A is square, least squares solution otherwise
843        */
844
845        public Matrix solve (Matrix B) {
846                return (m == n ? (new LUDecomposition(this)).solve(B) :
847                                                          (new QRDecomposition(this)).solve(B));
848        }
849
850        /** Solve X*A = B, which is also A'*X' = B'
851        @param B    right hand side
852        @return     solution if A is square, least squares solution otherwise.
853        */
854
855        public Matrix solveTranspose (Matrix B) {
856                return transpose().solve(B.transpose());
857        }
858
859        /** Matrix inverse or pseudoinverse
860        @return     inverse(A) if A is square, pseudoinverse otherwise.
861        */
862
863        public Matrix inverse () {
864                return solve(identity(m,m));
865        }
866
867        /** Matrix determinant
868        @return     determinant
869        */
870
871        public double det () {
872                return new LUDecomposition(this).det();
873        }
874
875        /** Matrix rank
876        @return     effective numerical rank, obtained from SVD.
877        */
878
879        public int rank () {
880                return new SingularValueDecomposition(this).rank();
881        }
882
883        /** Matrix condition (2 norm)
884        @return     ratio of largest to smallest singular value.
885        */
886
887        public double cond () {
888                return new SingularValueDecomposition(this).cond();
889        }
890
891        /** Matrix trace.
892        @return     sum of the diagonal elements.
893        */
894
895        public double trace () {
896                double t = 0;
897                for (int i = 0; i < Math.min(m,n); i++) {
898                        t += A[i][i];
899                }
900                return t;
901        }
902
903        /** Generate matrix with random elements
904        @param m    Number of rows.
905        @param n    Number of colums.
906        @return     An m-by-n matrix with uniformly distributed random elements.
907        */
908
909        public static Matrix random (int m, int n) {
910                Matrix A = new Matrix(m,n);
911                double[][] X = A.getArray();
912                for (int i = 0; i < m; i++) {
913                        for (int j = 0; j < n; j++) {
914                                X[i][j] = Math.random();
915                        }
916                }
917                return A;
918        }
919
920        /** Generate identity matrix
921        @param m    Number of rows.
922        @param n    Number of colums.
923        @return     An m-by-n matrix with ones on the diagonal and zeros elsewhere.
924        */
925
926        public static Matrix identity (int m, int n) {
927                Matrix A = new Matrix(m,n);
928                double[][] X = A.getArray();
929                for (int i = 0; i < m; i++) {
930                        for (int j = 0; j < n; j++) {
931                                X[i][j] = (i == j ? 1.0 : 0.0);
932                        }
933                }
934                return A;
935        }
936
937        @Override
938public String toString(){
939                StringWriter writer = new StringWriter();
940                PrintWriter printWriter = new PrintWriter(writer);
941                print(printWriter,getColumnDimension(),3);
942                return writer.toString();
943        }
944
945
946        /** Print the matrix to stdout.   Line the elements up in columns
947          * with a Fortran-like 'Fw.d' style format.
948        @param w    Column width.
949        @param d    Number of digits after the decimal.
950        */
951
952        public void print (int w, int d) {
953                print(new PrintWriter(System.out,true),w,d); }
954
955        /** Print the matrix to the output stream.   Line the elements up in
956          * columns with a Fortran-like 'Fw.d' style format.
957        @param output Output stream.
958        @param w      Column width.
959        @param d      Number of digits after the decimal.
960        */
961
962        public void print (PrintWriter output, int w, int d) {
963                DecimalFormat format = new DecimalFormat();
964                format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
965                format.setMinimumIntegerDigits(1);
966                format.setMaximumFractionDigits(d);
967                format.setMinimumFractionDigits(d);
968                format.setGroupingUsed(false);
969                print(output,format,w+2);
970        }
971
972        /** Print the matrix to stdout.  Line the elements up in columns.
973          * Use the format object, and right justify within columns of width
974          * characters.
975          * Note that is the matrix is to be read back in, you probably will want
976          * to use a NumberFormat that is set to US Locale.
977        @param format A  Formatting object for individual elements.
978        @param width     Field width for each column.
979        @see java.text.DecimalFormat#setDecimalFormatSymbols
980        */
981
982        public void print (NumberFormat format, int width) {
983                print(new PrintWriter(System.out,true),format,width); }
984
985        // DecimalFormat is a little disappointing coming from Fortran or C's printf.
986        // Since it doesn't pad on the left, the elements will come out different
987        // widths.  Consequently, we'll pass the desired column width in as an
988        // argument and do the extra padding ourselves.
989
990        /** Print the matrix to the output stream.  Line the elements up in columns.
991          * Use the format object, and right justify within columns of width
992          * characters.
993          * Note that is the matrix is to be read back in, you probably will want
994          * to use a NumberFormat that is set to US Locale.
995        @param output the output stream.
996        @param format A formatting object to format the matrix elements
997        @param width  Column width.
998        @see java.text.DecimalFormat#setDecimalFormatSymbols
999        */
1000
1001        public void print (PrintWriter output, NumberFormat format, int width) {
1002                output.println();  // start on new line.
1003                for (int i = 0; i < m; i++) {
1004                        for (int j = 0; j < n; j++) {
1005                                String s = format.format(A[i][j]); // format the number
1006                                int padding = Math.max(1,width-s.length()); // At _least_ 1 space
1007                                for (int k = 0; k < padding; k++)
1008                                        output.print(' ');
1009                                output.print(s);
1010                        }
1011                        output.println();
1012                }
1013                //output.println();   // end with blank line.
1014        }
1015
1016        /** Read a matrix from a stream.  The format is the same the print method,
1017          * so printed matrices can be read back in (provided they were printed using
1018          * US Locale).  Elements are separated by
1019          * whitespace, all the elements for each row appear on a single line,
1020          * the last row is followed by a blank line.
1021        @param input the input stream.
1022        @return a Matrix
1023        @throws java.io.IOException
1024        */
1025
1026        @SuppressWarnings({ "rawtypes", "unchecked" })
1027public static Matrix read (BufferedReader input) throws java.io.IOException {
1028                StreamTokenizer tokenizer= new StreamTokenizer(input);
1029
1030                // Although StreamTokenizer will parse numbers, it doesn't recognize
1031                // scientific notation (E or D); however, Double.valueOf does.
1032                // The strategy here is to disable StreamTokenizer's number parsing.
1033                // We'll only get whitespace delimited words, EOL's and EOF's.
1034                // These words should all be numbers, for Double.valueOf to parse.
1035
1036                tokenizer.resetSyntax();
1037                tokenizer.wordChars(0,255);
1038                tokenizer.whitespaceChars(0, ' ');
1039                tokenizer.eolIsSignificant(true);
1040                java.util.Vector v = new java.util.Vector();
1041
1042                // Ignore initial empty lines
1043                while (tokenizer.nextToken() == StreamTokenizer.TT_EOL);
1044                if (tokenizer.ttype == StreamTokenizer.TT_EOF)
1045        throw new java.io.IOException("Unexpected EOF on matrix read.");
1046                do {
1047                        v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
1048                } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
1049
1050                int n = v.size();  // Now we've got the number of columns!
1051                double[] row = new double[n];
1052                for (int j=0; j<n; j++)  // extract the elements of the 1st row.
1053                        row[j]=((Double)v.elementAt(j)).doubleValue();
1054                v.removeAllElements();
1055                v.addElement(row);  // Start storing rows instead of columns.
1056                while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
1057                        // While non-empty lines
1058                        v.addElement(row = new double[n]);
1059                        int j = 0;
1060                        do {
1061                                if (j >= n) throw new java.io.IOException
1062                                        ("Row " + v.size() + " is too long.");
1063                                row[j++] = Double.parseDouble(tokenizer.sval);
1064                        } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
1065                        if (j < n) throw new java.io.IOException
1066                                ("Row " + v.size() + " is too short.");
1067                }
1068                int m = v.size();  // Now we've got the number of rows.
1069                double[][] A = new double[m][];
1070                v.copyInto(A);  // copy the rows out of the vector
1071                return new Matrix(A);
1072        }
1073
1074
1075/* ------------------------
1076        Private Methods
1077 * ------------------------ */
1078
1079        /** Check if size(A) == size(B) **/
1080
1081        private void checkMatrixDimensions (Matrix B) {
1082                if (B.m != m || B.n != n) {
1083                        throw new IllegalArgumentException("Matrix dimensions must agree.");
1084                }
1085        }
1086
1087}