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