    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    import org.apache.commons.math.exception.util.LocalizedFormats;
    024
    025    /**
    026     * Implementation of {@link BigMatrix} using a BigDecimal[][] array to store entries
    027     * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
    028     * LU decompostion</a> to support linear system
    029     * solution and inverse.
    030     * <p>
    031     * The LU decompostion is performed as needed, to support the following operations: <ul>
    032     * <li>solve</li>
    033     * <li>isSingular</li>
    034     * <li>getDeterminant</li>
    035     * <li>inverse</li> </ul></p>
    036     * <p>
    037    * <strong>Usage notes</strong>:<br>
    038     * <ul><li>
    039     * The LU decomposition is stored and reused on subsequent calls.  If matrix
    040     * data are modified using any of the public setXxx methods, the saved
    041     * decomposition is discarded.  If data are modified via references to the
    042     * underlying array obtained using <code>getDataRef()</code>, then the stored
    043     * LU decomposition will not be discarded.  In this case, you need to
    044     * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
    045     * before using any of the methods above.</li>
    046     * <li>
    047     * As specified in the {@link BigMatrix} interface, matrix element indexing
    048     * is 0-based -- e.g., <code>getEntry(0, 0)</code>
    049     * returns the element in the first row, first column of the matrix.</li></ul></p>
    050     *
    051     * @deprecated as of 2.0, replaced by {@link Array2DRowFieldMatrix} with a {@link
    052     * org.apache.commons.math.util.BigReal} parameter
    053     * @version $Revision: 1042376 $ $Date: 2010-12-05 16:54:55 +0100 (dim. 05 d??c. 2010) $
    054     */
    055    @Deprecated
    056    public class BigMatrixImpl implements BigMatrix, Serializable {
    057
    058        /** BigDecimal 0 */
    059        static final BigDecimal ZERO = new BigDecimal(0);
    060
    061        /** BigDecimal 1 */
    062        static final BigDecimal ONE = new BigDecimal(1);
    063
    064        /** Bound to determine effective singularity in LU decomposition */
    065        private static final BigDecimal TOO_SMALL = new BigDecimal(10E-12);
    066
    067        /** Serialization id */
    068        private static final long serialVersionUID = -1011428905656140431L;
    069
    070        /** Entries of the matrix */
    071        protected BigDecimal data[][] = null;
    072
    073        /** Entries of cached LU decomposition.
    074         *  All updates to data (other than luDecompose()) *must* set this to null
    075         */
    076        protected BigDecimal lu[][] = null;
    077
    078        /** Permutation associated with LU decomposition */
    079        protected int[] permutation = null;
    080
    081        /** Parity of the permutation associated with the LU decomposition */
    082        protected int parity = 1;
    083
    084        /** Rounding mode for divisions **/
    085        private int roundingMode = BigDecimal.ROUND_HALF_UP;
    086
    087        /*** BigDecimal scale ***/
    088        private int scale = 64;
    089
    090        /**
    091         * Creates a matrix with no data
    092         */
    093        public BigMatrixImpl() {
    094        }
    095
    096        /**
    097         * Create a new BigMatrix with the supplied row and column dimensions.
    098         *
    099         * @param rowDimension      the number of rows in the new matrix
    100         * @param columnDimension   the number of columns in the new matrix
    101         * @throws IllegalArgumentException if row or column dimension is not
    102         *  positive
    103         */
    104        public BigMatrixImpl(int rowDimension, int columnDimension) {
    105            if (rowDimension < 1 ) {
    106                throw MathRuntimeException.createIllegalArgumentException(
    107                        LocalizedFormats.INSUFFICIENT_DIMENSION, rowDimension, 1);
    108            }
    109            if (columnDimension < 1) {
    110                throw MathRuntimeException.createIllegalArgumentException(
    111                        LocalizedFormats.INSUFFICIENT_DIMENSION, columnDimension, 1);
    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(LocalizedFormats.AT_LEAST_ONE_ROW);
    159                }
    160
    161                final int nCols = d[0].length;
    162                if (nCols == 0) {
    163                    throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
    164                }
    165                for (int r = 1; r < nRows; r++) {
    166                    if (d[r].length != nCols) {
    167                        throw MathRuntimeException.createIllegalArgumentException(
    168                              LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    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(LocalizedFormats.AT_LEAST_ONE_ROW);
    192            }
    193
    194            final int nCols = d[0].length;
    195            if (nCols == 0) {
    196                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
    197            }
    198            for (int row = 1; row < nRows; row++) {
    199                if (d[row].length != nCols) {
    200                    throw MathRuntimeException.createIllegalArgumentException(
    201                          LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    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(LocalizedFormats.AT_LEAST_ONE_ROW);
    222            }
    223
    224            final int nCols = d[0].length;
    225            if (nCols == 0) {
    226                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
    227            }
    228            for (int row = 1; row < nRows; row++) {
    229                if (d[row].length != nCols) {
    230                    throw MathRuntimeException.createIllegalArgumentException(
    231                          LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    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(LocalizedFormats.INITIAL_ROW_AFTER_FINAL_ROW,
    607                                               startRow, endRow);
    608            }
    609
    610            MatrixUtils.checkColumnIndex(this, startColumn);
    611            MatrixUtils.checkColumnIndex(this, endColumn);
    612            if (startColumn > endColumn) {
    613                throw new MatrixIndexException(LocalizedFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
    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(LocalizedFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
    646                }
    647                throw new MatrixIndexException(LocalizedFormats.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(LocalizedFormats.AT_LEAST_ONE_ROW);
    706            }
    707
    708            final int nCols = subMatrix[0].length;
    709            if (nCols == 0) {
    710                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.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                          LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    717                          nCols, subMatrix[r].length);
    718                }
    719            }
    720
    721            if (data == null) {
    722                if (row > 0) {
    723                    throw MathRuntimeException.createIllegalStateException(
    724                            LocalizedFormats.FIRST_ROWS_NOT_INITIALIZED_YET,
    725                            row);
    726                }
    727                if (column > 0) {
    728                    throw MathRuntimeException.createIllegalStateException(
    729                            LocalizedFormats.FIRST_COLUMNS_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                        LocalizedFormats.NO_SUCH_MATRIX_ENTRY,
    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 < 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 getColumnDimension() == 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                        LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    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                        LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    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                        LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    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                        LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    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            StringBuilder res = new StringBuilder();
    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    }