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