001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    import java.io.Serializable;
020    import java.math.BigDecimal;
021    
022    import org.apache.commons.math.MathRuntimeException;
023    
024    /**
025     * Implementation of {@link BigMatrix} using a BigDecimal[][] array to store entries
026     * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
027     * LU decompostion</a> to support linear system 
028     * solution and inverse.
029     * <p>
030     * The LU decompostion is performed as needed, to support the following operations: <ul>
031     * <li>solve</li>
032     * <li>isSingular</li>
033     * <li>getDeterminant</li>
034     * <li>inverse</li> </ul></p>
035     * <p>
036    * <strong>Usage notes</strong>:<br>
037     * <ul><li>
038     * The LU decomposition is stored and reused on subsequent calls.  If matrix
039     * data are modified using any of the public setXxx methods, the saved
040     * decomposition is discarded.  If data are modified via references to the
041     * underlying array obtained using <code>getDataRef()</code>, then the stored
042     * LU decomposition will not be discarded.  In this case, you need to
043     * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
044     * before using any of the methods above.</li>
045     * <li>
046     * As specified in the {@link BigMatrix} interface, matrix element indexing
047     * is 0-based -- e.g., <code>getEntry(0, 0)</code>
048     * returns the element in the first row, first column of the matrix.</li></ul></p>
049     * 
050     * @deprecated as of 2.0, replaced by {@link Array2DRowFieldMatrix} with a {@link
051     * org.apache.commons.math.util.BigReal} parameter
052     * @version $Revision: 783702 $ $Date: 2009-06-11 04:54:02 -0400 (Thu, 11 Jun 2009) $
053     */
054    @Deprecated
055    public class BigMatrixImpl implements BigMatrix, Serializable {
056        
057        /** Serialization id */
058        private static final long serialVersionUID = -1011428905656140431L;
059        
060        /** Entries of the matrix */
061        protected BigDecimal data[][] = null;
062        
063        /** Entries of cached LU decomposition.
064         *  All updates to data (other than luDecompose()) *must* set this to null
065         */
066        protected BigDecimal lu[][] = null;
067        
068        /** Permutation associated with LU decomposition */
069        protected int[] permutation = null;
070        
071        /** Parity of the permutation associated with the LU decomposition */
072        protected int parity = 1;
073        
074        /** Rounding mode for divisions **/
075        private int roundingMode = BigDecimal.ROUND_HALF_UP;
076        
077        /*** BigDecimal scale ***/
078        private int scale = 64;
079        
080        /** Bound to determine effective singularity in LU decomposition */
081        private static final BigDecimal TOO_SMALL = new BigDecimal(10E-12);
082        
083        /** BigDecimal 0 */
084        static final BigDecimal ZERO = new BigDecimal(0);
085        /** BigDecimal 1 */
086        static final BigDecimal ONE = new BigDecimal(1);
087        
088        /** 
089         * Creates a matrix with no data
090         */
091        public BigMatrixImpl() {
092        }
093        
094        /**
095         * Create a new BigMatrix with the supplied row and column dimensions.
096         *
097         * @param rowDimension      the number of rows in the new matrix
098         * @param columnDimension   the number of columns in the new matrix
099         * @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    }