View Javadoc

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