View Javadoc

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