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