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