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