View Javadoc

1   /*
2    * Copyright 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  import java.io.Serializable;
19  import java.math.BigDecimal;
20  
21  /***
22   * Implementation for {@link BigMatrix} using a BigDecimal[][] array to store entries
23   * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
24   * LU decompostion</a> to support linear system 
25   * solution and inverse.
26   * <p>
27   * The LU decompostion is performed as needed, to support the following operations: <ul>
28   * <li>solve</li>
29   * <li>isSingular</li>
30   * <li>getDeterminant</li>
31   * <li>inverse</li> </ul>
32   * <p>
33  * <strong>Usage notes</strong>:<br>
34   * <ul><li>
35   * The LU decomposition is stored and reused on subsequent calls.  If matrix
36   * data are modified using any of the public setXxx methods, the saved
37   * decomposition is discarded.  If data are modified via references to the
38   * underlying array obtained using <code>getDataRef()</code>, then the stored
39   * LU decomposition will not be discarded.  In this case, you need to
40   * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
41   * before using any of the methods above.</li>
42   * <li>
43   * As specified in the {@link BigMatrix} interface, matrix element indexing
44   * is 0-based -- e.g., <code>getEntry(0, 0)</code>
45   * returns the element in the first row, first column of the matrix.</li></ul>
46   * @version $Revision: 1.10 $ $Date: 2004/11/07 20:19:22 $
47   */
48  public class BigMatrixImpl implements BigMatrix, Serializable {
49      
50      /*** Serialization id */
51      static final long serialVersionUID = -1011428905656140431L;
52      
53      /*** The number zero. */
54      private static final BigDecimal ZERO = new BigDecimal(0);
55      
56      /*** The number one. */
57      private static final BigDecimal ONE = new BigDecimal(1);
58      
59      /*** Entries of the matrix */
60      private BigDecimal data[][] = null;
61      
62      /*** Entries of cached LU decomposition.
63       *  All updates to data (other than luDecompose()) *must* set this to null
64       */
65      private BigDecimal lu[][] = null;
66      
67      /*** Permutation associated with LU decompostion */
68      private int[] permutation = null;
69      
70      /*** Parity of the permutation associated with the LU decomposition */
71      private int parity = 1;
72      
73      /*** Rounding mode for divisions **/
74      private int roundingMode = BigDecimal.ROUND_HALF_UP;
75      
76      /**** BigDecimal scale ***/
77      private int scale = 64;
78      
79      /*** Bound to determine effective singularity in LU decomposition */
80      protected static BigDecimal TOO_SMALL = new BigDecimal(10E-12);
81      
82      /*** 
83       * Creates a matrix with no data
84       */
85      public BigMatrixImpl() {
86      }
87      
88      /***
89       * Create a new BigMatrix with the supplied row and column dimensions.
90       *
91       * @param rowDimension      the number of rows in the new matrix
92       * @param columnDimension   the number of columns in the new matrix
93       */
94      public BigMatrixImpl(int rowDimension, int columnDimension) {
95          data = new BigDecimal[rowDimension][columnDimension];
96          lu = null;
97      }
98      
99      /***
100      * Create a new BigMatrix using the <code>data</code> as the underlying
101      * data array.
102      * <p>
103      * The input array is copied, not referenced.
104      *
105      * @param d data for new matrix
106      * @throws IllegalArgumentException if <code>d</code> is not rectangular
107      *  (not all rows have the same length) or empty
108      * @throws NullPointerException if <code>d</code> is null
109      */
110     public BigMatrixImpl(BigDecimal[][] d) {
111         int nRows = d.length;
112         if (nRows == 0) {
113             throw new IllegalArgumentException(
114             "Matrix must have at least one row."); 
115         }
116         int nCols = d[0].length;
117         if (nCols == 0) {
118             throw new IllegalArgumentException(
119             "Matrix must have at least one column."); 
120         }
121         for (int row = 1; row < nRows; row++) {
122             if (d[row].length != nCols) {
123                 throw new IllegalArgumentException(
124                 "All input rows must have the same length.");
125             }
126         }
127         this.copyIn(d);
128         lu = null;
129     }
130     
131     /***
132      * Create a new BigMatrix using the <code>data</code> as the underlying
133      * data array.
134      * <p>
135      * The input array is copied, not referenced.
136      *
137      * @param d data for new matrix
138      * @throws IllegalArgumentException if <code>d</code> is not rectangular
139      *  (not all rows have the same length) or empty
140      * @throws NullPointerException if <code>d</code> is null
141      */
142     public BigMatrixImpl(double[][] d) {
143         int nRows = d.length;
144         if (nRows == 0) {
145             throw new IllegalArgumentException(
146             "Matrix must have at least one row."); 
147         }
148         int nCols = d[0].length;
149         if (nCols == 0) {
150             throw new IllegalArgumentException(
151             "Matrix must have at least one column."); 
152         }
153         for (int row = 1; row < nRows; row++) {
154             if (d[row].length != nCols) {
155                 throw new IllegalArgumentException(
156                 "All input rows must have the same length.");
157             }
158         }
159         this.copyIn(d);
160         lu = null;
161     }
162     
163     /***
164      * Create a new BigMatrix using the values represented by the strings in 
165      * <code>data</code> as the underlying data array.
166      *
167      * @param d data for new matrix
168      * @throws IllegalArgumentException if <code>d</code> is not rectangular
169      *  (not all rows have the same length) or empty
170      * @throws NullPointerException if <code>d</code> is null
171      */
172     public BigMatrixImpl(String[][] d) {
173         int nRows = d.length;
174         if (nRows == 0) {
175             throw new IllegalArgumentException(
176             "Matrix must have at least one row."); 
177         }
178         int nCols = d[0].length;
179         if (nCols == 0) {
180             throw new IllegalArgumentException(
181             "Matrix must have at least one column."); 
182         }
183         for (int row = 1; row < nRows; row++) {
184             if (d[row].length != nCols) {
185                 throw new IllegalArgumentException(
186                 "All input rows must have the same length.");
187             }
188         }
189         this.copyIn(d);
190         lu = null;
191     }
192     
193     /***
194      * Create a new (column) BigMatrix using <code>v</code> as the
195      * data for the unique column of the <code>v.length x 1</code> matrix 
196      * created.
197      * <p>
198      * The input array is copied, not referenced.
199      *
200      * @param v column vector holding data for new matrix
201      */
202     public BigMatrixImpl(BigDecimal[] v) {
203         int nRows = v.length;
204         data = new BigDecimal[nRows][1];
205         for (int row = 0; row < nRows; row++) {
206             data[row][0] = v[row];
207         }
208     }
209     
210     /***
211      * Create a new BigMatrix which is a copy of this.
212      *
213      * @return  the cloned matrix
214      */
215     public BigMatrix copy() {
216         return new BigMatrixImpl(this.copyOut());
217     }
218     
219     /***
220      * Compute the sum of this and <code>m</code>.
221      *
222      * @param m    matrix to be added
223      * @return     this + m
224      * @exception  IllegalArgumentException if m is not the same size as this
225      */
226     public BigMatrix add(BigMatrix m) throws IllegalArgumentException {
227         if (this.getColumnDimension() != m.getColumnDimension() ||
228                 this.getRowDimension() != m.getRowDimension()) {
229             throw new IllegalArgumentException("matrix dimension mismatch");
230         }
231         int rowCount = this.getRowDimension();
232         int columnCount = this.getColumnDimension();
233         BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
234         for (int row = 0; row < rowCount; row++) {
235             for (int col = 0; col < columnCount; col++) {
236                 outData[row][col] = data[row][col].add(m.getEntry(row, col));
237             }
238         }
239         return new BigMatrixImpl(outData);
240     }
241     
242     /***
243      * Compute  this minus <code>m</code>.
244      *
245      * @param m    matrix to be subtracted
246      * @return     this + m
247      * @exception  IllegalArgumentException if m is not the same size as *this
248      */
249     public BigMatrix subtract(BigMatrix m) throws IllegalArgumentException {
250         if (this.getColumnDimension() != m.getColumnDimension() ||
251                 this.getRowDimension() != m.getRowDimension()) {
252             throw new IllegalArgumentException("matrix dimension mismatch");
253         }
254         int rowCount = this.getRowDimension();
255         int columnCount = this.getColumnDimension();
256         BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
257         for (int row = 0; row < rowCount; row++) {
258             for (int col = 0; col < columnCount; col++) {
259                 outData[row][col] = data[row][col].subtract(m.getEntry(row, col));
260             }
261         }
262         return new BigMatrixImpl(outData);
263     }
264     
265     /***
266      * Returns the result of adding d to each entry of this.
267      *
268      * @param d    value to be added to each entry
269      * @return     d + this
270      */
271     public BigMatrix scalarAdd(BigDecimal d) {
272         int rowCount = this.getRowDimension();
273         int columnCount = this.getColumnDimension();
274         BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
275         for (int row = 0; row < rowCount; row++) {
276             for (int col = 0; col < columnCount; col++) {
277                 outData[row][col] = data[row][col].add(d);
278             }
279         }
280         return new BigMatrixImpl(outData);
281     }
282     
283     /***
284      * Returns the result multiplying each entry of this by <code>d</code>
285      * @param d  value to multiply all entries by
286      * @return d * this
287      */
288     public BigMatrix scalarMultiply(BigDecimal d) {
289         int rowCount = this.getRowDimension();
290         int columnCount = this.getColumnDimension();
291         BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
292         for (int row = 0; row < rowCount; row++) {
293             for (int col = 0; col < columnCount; col++) {
294                 outData[row][col] = data[row][col].multiply(d);
295             }
296         }
297         return new BigMatrixImpl(outData);
298     }
299     
300     /***
301      * Returns the result of postmultiplying this by <code>m</code>.
302      * @param m    matrix to postmultiply by
303      * @return     this*m
304      * @throws     IllegalArgumentException
305      *             if columnDimension(this) != rowDimension(m)
306      */
307     public BigMatrix multiply(BigMatrix m) throws IllegalArgumentException {
308         if (this.getColumnDimension() != m.getRowDimension()) {
309             throw new IllegalArgumentException("Matrices are not multiplication compatible.");
310         }
311         int nRows = this.getRowDimension();
312         int nCols = m.getColumnDimension();
313         int nSum = this.getColumnDimension();
314         BigDecimal[][] outData = new BigDecimal[nRows][nCols];
315         BigDecimal sum = ZERO;
316         for (int row = 0; row < nRows; row++) {
317             for (int col = 0; col < nCols; col++) {
318                 sum = ZERO;
319                 for (int i = 0; i < nSum; i++) {
320                     sum = sum.add(data[row][i].multiply(m.getEntry(i, col)));
321                 }
322                 outData[row][col] = sum;
323             }
324         }
325         return new BigMatrixImpl(outData);
326     }
327     
328     /***
329      * Returns the result premultiplying this by <code>m</code>.
330      * @param m    matrix to premultiply by
331      * @return     m * this
332      * @throws     IllegalArgumentException
333      *             if rowDimension(this) != columnDimension(m)
334      */
335     public BigMatrix preMultiply(BigMatrix m) throws IllegalArgumentException {
336         return m.multiply(this);
337     }
338     
339     /***
340      * Returns matrix entries as a two-dimensional array.
341      * <p>
342      * Makes a fresh copy of the underlying data.
343      *
344      * @return    2-dimensional array of entries
345      */
346     public BigDecimal[][] getData() {
347         return copyOut();
348     }
349     
350     /***
351      * Returns matrix entries as a two-dimensional array.
352      * <p>
353      * Makes a fresh copy of the underlying data converted to
354      * <code>double</code> values.
355      *
356      * @return    2-dimensional array of entries
357      */
358     public double[][] getDataAsDoubleArray() {
359         int nRows = getRowDimension();
360         int nCols = getColumnDimension();
361         double d[][] = new double[nRows][nCols];
362         for (int i = 0; i < nRows; i++) {
363             for (int j=0; j<nCols;j++) {
364                 d[i][j] = data[i][j].doubleValue();
365             }
366         }
367         return d;
368     }
369     
370     /***
371      * Returns a reference to the underlying data array.
372      * <p>
373      * Does not make a fresh copy of the underlying data.
374      *
375      * @return 2-dimensional array of entries
376      */
377     public BigDecimal[][] getDataRef() {
378         return data;
379     }
380     
381     /****
382      * Gets the rounding mode for division operations
383      * The default is {@link java.math.BigDecimal#ROUND_HALF_UP}
384      * @see BigDecimal
385      * @return the rounding mode.
386      */ 
387     public int getRoundingMode() {
388         return roundingMode;
389     }
390     
391     /****
392      * Sets the rounding mode for decimal divisions.
393      * @see BigDecimal
394      * @param roundingMode
395      */
396     public void setRoundingMode(int roundingMode) {
397         this.roundingMode = roundingMode;
398     }
399     
400     /****
401      * Sets the scale for division operations.
402      * The default is 64
403      * @see BigDecimal
404      * @return the scale
405      */
406     public int getScale() {
407         return scale;
408     }
409     
410     /****
411      * Sets the scale for division operations.
412      * @see BigDecimal
413      * @param scale
414      */
415     public void setScale(int scale) {
416         this.scale = scale;
417     }
418     
419     /***
420      * Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html">
421      * maximum absolute row sum norm</a> of the matrix.
422      *
423      * @return norm
424      */
425     public BigDecimal getNorm() {
426         BigDecimal maxColSum = ZERO;
427         for (int col = 0; col < this.getColumnDimension(); col++) {
428             BigDecimal sum = ZERO;
429             for (int row = 0; row < this.getRowDimension(); row++) {
430                 sum = sum.add(data[row][col].abs());
431             }
432             maxColSum = maxColSum.max(sum);
433         }
434         return maxColSum;
435     }
436     
437     /***
438      * Gets a submatrix. Rows and columns are indicated
439      * counting from 0 to n-1.
440      *
441      * @param startRow Initial row index
442      * @param endRow Final row index
443      * @param startColumn Initial column index
444      * @param endColumn Final column index
445      * @return The subMatrix containing the data of the
446      *         specified rows and columns
447      * @exception MatrixIndexException if row or column selections are not valid
448      */
449     public BigMatrix getSubMatrix(int startRow, int endRow, int startColumn,
450             int endColumn) throws MatrixIndexException {
451         if (startRow < 0 || startRow > endRow || endRow > data.length ||
452                 startColumn < 0 || startColumn > endColumn ||
453                 endColumn > data[0].length ) {
454             throw new MatrixIndexException(
455             "invalid row or column index selection");
456         }
457         BigMatrixImpl subMatrix = new BigMatrixImpl(endRow - startRow+1,
458                 endColumn - startColumn+1);
459         BigDecimal[][] subMatrixData = subMatrix.getDataRef();
460         for (int i = startRow; i <= endRow; i++) {
461             for (int j = startColumn; j <= endColumn; j++) {
462                 subMatrixData[i - startRow][j - startColumn] = data[i][j];
463             }
464         }
465         return subMatrix;
466     }
467     
468     /***
469      * Gets a submatrix. Rows and columns are indicated
470      * counting from 0 to n-1.
471      *
472      * @param selectedRows Array of row indices must be non-empty
473      * @param selectedColumns Array of column indices must be non-empty
474      * @return The subMatrix containing the data in the
475      *     specified rows and columns
476      * @exception MatrixIndexException  if supplied row or column index arrays
477      *     are not valid
478      */
479     public BigMatrix getSubMatrix(int[] selectedRows, int[] selectedColumns)
480     throws MatrixIndexException {
481         if (selectedRows.length * selectedColumns.length == 0) {
482             throw new MatrixIndexException(
483             "selected row and column index arrays must be non-empty");
484         }
485         BigMatrixImpl subMatrix = new BigMatrixImpl(selectedRows.length,
486                 selectedColumns.length);
487         BigDecimal[][] subMatrixData = subMatrix.getDataRef();
488         try  {
489             for (int i = 0; i < selectedRows.length; i++) {
490                 for (int j = 0; j < selectedColumns.length; j++) {
491                     subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
492                 }
493             }
494         }
495         catch (ArrayIndexOutOfBoundsException e) {
496             throw new MatrixIndexException("matrix dimension mismatch");
497         }
498         return subMatrix;
499     } 
500     
501     /***
502      * Returns the entries in row number <code>row</code>
503      * as a row matrix.  Row indices start at 0.
504      *
505      * @param row the row to be fetched
506      * @return row matrix
507      * @throws MatrixIndexException if the specified row index is invalid
508      */
509     public BigMatrix getRowMatrix(int row) throws MatrixIndexException {
510         if ( !isValidCoordinate( row, 0)) {
511             throw new MatrixIndexException("illegal row argument");
512         }
513         int ncols = this.getColumnDimension();
514         BigDecimal[][] out = new BigDecimal[1][ncols]; 
515         System.arraycopy(data[row], 0, out[0], 0, ncols);
516         return new BigMatrixImpl(out);
517     } 
518     
519     /***
520      * Returns the entries in column number <code>column</code>
521      * as a column matrix.  Column indices start at 0.
522      *
523      * @param column the column to be fetched
524      * @return column matrix
525      * @throws MatrixIndexException if the specified column index is invalid
526      */
527     public BigMatrix getColumnMatrix(int column) throws MatrixIndexException {
528         if ( !isValidCoordinate( 0, column)) {
529             throw new MatrixIndexException("illegal column argument");
530         }
531         int nRows = this.getRowDimension();
532         BigDecimal[][] out = new BigDecimal[nRows][1]; 
533         for (int row = 0; row < nRows; row++) {
534             out[row][0] = data[row][column];
535         }
536         return new BigMatrixImpl(out);
537     }
538     
539     /***
540      * Returns the entries in row number <code>row</code> as an array.
541      * <p>
542      * Row indices start at 0.  A <code>MatrixIndexException</code> is thrown
543      * unless <code>0 <= row < rowDimension.</code>
544      *
545      * @param row the row to be fetched
546      * @return array of entries in the row
547      * @throws MatrixIndexException if the specified row index is not valid
548      */
549     public BigDecimal[] getRow(int row) throws MatrixIndexException {
550         if ( !isValidCoordinate( row, 0 ) ) {
551             throw new MatrixIndexException("illegal row argument");
552         }
553         int ncols = this.getColumnDimension();
554         BigDecimal[] out = new BigDecimal[ncols];
555         System.arraycopy(data[row], 0, out, 0, ncols);
556         return out;
557     }
558     
559      /***
560      * Returns the entries in row number <code>row</code> as an array
561      * of double values.
562      * <p>
563      * Row indices start at 0.  A <code>MatrixIndexException</code> is thrown
564      * unless <code>0 <= row < rowDimension.</code>
565      *
566      * @param row the row to be fetched
567      * @return array of entries in the row
568      * @throws MatrixIndexException if the specified row index is not valid
569      */
570     public double[] getRowAsDoubleArray(int row) throws MatrixIndexException {
571         if ( !isValidCoordinate( row, 0 ) ) {
572             throw new MatrixIndexException("illegal row argument");
573         }
574         int ncols = this.getColumnDimension();
575         double[] out = new double[ncols];
576         for (int i=0;i<ncols;i++) {
577             out[i] = data[row][i].doubleValue();
578         }
579         return out;
580     }
581     
582      /***
583      * Returns the entries in column number <code>col</code> as an array.
584      * <p>
585      * Column indices start at 0.  A <code>MatrixIndexException</code> is thrown
586      * unless <code>0 <= column < columnDimension.</code>
587      *
588      * @param col the column to be fetched
589      * @return array of entries in the column
590      * @throws MatrixIndexException if the specified column index is not valid
591      */
592     public BigDecimal[] getColumn(int col) throws MatrixIndexException {
593         if ( !isValidCoordinate(0, col) ) {
594             throw new MatrixIndexException("illegal column argument");
595         }
596         int nRows = this.getRowDimension();
597         BigDecimal[] out = new BigDecimal[nRows];
598         for (int i = 0; i < nRows; i++) {
599             out[i] = data[i][col];
600         }
601         return out;
602     }
603     
604     /***
605      * Returns the entries in column number <code>col</code> as an array
606      * of double values.
607      * <p>
608      * Column indices start at 0.  A <code>MatrixIndexException</code> is thrown
609      * unless <code>0 <= column < columnDimension.</code>
610      *
611      * @param col the column to be fetched
612      * @return array of entries in the column
613      * @throws MatrixIndexException if the specified column index is not valid
614      */
615     public double[] getColumnAsDoubleArray(int col) throws MatrixIndexException {
616         if ( !isValidCoordinate( 0, col ) ) {
617             throw new MatrixIndexException("illegal column argument");
618         }
619         int nrows = this.getRowDimension();
620         double[] out = new double[nrows];
621         for (int i=0;i<nrows;i++) {
622             out[i] = data[i][col].doubleValue();
623         }
624         return out;
625     }
626     
627      /***
628      * Returns the entry in the specified row and column.
629      * <p>
630      * Row and column indices start at 0 and must satisfy 
631      * <ul>
632      * <li><code>0 <= row < rowDimension</code></li>
633      * <li><code> 0 <= column < columnDimension</code></li>
634      * </ul>
635      * otherwise a <code>MatrixIndexException</code> is thrown.
636      *
637      * @param row  row location of entry to be fetched  
638      * @param column  column location of entry to be fetched
639      * @return matrix entry in row,column
640      * @throws MatrixIndexException if the row or column index is not valid
641      */
642     public BigDecimal getEntry(int row, int column)
643     throws MatrixIndexException {
644         if (!isValidCoordinate(row,column)) {
645             throw new MatrixIndexException("matrix entry does not exist");
646         }
647         return data[row][column];
648     }
649     
650     /***
651      * Returns the entry in the specified row and column as a double.
652      * <p>
653      * Row and column indices start at 0 and must satisfy 
654      * <ul>
655      * <li><code>0 <= row < rowDimension</code></li>
656      * <li><code> 0 <= column < columnDimension</code></li>
657      * </ul>
658      * otherwise a <code>MatrixIndexException</code> is thrown.
659      *
660      * @param row  row location of entry to be fetched
661      * @param column  column location of entry to be fetched
662      * @return matrix entry in row,column
663      * @throws MatrixIndexException if the row
664      * or column index is not valid
665      */
666     public double getEntryAsDouble(int row, int column) throws MatrixIndexException {
667         return getEntry(row,column).doubleValue();
668     }
669     
670     /***
671      * Returns the transpose matrix.
672      *
673      * @return transpose matrix
674      */
675     public BigMatrix transpose() {
676         int nRows = this.getRowDimension();
677         int nCols = this.getColumnDimension();
678         BigMatrixImpl out = new BigMatrixImpl(nCols, nRows);
679         BigDecimal[][] outData = out.getDataRef();
680         for (int row = 0; row < nRows; row++) {
681             for (int col = 0; col < nCols; col++) {
682                 outData[col][row] = data[row][col];
683             }
684         }
685         return out;
686     }
687     
688     /***
689      * Returns the inverse matrix if this matrix is invertible.
690      * 
691      * @return inverse matrix
692      * @throws InvalidMatrixException if this is not invertible
693      */
694     public BigMatrix inverse() throws InvalidMatrixException {
695         return solve(getIdentity(this.getRowDimension()));
696     }
697     
698     /***
699      * Returns the determinant of this matrix.
700      *
701      * @return determinant
702      * @throws InvalidMatrixException if matrix is not square
703      */
704     public BigDecimal getDeterminant() throws InvalidMatrixException {
705         if (!isSquare()) {
706             throw new InvalidMatrixException("matrix is not square");
707         }
708         if (isSingular()) {   // note: this has side effect of attempting LU decomp if lu == null
709             return ZERO;
710         } else {
711             BigDecimal det = (parity == 1) ? ONE : ONE.negate();
712             for (int i = 0; i < this.getRowDimension(); i++) {
713                 det = det.multiply(lu[i][i]);
714             }
715             return det;
716         }
717     }
718     
719      /***
720      * Is this a square matrix?
721      * @return true if the matrix is square (rowDimension = columnDimension)
722      */
723     public boolean isSquare() {
724         return (this.getColumnDimension() == this.getRowDimension());
725     }
726     
727     /***
728      * Is this a singular matrix?
729      * @return true if the matrix is singular
730      */
731     public boolean isSingular() {
732         if (lu == null) {
733             try {
734                 luDecompose();
735                 return false;
736             } catch (InvalidMatrixException ex) {
737                 return true;
738             }
739         } else { // LU decomp must have been successfully performed
740             return false; // so the matrix is not singular
741         }
742     }
743     
744     /***
745      * Returns the number of rows in the matrix.
746      *
747      * @return rowDimension
748      */
749     public int getRowDimension() {
750         return data.length;
751     }
752     
753     /***
754      * Returns the number of columns in the matrix.
755      *
756      * @return columnDimension
757      */
758     public int getColumnDimension() {
759         return data[0].length;
760     }
761     
762      /***
763      * Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
764      * trace</a> of the matrix (the sum of the elements on the main diagonal).
765      *
766      * @return trace
767      * 
768      * @throws IllegalArgumentException if this matrix is not square.
769      */
770     public BigDecimal getTrace() throws IllegalArgumentException {
771         if (!isSquare()) {
772             throw new IllegalArgumentException("matrix is not square");
773         }
774         BigDecimal trace = data[0][0];
775         for (int i = 1; i < this.getRowDimension(); i++) {
776             trace = trace.add(data[i][i]);
777         }
778         return trace;
779     }
780     
781     /***
782      * Returns the result of multiplying this by the vector <code>v</code>.
783      *
784      * @param v the vector to operate on
785      * @return this*v
786      * @throws IllegalArgumentException if columnDimension != v.size()
787      */
788     public BigDecimal[] operate(BigDecimal[] v) throws IllegalArgumentException {
789         if (v.length != this.getColumnDimension()) {
790             throw new IllegalArgumentException("vector has wrong length");
791         }
792         int nRows = this.getRowDimension();
793         int nCols = this.getColumnDimension();
794         BigDecimal[] out = new BigDecimal[v.length];
795         for (int row = 0; row < nRows; row++) {
796             BigDecimal sum = ZERO;
797             for (int i = 0; i < nCols; i++) {
798                 sum = sum.add(data[row][i].multiply(v[i]));
799             }
800             out[row] = sum;
801         }
802         return out;
803     }
804     
805     /***
806      * Returns the result of multiplying this by the vector <code>v</code>.
807      *
808      * @param v the vector to operate on
809      * @return this*v
810      * @throws IllegalArgumentException if columnDimension != v.size()
811      */
812     public BigDecimal[] operate(double[] v) throws IllegalArgumentException {
813         BigDecimal bd[] = new BigDecimal[v.length];
814         for (int i=0;i<bd.length;i++) {
815             bd[i] = new BigDecimal(v[i]);
816         }
817         return operate(bd);
818     }
819     
820     /***
821      * Returns the (row) vector result of premultiplying this by the vector <code>v</code>.
822      *
823      * @param v the row vector to premultiply by
824      * @return v*this
825      * @throws IllegalArgumentException if rowDimension != v.size()
826      */
827     public BigDecimal[] preMultiply(BigDecimal[] v) throws IllegalArgumentException {
828         int nRows = this.getRowDimension();
829         if (v.length != nRows) {
830             throw new IllegalArgumentException("vector has wrong length");
831         }
832         int nCols = this.getColumnDimension();
833         BigDecimal[] out = new BigDecimal[nCols];
834         for (int col = 0; col < nCols; col++) {
835             BigDecimal sum = ZERO;
836             for (int i = 0; i < nRows; i++) {
837                 sum = sum.add(data[i][col].multiply(v[i]));
838             }
839             out[col] = sum;
840         }
841         return out;
842     }
843     
844     /***
845      * Returns a matrix of (column) solution vectors for linear systems with
846      * coefficient matrix = this and constant vectors = columns of
847      * <code>b</code>. 
848      *
849      * @param b  array of constants forming RHS of linear systems to
850      * to solve
851      * @return solution array
852      * @throws IllegalArgumentException if this.rowDimension != row dimension
853      * @throws InvalidMatrixException if this matrix is not square or is singular
854      */
855     public BigDecimal[] solve(BigDecimal[] b) throws IllegalArgumentException, InvalidMatrixException {
856         int nRows = this.getRowDimension();
857         if (b.length != nRows) {
858             throw new IllegalArgumentException("constant vector has wrong length");
859         }
860         BigMatrix bMatrix = new BigMatrixImpl(b);
861         BigDecimal[][] solution = ((BigMatrixImpl) (solve(bMatrix))).getDataRef();
862         BigDecimal[] out = new BigDecimal[nRows];
863         for (int row = 0; row < nRows; row++) {
864             out[row] = solution[row][0];
865         }
866         return out;
867     }
868     
869     /***
870      * Returns a matrix of (column) solution vectors for linear systems with
871      * coefficient matrix = this and constant vectors = columns of
872      * <code>b</code>. 
873      *
874      * @param b  array of constants forming RHS of linear systems to
875      * to solve
876      * @return solution array
877      * @throws IllegalArgumentException if this.rowDimension != row dimension
878      * @throws InvalidMatrixException if this matrix is not square or is singular
879      */
880     public BigDecimal[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
881         BigDecimal bd[] = new BigDecimal[b.length];
882         for (int i=0;i<bd.length;i++) {
883             bd[i] = new BigDecimal(b[i]);
884         }
885         return solve(bd);
886     }
887     
888     /***
889      * Returns a matrix of (column) solution vectors for linear systems with
890      * coefficient matrix = this and constant vectors = columns of
891      * <code>b</code>. 
892      *
893      * @param b  matrix of constant vectors forming RHS of linear systems to
894      * to solve
895      * @return matrix of solution vectors
896      * @throws IllegalArgumentException if this.rowDimension != row dimension
897      * @throws InvalidMatrixException if this matrix is not square or is singular
898      */
899     public BigMatrix solve(BigMatrix b) throws IllegalArgumentException, InvalidMatrixException  {
900         if (b.getRowDimension() != this.getRowDimension()) {
901             throw new IllegalArgumentException("Incorrect row dimension");
902         }
903         if (!this.isSquare()) {
904             throw new InvalidMatrixException("coefficient matrix is not square");
905         }
906         if (this.isSingular()) { // side effect: compute LU decomp
907             throw new InvalidMatrixException("Matrix is singular.");
908         }
909         
910         int nCol = this.getColumnDimension();
911         int nColB = b.getColumnDimension();
912         int nRowB = b.getRowDimension();
913         
914         // Apply permutations to b
915         BigDecimal[][] bp = new BigDecimal[nRowB][nColB];
916         for (int row = 0; row < nRowB; row++) {
917             for (int col = 0; col < nColB; col++) {
918                 bp[row][col] = b.getEntry(permutation[row], col);
919             }
920         }
921         
922         // Solve LY = b
923         for (int col = 0; col < nCol; col++) {
924             for (int i = col + 1; i < nCol; i++) {
925                 for (int j = 0; j < nColB; j++) {
926                     bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
927                 }
928             }
929         }
930         
931         // Solve UX = Y
932         for (int col = nCol - 1; col >= 0; col--) {
933             for (int j = 0; j < nColB; j++) {
934                 bp[col][j] = bp[col][j].divide(lu[col][col], scale, roundingMode);
935             }
936             for (int i = 0; i < col; i++) {
937                 for (int j = 0; j < nColB; j++) {
938                     bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
939                 }
940             }
941         }
942         
943         BigMatrixImpl outMat = new BigMatrixImpl(bp);
944         return outMat;
945     }
946     
947     /***
948      * Computes a new 
949      * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
950      * LU decompostion</a> for this matrix, storing the result for use by other methods. 
951      * <p>
952      * <strong>Implementation Note</strong>:<br>
953      * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
954      * Crout's algortithm</a>, with partial pivoting.
955      * <p>
956      * <strong>Usage Note</strong>:<br>
957      * This method should rarely be invoked directly. Its only use is
958      * to force recomputation of the LU decomposition when changes have been
959      * made to the underlying data using direct array references. Changes
960      * made using setXxx methods will trigger recomputation when needed
961      * automatically.
962      *
963      * @throws InvalidMatrixException if the matrix is non-square or singular.
964      */
965     public void luDecompose() throws InvalidMatrixException {
966         
967         int nRows = this.getRowDimension();
968         int nCols = this.getColumnDimension();
969         if (nRows != nCols) {
970             throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
971         }
972         lu = this.getData();
973         
974         // Initialize permutation array and parity
975         permutation = new int[nRows];
976         for (int row = 0; row < nRows; row++) {
977             permutation[row] = row;
978         }
979         parity = 1;
980         
981         // Loop over columns
982         for (int col = 0; col < nCols; col++) {
983             
984             BigDecimal sum = ZERO;
985             
986             // upper
987             for (int row = 0; row < col; row++) {
988                 sum = lu[row][col];
989                 for (int i = 0; i < row; i++) {
990                     sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
991                 }
992                 lu[row][col] = sum;
993             }
994             
995             // lower
996             int max = col; // permutation row
997             BigDecimal largest = ZERO;
998             for (int row = col; row < nRows; row++) {
999                 sum = lu[row][col];
1000                 for (int i = 0; i < col; i++) {
1001                     sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
1002                 }
1003                 lu[row][col] = sum;
1004                 
1005                 // maintain best permutation choice
1006                 if (sum.abs().compareTo(largest) == 1) {
1007                     largest = sum.abs();
1008                     max = row;
1009                 }
1010             }
1011             
1012             // Singularity check
1013             if (lu[max][col].abs().compareTo(TOO_SMALL) <= 0) {
1014                 lu = null;
1015                 throw new InvalidMatrixException("matrix is singular");
1016             }
1017             
1018             // Pivot if necessary
1019             if (max != col) {
1020                 BigDecimal tmp = ZERO;
1021                 for (int i = 0; i < nCols; i++) {
1022                     tmp = lu[max][i];
1023                     lu[max][i] = lu[col][i];
1024                     lu[col][i] = tmp;
1025                 }
1026                 int temp = permutation[max];
1027                 permutation[max] = permutation[col];
1028                 permutation[col] = temp;
1029                 parity = -parity;
1030             }
1031             
1032             //Divide the lower elements by the "winning" diagonal elt.
1033             for (int row = col + 1; row < nRows; row++) {
1034                 lu[row][col] = lu[row][col].divide(lu[col][col], scale, roundingMode);
1035             }
1036             
1037         }
1038         
1039     }
1040     
1041     /***
1042      * 
1043      * @see Object#toString()
1044      */
1045     public String toString() {
1046         StringBuffer res = new StringBuffer();
1047         res.append("BigMatrixImpl{");
1048         if (data != null) {
1049             for (int i = 0; i < data.length; i++) {
1050                 if (i > 0)
1051                     res.append(",");
1052                 res.append("{");
1053                 for (int j = 0; j < data[0].length; j++) {
1054                     if (j > 0)
1055                         res.append(",");
1056                     res.append(data[i][j]);
1057                 } 
1058                 res.append("}");
1059             } 
1060         }
1061         res.append("}");
1062         return res.toString();
1063     } 
1064     
1065     /***
1066      * Returns true iff <code>object</code> is a 
1067      * <code>BigMatrixImpl</code> instance with the same dimensions as this
1068      * and all corresponding matrix entries are equal.  BigDecimal.equals
1069      * is used to compare corresponding entries.
1070      * 
1071      * @param object the object to test equality against.
1072      * @return true if object equals this
1073      */
1074     public boolean equals(Object object) {
1075         if (object == this ) {
1076             return true;
1077         }
1078         if (object instanceof BigMatrixImpl == false) {
1079             return false;
1080         }
1081         BigMatrix m = (BigMatrix) object;
1082         int nRows = getRowDimension();
1083         int nCols = getColumnDimension();
1084         if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
1085             return false;
1086         }
1087         for (int row = 0; row < nRows; row++) {
1088             for (int col = 0; col < nCols; col++) {
1089                 if (!data[row][col].equals(m.getEntry(row, col))) {
1090                     return false;
1091                 }
1092             }
1093         }
1094         return true;
1095     }
1096     
1097     /***
1098      * Computes a hashcode for the matrix.
1099      * 
1100      * @return hashcode for matrix
1101      */
1102     public int hashCode() {
1103         int ret = 7;
1104         int nRows = getRowDimension();
1105         int nCols = getColumnDimension();
1106         ret = ret * 31 + nRows;
1107         ret = ret * 31 + nCols;
1108         for (int row = 0; row < nRows; row++) {
1109             for (int col = 0; col < nCols; col++) {
1110                 ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) * 
1111                 data[row][col].hashCode();
1112             }
1113         }   
1114         return ret;
1115     }
1116     
1117     //------------------------ Protected methods
1118     
1119     /***
1120      * Returns <code>dimension x dimension</code> identity matrix.
1121      *
1122      * @param dimension dimension of identity matrix to generate
1123      * @return identity matrix
1124      */
1125     protected BigMatrix getIdentity(int dimension) {
1126         BigMatrixImpl out = new BigMatrixImpl(dimension, dimension);
1127         BigDecimal[][] d = out.getDataRef();
1128         for (int row = 0; row < dimension; row++) {
1129             for (int col = 0; col < dimension; col++) {
1130                 d[row][col] = row == col ? ONE : ZERO;
1131             }
1132         }
1133         return out;
1134     }
1135     
1136     /***
1137      *  Returns the LU decomposition as a BigMatrix.
1138      *  Returns a fresh copy of the cached LU matrix if this has been computed; 
1139      *  otherwise the composition is computed and cached for use by other methods.   
1140      *  Since a copy is returned in either case, changes to the returned matrix do not 
1141      *  affect the LU decomposition property. 
1142      * <p>
1143      * The matrix returned is a compact representation of the LU decomposition. 
1144      * Elements below the main diagonal correspond to entries of the "L" matrix;   
1145      * elements on and above the main diagonal correspond to entries of the "U"
1146      * matrix.
1147      * <p>
1148      * Example: <pre>
1149      * 
1150      *     Returned matrix                L                  U
1151      *         2  3  1                   1  0  0            2  3  1          
1152      *         5  4  6                   5  1  0            0  4  6
1153      *         1  7  8                   1  7  1            0  0  8          
1154      * </pre>
1155      * 
1156      * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
1157      *  where permuteRows reorders the rows of the matrix to follow the order determined
1158      *  by the <a href=#getPermutation()>permutation</a> property.
1159      * 
1160      * @return LU decomposition matrix
1161      * @throws InvalidMatrixException if the matrix is non-square or singular.
1162      */
1163     protected BigMatrix getLUMatrix() throws InvalidMatrixException {
1164         if (lu == null) {
1165             luDecompose();
1166         }
1167         return new BigMatrixImpl(lu);
1168     }
1169     
1170     /***
1171      * Returns the permutation associated with the lu decomposition.
1172      * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
1173      * <p>
1174      * Example:
1175      * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
1176      * and current first row is last.
1177      * <p>
1178      * Returns a fresh copy of the array.
1179      * 
1180      * @return the permutation
1181      */
1182     protected int[] getPermutation() {
1183         int[] out = new int[permutation.length];
1184         System.arraycopy(permutation, 0, out, 0, permutation.length);
1185         return out;
1186     }
1187     
1188     //------------------------ Private methods
1189     
1190     /***
1191      * Returns a fresh copy of the underlying data array.
1192      *
1193      * @return a copy of the underlying data array.
1194      */
1195     private BigDecimal[][] copyOut() {
1196         int nRows = this.getRowDimension();
1197         BigDecimal[][] out = new BigDecimal[nRows][this.getColumnDimension()];
1198         // can't copy 2-d array in one shot, otherwise get row references
1199         for (int i = 0; i < nRows; i++) {
1200             System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1201         }
1202         return out;
1203     }
1204     
1205     /***
1206      * Replaces data with a fresh copy of the input array.
1207      *
1208      * @param in data to copy in
1209      */
1210     private void copyIn(BigDecimal[][] in) {
1211         int nRows = in.length;
1212         int nCols = in[0].length;
1213         data = new BigDecimal[nRows][nCols];
1214         System.arraycopy(in, 0, data, 0, in.length);
1215         for (int i = 0; i < nRows; i++) {
1216             System.arraycopy(in[i], 0, data[i], 0, nCols);
1217         }
1218         lu = null;
1219     }
1220     
1221     /***
1222      * Replaces data with a fresh copy of the input array.
1223      *
1224      * @param in data to copy in
1225      */
1226     private void copyIn(double[][] in) {
1227         int nRows = in.length;
1228         int nCols = in[0].length;
1229         data = new BigDecimal[nRows][nCols];
1230         for (int i = 0; i < nRows; i++) {
1231             for (int j=0; j < nCols; j++) {
1232                 data[i][j] = new BigDecimal(in[i][j]);
1233             }
1234         }
1235         lu = null;
1236     }
1237     
1238     /***
1239      * Replaces data with BigDecimals represented by the strings in the input
1240      * array.
1241      *
1242      * @param in data to copy in
1243      */
1244     private void copyIn(String[][] in) {
1245         int nRows = in.length;
1246         int nCols = in[0].length;
1247         data = new BigDecimal[nRows][nCols];
1248         for (int i = 0; i < nRows; i++) {
1249             for (int j=0; j < nCols; j++) {
1250                 data[i][j] = new BigDecimal(in[i][j]);
1251             }
1252         }
1253         lu = null;
1254     }
1255     
1256     /***
1257      * Tests a given coordinate as being valid or invalid
1258      *
1259      * @param row the row index.
1260      * @param col the column index.
1261      * @return true if the coordinate is with the current dimensions
1262      */
1263     private boolean isValidCoordinate(int row, int col) {
1264         int nRows = this.getRowDimension();
1265         int nCols = this.getColumnDimension();
1266         
1267         return !(row < 0 || row >= nRows || col < 0 || col >= nCols);
1268     }
1269     
1270 }