View Javadoc

1   /*
2    * Copyright 2003-2005 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 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()) {   // note: this has side effect of attempting LU decomp if lu == null
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 { // LU decomp must have been successfully performed
594             return false; // so the matrix is not singular
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()) { // side effect: compute LU decomp
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         // Apply permutations to b
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         // Solve LY = b
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         // Solve UX = Y
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         // Initialize permutation array and parity
783         permutation = new int[nRows];
784         for (int row = 0; row < nRows; row++) {
785             permutation[row] = row;
786         }
787         parity = 1;
788 
789         // Loop over columns
790         for (int col = 0; col < nCols; col++) {
791 
792             double sum = 0;
793 
794             // upper
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             // lower
804             int max = col; // permutation row
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                 // maintain best permutation choice
814                 if (Math.abs(sum) > largest) {
815                     largest = Math.abs(sum);
816                     max = row;
817                 }
818             }
819 
820             // Singularity check
821             if (Math.abs(lu[max][col]) < TOO_SMALL) {
822                 lu = null;
823                 throw new InvalidMatrixException("matrix is singular");
824             }
825 
826             // Pivot if necessary
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             //Divide the lower elements by the "winning" diagonal elt.
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     //------------------------ Protected methods
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     //------------------------ Private methods
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         // can't copy 2-d array in one shot, otherwise get row references
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 }