View Javadoc

1   /*
2    * Copyright 2003-2004 The Apache Software Foundation.
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *      http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
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()) {   // note: this has side effect of attempting LU decomp if lu == null
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 { // LU decomp must have been successfully performed
538             return false; // so the matrix is not singular
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()) { // side effect: compute LU decomp
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         // Apply permutations to b
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         // Solve LY = b
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         // Solve UX = Y
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         // Initialize permutation array and parity
727         permutation = new int[nRows];
728         for (int row = 0; row < nRows; row++) {
729             permutation[row] = row;
730         }
731         parity = 1;
732 
733         // Loop over columns
734         for (int col = 0; col < nCols; col++) {
735 
736             double sum = 0;
737 
738             // upper
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             // lower
748             int max = col; // permutation row
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                 // maintain best permutation choice
758                 if (Math.abs(sum) > largest) {
759                     largest = Math.abs(sum);
760                     max = row;
761                 }
762             }
763 
764             // Singularity check
765             if (Math.abs(lu[max][col]) < TOO_SMALL) {
766                 lu = null;
767                 throw new InvalidMatrixException("matrix is singular");
768             }
769 
770             // Pivot if necessary
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             //Divide the lower elements by the "winning" diagonal elt.
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     //------------------------ Protected methods
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     //------------------------ Private methods
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         // can't copy 2-d array in one shot, otherwise get row references
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 }