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