1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.linear;
18 import java.io.Serializable;
19 import java.math.BigDecimal;
20
21 /***
22 * Implementation for {@link BigMatrix} using a BigDecimal[][] array to store entries
23 * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
24 * LU decompostion</a> to support linear system
25 * solution and inverse.
26 * <p>
27 * The LU decompostion is performed as needed, to support the following operations: <ul>
28 * <li>solve</li>
29 * <li>isSingular</li>
30 * <li>getDeterminant</li>
31 * <li>inverse</li> </ul>
32 * <p>
33 * <strong>Usage notes</strong>:<br>
34 * <ul><li>
35 * The LU decomposition is stored and reused on subsequent calls. If matrix
36 * data are modified using any of the public setXxx methods, the saved
37 * decomposition is discarded. If data are modified via references to the
38 * underlying array obtained using <code>getDataRef()</code>, then the stored
39 * LU decomposition will not be discarded. In this case, you need to
40 * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
41 * before using any of the methods above.</li>
42 * <li>
43 * As specified in the {@link BigMatrix} interface, matrix element indexing
44 * is 0-based -- e.g., <code>getEntry(0, 0)</code>
45 * returns the element in the first row, first column of the matrix.</li></ul>
46 * @version $Revision: 1.10 $ $Date: 2004/11/07 20:19:22 $
47 */
48 public class BigMatrixImpl implements BigMatrix, Serializable {
49
50 /*** Serialization id */
51 static final long serialVersionUID = -1011428905656140431L;
52
53 /*** The number zero. */
54 private static final BigDecimal ZERO = new BigDecimal(0);
55
56 /*** The number one. */
57 private static final BigDecimal ONE = new BigDecimal(1);
58
59 /*** Entries of the matrix */
60 private BigDecimal data[][] = null;
61
62 /*** Entries of cached LU decomposition.
63 * All updates to data (other than luDecompose()) *must* set this to null
64 */
65 private BigDecimal lu[][] = null;
66
67 /*** Permutation associated with LU decompostion */
68 private int[] permutation = null;
69
70 /*** Parity of the permutation associated with the LU decomposition */
71 private int parity = 1;
72
73 /*** Rounding mode for divisions **/
74 private int roundingMode = BigDecimal.ROUND_HALF_UP;
75
76 /**** BigDecimal scale ***/
77 private int scale = 64;
78
79 /*** Bound to determine effective singularity in LU decomposition */
80 protected static BigDecimal TOO_SMALL = new BigDecimal(10E-12);
81
82 /***
83 * Creates a matrix with no data
84 */
85 public BigMatrixImpl() {
86 }
87
88 /***
89 * Create a new BigMatrix with the supplied row and column dimensions.
90 *
91 * @param rowDimension the number of rows in the new matrix
92 * @param columnDimension the number of columns in the new matrix
93 */
94 public BigMatrixImpl(int rowDimension, int columnDimension) {
95 data = new BigDecimal[rowDimension][columnDimension];
96 lu = null;
97 }
98
99 /***
100 * Create a new BigMatrix using the <code>data</code> as the underlying
101 * data array.
102 * <p>
103 * The input array is copied, not referenced.
104 *
105 * @param d data for new matrix
106 * @throws IllegalArgumentException if <code>d</code> is not rectangular
107 * (not all rows have the same length) or empty
108 * @throws NullPointerException if <code>d</code> is null
109 */
110 public BigMatrixImpl(BigDecimal[][] d) {
111 int nRows = d.length;
112 if (nRows == 0) {
113 throw new IllegalArgumentException(
114 "Matrix must have at least one row.");
115 }
116 int nCols = d[0].length;
117 if (nCols == 0) {
118 throw new IllegalArgumentException(
119 "Matrix must have at least one column.");
120 }
121 for (int row = 1; row < nRows; row++) {
122 if (d[row].length != nCols) {
123 throw new IllegalArgumentException(
124 "All input rows must have the same length.");
125 }
126 }
127 this.copyIn(d);
128 lu = null;
129 }
130
131 /***
132 * Create a new BigMatrix using the <code>data</code> as the underlying
133 * data array.
134 * <p>
135 * The input array is copied, not referenced.
136 *
137 * @param d data for new matrix
138 * @throws IllegalArgumentException if <code>d</code> is not rectangular
139 * (not all rows have the same length) or empty
140 * @throws NullPointerException if <code>d</code> is null
141 */
142 public BigMatrixImpl(double[][] d) {
143 int nRows = d.length;
144 if (nRows == 0) {
145 throw new IllegalArgumentException(
146 "Matrix must have at least one row.");
147 }
148 int nCols = d[0].length;
149 if (nCols == 0) {
150 throw new IllegalArgumentException(
151 "Matrix must have at least one column.");
152 }
153 for (int row = 1; row < nRows; row++) {
154 if (d[row].length != nCols) {
155 throw new IllegalArgumentException(
156 "All input rows must have the same length.");
157 }
158 }
159 this.copyIn(d);
160 lu = null;
161 }
162
163 /***
164 * Create a new BigMatrix using the values represented by the strings in
165 * <code>data</code> as the underlying data array.
166 *
167 * @param d data for new matrix
168 * @throws IllegalArgumentException if <code>d</code> is not rectangular
169 * (not all rows have the same length) or empty
170 * @throws NullPointerException if <code>d</code> is null
171 */
172 public BigMatrixImpl(String[][] d) {
173 int nRows = d.length;
174 if (nRows == 0) {
175 throw new IllegalArgumentException(
176 "Matrix must have at least one row.");
177 }
178 int nCols = d[0].length;
179 if (nCols == 0) {
180 throw new IllegalArgumentException(
181 "Matrix must have at least one column.");
182 }
183 for (int row = 1; row < nRows; row++) {
184 if (d[row].length != nCols) {
185 throw new IllegalArgumentException(
186 "All input rows must have the same length.");
187 }
188 }
189 this.copyIn(d);
190 lu = null;
191 }
192
193 /***
194 * Create a new (column) BigMatrix using <code>v</code> as the
195 * data for the unique column of the <code>v.length x 1</code> matrix
196 * created.
197 * <p>
198 * The input array is copied, not referenced.
199 *
200 * @param v column vector holding data for new matrix
201 */
202 public BigMatrixImpl(BigDecimal[] v) {
203 int nRows = v.length;
204 data = new BigDecimal[nRows][1];
205 for (int row = 0; row < nRows; row++) {
206 data[row][0] = v[row];
207 }
208 }
209
210 /***
211 * Create a new BigMatrix which is a copy of this.
212 *
213 * @return the cloned matrix
214 */
215 public BigMatrix copy() {
216 return new BigMatrixImpl(this.copyOut());
217 }
218
219 /***
220 * Compute the sum of this and <code>m</code>.
221 *
222 * @param m matrix to be added
223 * @return this + m
224 * @exception IllegalArgumentException if m is not the same size as this
225 */
226 public BigMatrix add(BigMatrix m) throws IllegalArgumentException {
227 if (this.getColumnDimension() != m.getColumnDimension() ||
228 this.getRowDimension() != m.getRowDimension()) {
229 throw new IllegalArgumentException("matrix dimension mismatch");
230 }
231 int rowCount = this.getRowDimension();
232 int columnCount = this.getColumnDimension();
233 BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
234 for (int row = 0; row < rowCount; row++) {
235 for (int col = 0; col < columnCount; col++) {
236 outData[row][col] = data[row][col].add(m.getEntry(row, col));
237 }
238 }
239 return new BigMatrixImpl(outData);
240 }
241
242 /***
243 * Compute this minus <code>m</code>.
244 *
245 * @param m matrix to be subtracted
246 * @return this + m
247 * @exception IllegalArgumentException if m is not the same size as *this
248 */
249 public BigMatrix subtract(BigMatrix m) throws IllegalArgumentException {
250 if (this.getColumnDimension() != m.getColumnDimension() ||
251 this.getRowDimension() != m.getRowDimension()) {
252 throw new IllegalArgumentException("matrix dimension mismatch");
253 }
254 int rowCount = this.getRowDimension();
255 int columnCount = this.getColumnDimension();
256 BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
257 for (int row = 0; row < rowCount; row++) {
258 for (int col = 0; col < columnCount; col++) {
259 outData[row][col] = data[row][col].subtract(m.getEntry(row, col));
260 }
261 }
262 return new BigMatrixImpl(outData);
263 }
264
265 /***
266 * Returns the result of adding d to each entry of this.
267 *
268 * @param d value to be added to each entry
269 * @return d + this
270 */
271 public BigMatrix scalarAdd(BigDecimal d) {
272 int rowCount = this.getRowDimension();
273 int columnCount = this.getColumnDimension();
274 BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
275 for (int row = 0; row < rowCount; row++) {
276 for (int col = 0; col < columnCount; col++) {
277 outData[row][col] = data[row][col].add(d);
278 }
279 }
280 return new BigMatrixImpl(outData);
281 }
282
283 /***
284 * Returns the result multiplying each entry of this by <code>d</code>
285 * @param d value to multiply all entries by
286 * @return d * this
287 */
288 public BigMatrix scalarMultiply(BigDecimal d) {
289 int rowCount = this.getRowDimension();
290 int columnCount = this.getColumnDimension();
291 BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
292 for (int row = 0; row < rowCount; row++) {
293 for (int col = 0; col < columnCount; col++) {
294 outData[row][col] = data[row][col].multiply(d);
295 }
296 }
297 return new BigMatrixImpl(outData);
298 }
299
300 /***
301 * Returns the result of postmultiplying this by <code>m</code>.
302 * @param m matrix to postmultiply by
303 * @return this*m
304 * @throws IllegalArgumentException
305 * if columnDimension(this) != rowDimension(m)
306 */
307 public BigMatrix multiply(BigMatrix m) throws IllegalArgumentException {
308 if (this.getColumnDimension() != m.getRowDimension()) {
309 throw new IllegalArgumentException("Matrices are not multiplication compatible.");
310 }
311 int nRows = this.getRowDimension();
312 int nCols = m.getColumnDimension();
313 int nSum = this.getColumnDimension();
314 BigDecimal[][] outData = new BigDecimal[nRows][nCols];
315 BigDecimal sum = ZERO;
316 for (int row = 0; row < nRows; row++) {
317 for (int col = 0; col < nCols; col++) {
318 sum = ZERO;
319 for (int i = 0; i < nSum; i++) {
320 sum = sum.add(data[row][i].multiply(m.getEntry(i, col)));
321 }
322 outData[row][col] = sum;
323 }
324 }
325 return new BigMatrixImpl(outData);
326 }
327
328 /***
329 * Returns the result premultiplying this by <code>m</code>.
330 * @param m matrix to premultiply by
331 * @return m * this
332 * @throws IllegalArgumentException
333 * if rowDimension(this) != columnDimension(m)
334 */
335 public BigMatrix preMultiply(BigMatrix m) throws IllegalArgumentException {
336 return m.multiply(this);
337 }
338
339 /***
340 * Returns matrix entries as a two-dimensional array.
341 * <p>
342 * Makes a fresh copy of the underlying data.
343 *
344 * @return 2-dimensional array of entries
345 */
346 public BigDecimal[][] getData() {
347 return copyOut();
348 }
349
350 /***
351 * Returns matrix entries as a two-dimensional array.
352 * <p>
353 * Makes a fresh copy of the underlying data converted to
354 * <code>double</code> values.
355 *
356 * @return 2-dimensional array of entries
357 */
358 public double[][] getDataAsDoubleArray() {
359 int nRows = getRowDimension();
360 int nCols = getColumnDimension();
361 double d[][] = new double[nRows][nCols];
362 for (int i = 0; i < nRows; i++) {
363 for (int j=0; j<nCols;j++) {
364 d[i][j] = data[i][j].doubleValue();
365 }
366 }
367 return d;
368 }
369
370 /***
371 * Returns a reference to the underlying data array.
372 * <p>
373 * Does not make a fresh copy of the underlying data.
374 *
375 * @return 2-dimensional array of entries
376 */
377 public BigDecimal[][] getDataRef() {
378 return data;
379 }
380
381 /****
382 * Gets the rounding mode for division operations
383 * The default is {@link java.math.BigDecimal#ROUND_HALF_UP}
384 * @see BigDecimal
385 * @return the rounding mode.
386 */
387 public int getRoundingMode() {
388 return roundingMode;
389 }
390
391 /****
392 * Sets the rounding mode for decimal divisions.
393 * @see BigDecimal
394 * @param roundingMode
395 */
396 public void setRoundingMode(int roundingMode) {
397 this.roundingMode = roundingMode;
398 }
399
400 /****
401 * Sets the scale for division operations.
402 * The default is 64
403 * @see BigDecimal
404 * @return the scale
405 */
406 public int getScale() {
407 return scale;
408 }
409
410 /****
411 * Sets the scale for division operations.
412 * @see BigDecimal
413 * @param scale
414 */
415 public void setScale(int scale) {
416 this.scale = scale;
417 }
418
419 /***
420 * Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html">
421 * maximum absolute row sum norm</a> of the matrix.
422 *
423 * @return norm
424 */
425 public BigDecimal getNorm() {
426 BigDecimal maxColSum = ZERO;
427 for (int col = 0; col < this.getColumnDimension(); col++) {
428 BigDecimal sum = ZERO;
429 for (int row = 0; row < this.getRowDimension(); row++) {
430 sum = sum.add(data[row][col].abs());
431 }
432 maxColSum = maxColSum.max(sum);
433 }
434 return maxColSum;
435 }
436
437 /***
438 * Gets a submatrix. Rows and columns are indicated
439 * counting from 0 to n-1.
440 *
441 * @param startRow Initial row index
442 * @param endRow Final row index
443 * @param startColumn Initial column index
444 * @param endColumn Final column index
445 * @return The subMatrix containing the data of the
446 * specified rows and columns
447 * @exception MatrixIndexException if row or column selections are not valid
448 */
449 public BigMatrix getSubMatrix(int startRow, int endRow, int startColumn,
450 int endColumn) throws MatrixIndexException {
451 if (startRow < 0 || startRow > endRow || endRow > data.length ||
452 startColumn < 0 || startColumn > endColumn ||
453 endColumn > data[0].length ) {
454 throw new MatrixIndexException(
455 "invalid row or column index selection");
456 }
457 BigMatrixImpl subMatrix = new BigMatrixImpl(endRow - startRow+1,
458 endColumn - startColumn+1);
459 BigDecimal[][] subMatrixData = subMatrix.getDataRef();
460 for (int i = startRow; i <= endRow; i++) {
461 for (int j = startColumn; j <= endColumn; j++) {
462 subMatrixData[i - startRow][j - startColumn] = data[i][j];
463 }
464 }
465 return subMatrix;
466 }
467
468 /***
469 * Gets a submatrix. Rows and columns are indicated
470 * counting from 0 to n-1.
471 *
472 * @param selectedRows Array of row indices must be non-empty
473 * @param selectedColumns Array of column indices must be non-empty
474 * @return The subMatrix containing the data in the
475 * specified rows and columns
476 * @exception MatrixIndexException if supplied row or column index arrays
477 * are not valid
478 */
479 public BigMatrix getSubMatrix(int[] selectedRows, int[] selectedColumns)
480 throws MatrixIndexException {
481 if (selectedRows.length * selectedColumns.length == 0) {
482 throw new MatrixIndexException(
483 "selected row and column index arrays must be non-empty");
484 }
485 BigMatrixImpl subMatrix = new BigMatrixImpl(selectedRows.length,
486 selectedColumns.length);
487 BigDecimal[][] subMatrixData = subMatrix.getDataRef();
488 try {
489 for (int i = 0; i < selectedRows.length; i++) {
490 for (int j = 0; j < selectedColumns.length; j++) {
491 subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
492 }
493 }
494 }
495 catch (ArrayIndexOutOfBoundsException e) {
496 throw new MatrixIndexException("matrix dimension mismatch");
497 }
498 return subMatrix;
499 }
500
501 /***
502 * Returns the entries in row number <code>row</code>
503 * as a row matrix. Row indices start at 0.
504 *
505 * @param row the row to be fetched
506 * @return row matrix
507 * @throws MatrixIndexException if the specified row index is invalid
508 */
509 public BigMatrix getRowMatrix(int row) throws MatrixIndexException {
510 if ( !isValidCoordinate( row, 0)) {
511 throw new MatrixIndexException("illegal row argument");
512 }
513 int ncols = this.getColumnDimension();
514 BigDecimal[][] out = new BigDecimal[1][ncols];
515 System.arraycopy(data[row], 0, out[0], 0, ncols);
516 return new BigMatrixImpl(out);
517 }
518
519 /***
520 * Returns the entries in column number <code>column</code>
521 * as a column matrix. Column indices start at 0.
522 *
523 * @param column the column to be fetched
524 * @return column matrix
525 * @throws MatrixIndexException if the specified column index is invalid
526 */
527 public BigMatrix getColumnMatrix(int column) throws MatrixIndexException {
528 if ( !isValidCoordinate( 0, column)) {
529 throw new MatrixIndexException("illegal column argument");
530 }
531 int nRows = this.getRowDimension();
532 BigDecimal[][] out = new BigDecimal[nRows][1];
533 for (int row = 0; row < nRows; row++) {
534 out[row][0] = data[row][column];
535 }
536 return new BigMatrixImpl(out);
537 }
538
539 /***
540 * Returns the entries in row number <code>row</code> as an array.
541 * <p>
542 * Row indices start at 0. A <code>MatrixIndexException</code> is thrown
543 * unless <code>0 <= row < rowDimension.</code>
544 *
545 * @param row the row to be fetched
546 * @return array of entries in the row
547 * @throws MatrixIndexException if the specified row index is not valid
548 */
549 public BigDecimal[] getRow(int row) throws MatrixIndexException {
550 if ( !isValidCoordinate( row, 0 ) ) {
551 throw new MatrixIndexException("illegal row argument");
552 }
553 int ncols = this.getColumnDimension();
554 BigDecimal[] out = new BigDecimal[ncols];
555 System.arraycopy(data[row], 0, out, 0, ncols);
556 return out;
557 }
558
559 /***
560 * Returns the entries in row number <code>row</code> as an array
561 * of double values.
562 * <p>
563 * Row indices start at 0. A <code>MatrixIndexException</code> is thrown
564 * unless <code>0 <= row < rowDimension.</code>
565 *
566 * @param row the row to be fetched
567 * @return array of entries in the row
568 * @throws MatrixIndexException if the specified row index is not valid
569 */
570 public double[] getRowAsDoubleArray(int row) throws MatrixIndexException {
571 if ( !isValidCoordinate( row, 0 ) ) {
572 throw new MatrixIndexException("illegal row argument");
573 }
574 int ncols = this.getColumnDimension();
575 double[] out = new double[ncols];
576 for (int i=0;i<ncols;i++) {
577 out[i] = data[row][i].doubleValue();
578 }
579 return out;
580 }
581
582 /***
583 * Returns the entries in column number <code>col</code> as an array.
584 * <p>
585 * Column indices start at 0. A <code>MatrixIndexException</code> is thrown
586 * unless <code>0 <= column < columnDimension.</code>
587 *
588 * @param col the column to be fetched
589 * @return array of entries in the column
590 * @throws MatrixIndexException if the specified column index is not valid
591 */
592 public BigDecimal[] getColumn(int col) throws MatrixIndexException {
593 if ( !isValidCoordinate(0, col) ) {
594 throw new MatrixIndexException("illegal column argument");
595 }
596 int nRows = this.getRowDimension();
597 BigDecimal[] out = new BigDecimal[nRows];
598 for (int i = 0; i < nRows; i++) {
599 out[i] = data[i][col];
600 }
601 return out;
602 }
603
604 /***
605 * Returns the entries in column number <code>col</code> as an array
606 * of double values.
607 * <p>
608 * Column indices start at 0. A <code>MatrixIndexException</code> is thrown
609 * unless <code>0 <= column < columnDimension.</code>
610 *
611 * @param col the column to be fetched
612 * @return array of entries in the column
613 * @throws MatrixIndexException if the specified column index is not valid
614 */
615 public double[] getColumnAsDoubleArray(int col) throws MatrixIndexException {
616 if ( !isValidCoordinate( 0, col ) ) {
617 throw new MatrixIndexException("illegal column argument");
618 }
619 int nrows = this.getRowDimension();
620 double[] out = new double[nrows];
621 for (int i=0;i<nrows;i++) {
622 out[i] = data[i][col].doubleValue();
623 }
624 return out;
625 }
626
627 /***
628 * Returns the entry in the specified row and column.
629 * <p>
630 * Row and column indices start at 0 and must satisfy
631 * <ul>
632 * <li><code>0 <= row < rowDimension</code></li>
633 * <li><code> 0 <= column < columnDimension</code></li>
634 * </ul>
635 * otherwise a <code>MatrixIndexException</code> is thrown.
636 *
637 * @param row row location of entry to be fetched
638 * @param column column location of entry to be fetched
639 * @return matrix entry in row,column
640 * @throws MatrixIndexException if the row or column index is not valid
641 */
642 public BigDecimal getEntry(int row, int column)
643 throws MatrixIndexException {
644 if (!isValidCoordinate(row,column)) {
645 throw new MatrixIndexException("matrix entry does not exist");
646 }
647 return data[row][column];
648 }
649
650 /***
651 * Returns the entry in the specified row and column as a double.
652 * <p>
653 * Row and column indices start at 0 and must satisfy
654 * <ul>
655 * <li><code>0 <= row < rowDimension</code></li>
656 * <li><code> 0 <= column < columnDimension</code></li>
657 * </ul>
658 * otherwise a <code>MatrixIndexException</code> is thrown.
659 *
660 * @param row row location of entry to be fetched
661 * @param column column location of entry to be fetched
662 * @return matrix entry in row,column
663 * @throws MatrixIndexException if the row
664 * or column index is not valid
665 */
666 public double getEntryAsDouble(int row, int column) throws MatrixIndexException {
667 return getEntry(row,column).doubleValue();
668 }
669
670 /***
671 * Returns the transpose matrix.
672 *
673 * @return transpose matrix
674 */
675 public BigMatrix transpose() {
676 int nRows = this.getRowDimension();
677 int nCols = this.getColumnDimension();
678 BigMatrixImpl out = new BigMatrixImpl(nCols, nRows);
679 BigDecimal[][] outData = out.getDataRef();
680 for (int row = 0; row < nRows; row++) {
681 for (int col = 0; col < nCols; col++) {
682 outData[col][row] = data[row][col];
683 }
684 }
685 return out;
686 }
687
688 /***
689 * Returns the inverse matrix if this matrix is invertible.
690 *
691 * @return inverse matrix
692 * @throws InvalidMatrixException if this is not invertible
693 */
694 public BigMatrix inverse() throws InvalidMatrixException {
695 return solve(getIdentity(this.getRowDimension()));
696 }
697
698 /***
699 * Returns the determinant of this matrix.
700 *
701 * @return determinant
702 * @throws InvalidMatrixException if matrix is not square
703 */
704 public BigDecimal getDeterminant() throws InvalidMatrixException {
705 if (!isSquare()) {
706 throw new InvalidMatrixException("matrix is not square");
707 }
708 if (isSingular()) {
709 return ZERO;
710 } else {
711 BigDecimal det = (parity == 1) ? ONE : ONE.negate();
712 for (int i = 0; i < this.getRowDimension(); i++) {
713 det = det.multiply(lu[i][i]);
714 }
715 return det;
716 }
717 }
718
719 /***
720 * Is this a square matrix?
721 * @return true if the matrix is square (rowDimension = columnDimension)
722 */
723 public boolean isSquare() {
724 return (this.getColumnDimension() == this.getRowDimension());
725 }
726
727 /***
728 * Is this a singular matrix?
729 * @return true if the matrix is singular
730 */
731 public boolean isSingular() {
732 if (lu == null) {
733 try {
734 luDecompose();
735 return false;
736 } catch (InvalidMatrixException ex) {
737 return true;
738 }
739 } else {
740 return false;
741 }
742 }
743
744 /***
745 * Returns the number of rows in the matrix.
746 *
747 * @return rowDimension
748 */
749 public int getRowDimension() {
750 return data.length;
751 }
752
753 /***
754 * Returns the number of columns in the matrix.
755 *
756 * @return columnDimension
757 */
758 public int getColumnDimension() {
759 return data[0].length;
760 }
761
762 /***
763 * Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
764 * trace</a> of the matrix (the sum of the elements on the main diagonal).
765 *
766 * @return trace
767 *
768 * @throws IllegalArgumentException if this matrix is not square.
769 */
770 public BigDecimal getTrace() throws IllegalArgumentException {
771 if (!isSquare()) {
772 throw new IllegalArgumentException("matrix is not square");
773 }
774 BigDecimal trace = data[0][0];
775 for (int i = 1; i < this.getRowDimension(); i++) {
776 trace = trace.add(data[i][i]);
777 }
778 return trace;
779 }
780
781 /***
782 * Returns the result of multiplying this by the vector <code>v</code>.
783 *
784 * @param v the vector to operate on
785 * @return this*v
786 * @throws IllegalArgumentException if columnDimension != v.size()
787 */
788 public BigDecimal[] operate(BigDecimal[] v) throws IllegalArgumentException {
789 if (v.length != this.getColumnDimension()) {
790 throw new IllegalArgumentException("vector has wrong length");
791 }
792 int nRows = this.getRowDimension();
793 int nCols = this.getColumnDimension();
794 BigDecimal[] out = new BigDecimal[v.length];
795 for (int row = 0; row < nRows; row++) {
796 BigDecimal sum = ZERO;
797 for (int i = 0; i < nCols; i++) {
798 sum = sum.add(data[row][i].multiply(v[i]));
799 }
800 out[row] = sum;
801 }
802 return out;
803 }
804
805 /***
806 * Returns the result of multiplying this by the vector <code>v</code>.
807 *
808 * @param v the vector to operate on
809 * @return this*v
810 * @throws IllegalArgumentException if columnDimension != v.size()
811 */
812 public BigDecimal[] operate(double[] v) throws IllegalArgumentException {
813 BigDecimal bd[] = new BigDecimal[v.length];
814 for (int i=0;i<bd.length;i++) {
815 bd[i] = new BigDecimal(v[i]);
816 }
817 return operate(bd);
818 }
819
820 /***
821 * Returns the (row) vector result of premultiplying this by the vector <code>v</code>.
822 *
823 * @param v the row vector to premultiply by
824 * @return v*this
825 * @throws IllegalArgumentException if rowDimension != v.size()
826 */
827 public BigDecimal[] preMultiply(BigDecimal[] v) throws IllegalArgumentException {
828 int nRows = this.getRowDimension();
829 if (v.length != nRows) {
830 throw new IllegalArgumentException("vector has wrong length");
831 }
832 int nCols = this.getColumnDimension();
833 BigDecimal[] out = new BigDecimal[nCols];
834 for (int col = 0; col < nCols; col++) {
835 BigDecimal sum = ZERO;
836 for (int i = 0; i < nRows; i++) {
837 sum = sum.add(data[i][col].multiply(v[i]));
838 }
839 out[col] = sum;
840 }
841 return out;
842 }
843
844 /***
845 * Returns a matrix of (column) solution vectors for linear systems with
846 * coefficient matrix = this and constant vectors = columns of
847 * <code>b</code>.
848 *
849 * @param b array of constants forming RHS of linear systems to
850 * to solve
851 * @return solution array
852 * @throws IllegalArgumentException if this.rowDimension != row dimension
853 * @throws InvalidMatrixException if this matrix is not square or is singular
854 */
855 public BigDecimal[] solve(BigDecimal[] b) throws IllegalArgumentException, InvalidMatrixException {
856 int nRows = this.getRowDimension();
857 if (b.length != nRows) {
858 throw new IllegalArgumentException("constant vector has wrong length");
859 }
860 BigMatrix bMatrix = new BigMatrixImpl(b);
861 BigDecimal[][] solution = ((BigMatrixImpl) (solve(bMatrix))).getDataRef();
862 BigDecimal[] out = new BigDecimal[nRows];
863 for (int row = 0; row < nRows; row++) {
864 out[row] = solution[row][0];
865 }
866 return out;
867 }
868
869 /***
870 * Returns a matrix of (column) solution vectors for linear systems with
871 * coefficient matrix = this and constant vectors = columns of
872 * <code>b</code>.
873 *
874 * @param b array of constants forming RHS of linear systems to
875 * to solve
876 * @return solution array
877 * @throws IllegalArgumentException if this.rowDimension != row dimension
878 * @throws InvalidMatrixException if this matrix is not square or is singular
879 */
880 public BigDecimal[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
881 BigDecimal bd[] = new BigDecimal[b.length];
882 for (int i=0;i<bd.length;i++) {
883 bd[i] = new BigDecimal(b[i]);
884 }
885 return solve(bd);
886 }
887
888 /***
889 * Returns a matrix of (column) solution vectors for linear systems with
890 * coefficient matrix = this and constant vectors = columns of
891 * <code>b</code>.
892 *
893 * @param b matrix of constant vectors forming RHS of linear systems to
894 * to solve
895 * @return matrix of solution vectors
896 * @throws IllegalArgumentException if this.rowDimension != row dimension
897 * @throws InvalidMatrixException if this matrix is not square or is singular
898 */
899 public BigMatrix solve(BigMatrix b) throws IllegalArgumentException, InvalidMatrixException {
900 if (b.getRowDimension() != this.getRowDimension()) {
901 throw new IllegalArgumentException("Incorrect row dimension");
902 }
903 if (!this.isSquare()) {
904 throw new InvalidMatrixException("coefficient matrix is not square");
905 }
906 if (this.isSingular()) {
907 throw new InvalidMatrixException("Matrix is singular.");
908 }
909
910 int nCol = this.getColumnDimension();
911 int nColB = b.getColumnDimension();
912 int nRowB = b.getRowDimension();
913
914
915 BigDecimal[][] bp = new BigDecimal[nRowB][nColB];
916 for (int row = 0; row < nRowB; row++) {
917 for (int col = 0; col < nColB; col++) {
918 bp[row][col] = b.getEntry(permutation[row], col);
919 }
920 }
921
922
923 for (int col = 0; col < nCol; col++) {
924 for (int i = col + 1; i < nCol; i++) {
925 for (int j = 0; j < nColB; j++) {
926 bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
927 }
928 }
929 }
930
931
932 for (int col = nCol - 1; col >= 0; col--) {
933 for (int j = 0; j < nColB; j++) {
934 bp[col][j] = bp[col][j].divide(lu[col][col], scale, roundingMode);
935 }
936 for (int i = 0; i < col; i++) {
937 for (int j = 0; j < nColB; j++) {
938 bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
939 }
940 }
941 }
942
943 BigMatrixImpl outMat = new BigMatrixImpl(bp);
944 return outMat;
945 }
946
947 /***
948 * Computes a new
949 * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
950 * LU decompostion</a> for this matrix, storing the result for use by other methods.
951 * <p>
952 * <strong>Implementation Note</strong>:<br>
953 * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
954 * Crout's algortithm</a>, with partial pivoting.
955 * <p>
956 * <strong>Usage Note</strong>:<br>
957 * This method should rarely be invoked directly. Its only use is
958 * to force recomputation of the LU decomposition when changes have been
959 * made to the underlying data using direct array references. Changes
960 * made using setXxx methods will trigger recomputation when needed
961 * automatically.
962 *
963 * @throws InvalidMatrixException if the matrix is non-square or singular.
964 */
965 public void luDecompose() throws InvalidMatrixException {
966
967 int nRows = this.getRowDimension();
968 int nCols = this.getColumnDimension();
969 if (nRows != nCols) {
970 throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
971 }
972 lu = this.getData();
973
974
975 permutation = new int[nRows];
976 for (int row = 0; row < nRows; row++) {
977 permutation[row] = row;
978 }
979 parity = 1;
980
981
982 for (int col = 0; col < nCols; col++) {
983
984 BigDecimal sum = ZERO;
985
986
987 for (int row = 0; row < col; row++) {
988 sum = lu[row][col];
989 for (int i = 0; i < row; i++) {
990 sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
991 }
992 lu[row][col] = sum;
993 }
994
995
996 int max = col;
997 BigDecimal largest = ZERO;
998 for (int row = col; row < nRows; row++) {
999 sum = lu[row][col];
1000 for (int i = 0; i < col; i++) {
1001 sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
1002 }
1003 lu[row][col] = sum;
1004
1005
1006 if (sum.abs().compareTo(largest) == 1) {
1007 largest = sum.abs();
1008 max = row;
1009 }
1010 }
1011
1012
1013 if (lu[max][col].abs().compareTo(TOO_SMALL) <= 0) {
1014 lu = null;
1015 throw new InvalidMatrixException("matrix is singular");
1016 }
1017
1018
1019 if (max != col) {
1020 BigDecimal tmp = ZERO;
1021 for (int i = 0; i < nCols; i++) {
1022 tmp = lu[max][i];
1023 lu[max][i] = lu[col][i];
1024 lu[col][i] = tmp;
1025 }
1026 int temp = permutation[max];
1027 permutation[max] = permutation[col];
1028 permutation[col] = temp;
1029 parity = -parity;
1030 }
1031
1032
1033 for (int row = col + 1; row < nRows; row++) {
1034 lu[row][col] = lu[row][col].divide(lu[col][col], scale, roundingMode);
1035 }
1036
1037 }
1038
1039 }
1040
1041 /***
1042 *
1043 * @see Object#toString()
1044 */
1045 public String toString() {
1046 StringBuffer res = new StringBuffer();
1047 res.append("BigMatrixImpl{");
1048 if (data != null) {
1049 for (int i = 0; i < data.length; i++) {
1050 if (i > 0)
1051 res.append(",");
1052 res.append("{");
1053 for (int j = 0; j < data[0].length; j++) {
1054 if (j > 0)
1055 res.append(",");
1056 res.append(data[i][j]);
1057 }
1058 res.append("}");
1059 }
1060 }
1061 res.append("}");
1062 return res.toString();
1063 }
1064
1065 /***
1066 * Returns true iff <code>object</code> is a
1067 * <code>BigMatrixImpl</code> instance with the same dimensions as this
1068 * and all corresponding matrix entries are equal. BigDecimal.equals
1069 * is used to compare corresponding entries.
1070 *
1071 * @param object the object to test equality against.
1072 * @return true if object equals this
1073 */
1074 public boolean equals(Object object) {
1075 if (object == this ) {
1076 return true;
1077 }
1078 if (object instanceof BigMatrixImpl == false) {
1079 return false;
1080 }
1081 BigMatrix m = (BigMatrix) object;
1082 int nRows = getRowDimension();
1083 int nCols = getColumnDimension();
1084 if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
1085 return false;
1086 }
1087 for (int row = 0; row < nRows; row++) {
1088 for (int col = 0; col < nCols; col++) {
1089 if (!data[row][col].equals(m.getEntry(row, col))) {
1090 return false;
1091 }
1092 }
1093 }
1094 return true;
1095 }
1096
1097 /***
1098 * Computes a hashcode for the matrix.
1099 *
1100 * @return hashcode for matrix
1101 */
1102 public int hashCode() {
1103 int ret = 7;
1104 int nRows = getRowDimension();
1105 int nCols = getColumnDimension();
1106 ret = ret * 31 + nRows;
1107 ret = ret * 31 + nCols;
1108 for (int row = 0; row < nRows; row++) {
1109 for (int col = 0; col < nCols; col++) {
1110 ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) *
1111 data[row][col].hashCode();
1112 }
1113 }
1114 return ret;
1115 }
1116
1117
1118
1119 /***
1120 * Returns <code>dimension x dimension</code> identity matrix.
1121 *
1122 * @param dimension dimension of identity matrix to generate
1123 * @return identity matrix
1124 */
1125 protected BigMatrix getIdentity(int dimension) {
1126 BigMatrixImpl out = new BigMatrixImpl(dimension, dimension);
1127 BigDecimal[][] d = out.getDataRef();
1128 for (int row = 0; row < dimension; row++) {
1129 for (int col = 0; col < dimension; col++) {
1130 d[row][col] = row == col ? ONE : ZERO;
1131 }
1132 }
1133 return out;
1134 }
1135
1136 /***
1137 * Returns the LU decomposition as a BigMatrix.
1138 * Returns a fresh copy of the cached LU matrix if this has been computed;
1139 * otherwise the composition is computed and cached for use by other methods.
1140 * Since a copy is returned in either case, changes to the returned matrix do not
1141 * affect the LU decomposition property.
1142 * <p>
1143 * The matrix returned is a compact representation of the LU decomposition.
1144 * Elements below the main diagonal correspond to entries of the "L" matrix;
1145 * elements on and above the main diagonal correspond to entries of the "U"
1146 * matrix.
1147 * <p>
1148 * Example: <pre>
1149 *
1150 * Returned matrix L U
1151 * 2 3 1 1 0 0 2 3 1
1152 * 5 4 6 5 1 0 0 4 6
1153 * 1 7 8 1 7 1 0 0 8
1154 * </pre>
1155 *
1156 * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
1157 * where permuteRows reorders the rows of the matrix to follow the order determined
1158 * by the <a href=#getPermutation()>permutation</a> property.
1159 *
1160 * @return LU decomposition matrix
1161 * @throws InvalidMatrixException if the matrix is non-square or singular.
1162 */
1163 protected BigMatrix getLUMatrix() throws InvalidMatrixException {
1164 if (lu == null) {
1165 luDecompose();
1166 }
1167 return new BigMatrixImpl(lu);
1168 }
1169
1170 /***
1171 * Returns the permutation associated with the lu decomposition.
1172 * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
1173 * <p>
1174 * Example:
1175 * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
1176 * and current first row is last.
1177 * <p>
1178 * Returns a fresh copy of the array.
1179 *
1180 * @return the permutation
1181 */
1182 protected int[] getPermutation() {
1183 int[] out = new int[permutation.length];
1184 System.arraycopy(permutation, 0, out, 0, permutation.length);
1185 return out;
1186 }
1187
1188
1189
1190 /***
1191 * Returns a fresh copy of the underlying data array.
1192 *
1193 * @return a copy of the underlying data array.
1194 */
1195 private BigDecimal[][] copyOut() {
1196 int nRows = this.getRowDimension();
1197 BigDecimal[][] out = new BigDecimal[nRows][this.getColumnDimension()];
1198
1199 for (int i = 0; i < nRows; i++) {
1200 System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1201 }
1202 return out;
1203 }
1204
1205 /***
1206 * Replaces data with a fresh copy of the input array.
1207 *
1208 * @param in data to copy in
1209 */
1210 private void copyIn(BigDecimal[][] in) {
1211 int nRows = in.length;
1212 int nCols = in[0].length;
1213 data = new BigDecimal[nRows][nCols];
1214 System.arraycopy(in, 0, data, 0, in.length);
1215 for (int i = 0; i < nRows; i++) {
1216 System.arraycopy(in[i], 0, data[i], 0, nCols);
1217 }
1218 lu = null;
1219 }
1220
1221 /***
1222 * Replaces data with a fresh copy of the input array.
1223 *
1224 * @param in data to copy in
1225 */
1226 private void copyIn(double[][] in) {
1227 int nRows = in.length;
1228 int nCols = in[0].length;
1229 data = new BigDecimal[nRows][nCols];
1230 for (int i = 0; i < nRows; i++) {
1231 for (int j=0; j < nCols; j++) {
1232 data[i][j] = new BigDecimal(in[i][j]);
1233 }
1234 }
1235 lu = null;
1236 }
1237
1238 /***
1239 * Replaces data with BigDecimals represented by the strings in the input
1240 * array.
1241 *
1242 * @param in data to copy in
1243 */
1244 private void copyIn(String[][] in) {
1245 int nRows = in.length;
1246 int nCols = in[0].length;
1247 data = new BigDecimal[nRows][nCols];
1248 for (int i = 0; i < nRows; i++) {
1249 for (int j=0; j < nCols; j++) {
1250 data[i][j] = new BigDecimal(in[i][j]);
1251 }
1252 }
1253 lu = null;
1254 }
1255
1256 /***
1257 * Tests a given coordinate as being valid or invalid
1258 *
1259 * @param row the row index.
1260 * @param col the column index.
1261 * @return true if the coordinate is with the current dimensions
1262 */
1263 private boolean isValidCoordinate(int row, int col) {
1264 int nRows = this.getRowDimension();
1265 int nCols = this.getColumnDimension();
1266
1267 return !(row < 0 || row >= nRows || col < 0 || col >= nCols);
1268 }
1269
1270 }