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