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