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