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