View Javadoc

    1   /*
    2    * Copyright 2003-2005 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
    19  import java.io.Serializable;
    20  import org.apache.commons.math.util.MathUtils;
    21
    22
    23  /***
    24   * Implementation of RealMatrix using a double[][] array to store entries and
    25   * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
    26   * LU decompostion</a> to support linear system
    27   * solution and inverse.
    28   * <p>
    29   * The LU decompostion is performed as needed, to support the following operations: <ul>
    30   * <li>solve</li>
    31   * <li>isSingular</li>
    32   * <li>getDeterminant</li>
    33   * <li>inverse</li> </ul>
    34   * <p>
    35   * <strong>Usage notes</strong>:<br>
    36   * <ul><li>
    37   * The LU decomposition is cached and reused on subsequent calls.
    38   * If data are modified via references to the underlying array obtained using
    39   * <code>getDataRef()</code>, then the stored LU decomposition will not be
    40   * discarded.  In this case, you need to explicitly invoke
    41   * <code>LUDecompose()</code> to recompute the decomposition
    42   * before using any of the methods above.</li>
    43   * <li>
    44   * As specified in the {@link RealMatrix} interface, matrix element indexing
    45   * is 0-based -- e.g., <code>getEntry(0, 0)</code>
    46   * returns the element in the first row, first column of the matrix.</li></ul>
    47   *
    48   * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
    49   */
    50  public class RealMatrixImpl implements RealMatrix, Serializable {
    51
    52      /*** Serializable version identifier */
    53      private static final long serialVersionUID = 4237564493130426188L;
    54
    55      /*** Entries of the matrix */
    56      private double data[][] = null;
    57
    58      /*** Entries of cached LU decomposition.
    59       *  All updates to data (other than luDecompose()) *must* set this to null
    60       */
    61      private double lu[][] = null;
    62
    63      /*** Permutation associated with LU decompostion */
    64      private int[] permutation = null;
    65
    66      /*** Parity of the permutation associated with the LU decomposition */
    67      private int parity = 1;
    68
    69      /*** Bound to determine effective singularity in LU decomposition */
    70      protected static double TOO_SMALL = 10E-12;
    71
    72      /***
    73       * Creates a matrix with no data
    74       */
    75      public RealMatrixImpl() {
    76      }
    77
    78      /***
    79       * Create a new RealMatrix with the supplied row and column dimensions.
    80       *
    81       * @param rowDimension  the number of rows in the new matrix
    82       * @param columnDimension  the number of columns in the new matrix
    83       * @throws IllegalArgumentException if row or column dimension is not
    84       *  positive
    85       */
    86      public RealMatrixImpl(int rowDimension, int columnDimension) {
    87          if (rowDimension <= 0 || columnDimension <= 0) {
    88              throw new IllegalArgumentException(
    89                      "row and column dimensions must be postive");
    90          }
    91          data = new double[rowDimension][columnDimension];
    92          lu = null;
    93      }
    94
    95      /***
    96       * Create a new RealMatrix using the input array as the underlying
    97       * data array.
    98       * <p>
    99       * The input array is copied, not referenced.
    100      *
    101      * @param d data for new matrix
    102      * @throws IllegalArgumentException if <code>data</code> is not rectangular
    103      *  (not all rows have the same length) or empty
    104      * @throws NullPointerException if <code>data</code> is null
    105      */
    106     public RealMatrixImpl(double[][] d) {
    107         this.copyIn(d);
    108         lu = null;
    109     }
    110
    111     /***
    112      * Create a new (column) RealMatrix using <code>v</code> as the
    113      * data for the unique column of the <code>v.length x 1</code> matrix
    114      * created.
    115      * <p>
    116      * The input array is copied, not referenced.
    117      *
    118      * @param v column vector holding data for new matrix
    119      */
    120     public RealMatrixImpl(double[] v) {
    121         int nRows = v.length;
    122         data = new double[nRows][1];
    123         for (int row = 0; row < nRows; row++) {
    124             data[row][0] = v[row];
    125         }
    126     }
    127
    128     /***
    129      * Create a new RealMatrix which is a copy of this.
    130      *
    131      * @return  the cloned matrix
    132      */
    133     public RealMatrix copy() {
    134         return new RealMatrixImpl(this.copyOut());
    135     }
    136
    137     /***
    138      * Compute the sum of this and <code>m</code>.
    139      *
    140      * @param m    matrix to be added
    141      * @return     this + m
    142      * @throws  IllegalArgumentException if m is not the same size as this
    143      */
    144     public RealMatrix add(RealMatrix m) throws IllegalArgumentException {
    145         if (this.getColumnDimension() != m.getColumnDimension() ||
    146                 this.getRowDimension() != m.getRowDimension()) {
    147             throw new IllegalArgumentException("matrix dimension mismatch");
    148         }
    149         int rowCount = this.getRowDimension();
    150         int columnCount = this.getColumnDimension();
    151         double[][] outData = new double[rowCount][columnCount];
    152         for (int row = 0; row < rowCount; row++) {
    153             for (int col = 0; col < columnCount; col++) {
    154                 outData[row][col] = data[row][col] + m.getEntry(row, col);
    155             }
    156         }
    157         return new RealMatrixImpl(outData);
    158     }
    159
    160     /***
    161      * Compute  this minus <code>m</code>.
    162      *
    163      * @param m    matrix to be subtracted
    164      * @return     this + m
    165      * @throws  IllegalArgumentException if m is not the same size as *this
    166      */
    167     public RealMatrix subtract(RealMatrix m) throws IllegalArgumentException {
    168         if (this.getColumnDimension() != m.getColumnDimension() ||
    169                 this.getRowDimension() != m.getRowDimension()) {
    170             throw new IllegalArgumentException("matrix dimension mismatch");
    171         }
    172         int rowCount = this.getRowDimension();
    173         int columnCount = this.getColumnDimension();
    174         double[][] outData = new double[rowCount][columnCount];
    175         for (int row = 0; row < rowCount; row++) {
    176             for (int col = 0; col < columnCount; col++) {
    177                 outData[row][col] = data[row][col] - m.getEntry(row, col);
    178             }
    179         }
    180         return new RealMatrixImpl(outData);
    181     }
    182
    183     /***
    184      * Returns the result of adding d to each entry of this.
    185      *
    186      * @param d    value to be added to each entry
    187      * @return     d + this
    188      */
    189     public RealMatrix scalarAdd(double d) {
    190         int rowCount = this.getRowDimension();
    191         int columnCount = this.getColumnDimension();
    192         double[][] outData = new double[rowCount][columnCount];
    193         for (int row = 0; row < rowCount; row++) {
    194             for (int col = 0; col < columnCount; col++) {
    195                 outData[row][col] = data[row][col] + d;
    196             }
    197         }
    198         return new RealMatrixImpl(outData);
    199     }
    200
    201     /***
    202      * Returns the result multiplying each entry of this by <code>d</code>
    203      * @param d  value to multiply all entries by
    204      * @return d * this
    205      */
    206     public RealMatrix scalarMultiply(double d) {
    207         int rowCount = this.getRowDimension();
    208         int columnCount = this.getColumnDimension();
    209         double[][] outData = new double[rowCount][columnCount];
    210         for (int row = 0; row < rowCount; row++) {
    211             for (int col = 0; col < columnCount; col++) {
    212                 outData[row][col] = data[row][col] * d;
    213             }
    214         }
    215         return new RealMatrixImpl(outData);
    216     }
    217
    218     /***
    219      * Returns the result of postmultiplying this by <code>m</code>.
    220      * @param m    matrix to postmultiply by
    221      * @return     this*m
    222      * @throws     IllegalArgumentException
    223      *             if columnDimension(this) != rowDimension(m)
    224      */
    225     public RealMatrix multiply(RealMatrix m) throws IllegalArgumentException {
    226         if (this.getColumnDimension() != m.getRowDimension()) {
    227             throw new IllegalArgumentException("Matrices are not multiplication compatible.");
    228         }
    229         int nRows = this.getRowDimension();
    230         int nCols = m.getColumnDimension();
    231         int nSum = this.getColumnDimension();
    232         double[][] outData = new double[nRows][nCols];
    233         double sum = 0;
    234         for (int row = 0; row < nRows; row++) {
    235             for (int col = 0; col < nCols; col++) {
    236                 sum = 0;
    237                 for (int i = 0; i < nSum; i++) {
    238                     sum += data[row][i] * m.getEntry(i, col);
    239                 }
    240                 outData[row][col] = sum;
    241             }
    242         }
    243         return new RealMatrixImpl(outData);
    244     }
    245
    246     /***
    247      * Returns the result premultiplying this by <code>m</code>.
    248      * @param m    matrix to premultiply by
    249      * @return     m * this
    250      * @throws     IllegalArgumentException
    251      *             if rowDimension(this) != columnDimension(m)
    252      */
    253     public RealMatrix preMultiply(RealMatrix m) throws IllegalArgumentException {
    254         return m.multiply(this);
    255     }
    256
    257     /***
    258      * Returns matrix entries as a two-dimensional array.
    259      * <p>
    260      * Makes a fresh copy of the underlying data.
    261      *
    262      * @return    2-dimensional array of entries
    263      */
    264     public double[][] getData() {
    265         return copyOut();
    266     }
    267
    268     /***
    269      * Returns a reference to the underlying data array.
    270      * <p>
    271      * Does not make a fresh copy of the underlying data.
    272      *
    273      * @return 2-dimensional array of entries
    274      */
    275     public double[][] getDataRef() {
    276         return data;
    277     }
    278
    279     /***
    280      *
    281      * @return norm
    282      */
    283     public double getNorm() {
    284         double maxColSum = 0;
    285         for (int col = 0; col < this.getColumnDimension(); col++) {
    286             double sum = 0;
    287             for (int row = 0; row < this.getRowDimension(); row++) {
    288                 sum += Math.abs(data[row][col]);
    289             }
    290             maxColSum = Math.max(maxColSum, sum);
    291         }
    292         return maxColSum;
    293     }
    294
    295     /***
    296      * Gets a submatrix. Rows and columns are indicated
    297      * counting from 0 to n-1.
    298      *
    299      * @param startRow Initial row index
    300      * @param endRow Final row index
    301      * @param startColumn Initial column index
    302      * @param endColumn Final column index
    303      * @return The subMatrix containing the data of the
    304      *         specified rows and columns
    305      * @exception MatrixIndexException if row or column selections are not valid
    306      */
    307     public RealMatrix getSubMatrix(int startRow, int endRow, int startColumn,
    308             int endColumn) throws MatrixIndexException {
    309         if (startRow < 0 || startRow > endRow || endRow > data.length ||
    310              startColumn < 0 || startColumn > endColumn ||
    311              endColumn > data[0].length ) {
    312             throw new MatrixIndexException(
    313                     "invalid row or column index selection");
    314         }
    315         RealMatrixImpl subMatrix = new RealMatrixImpl(endRow - startRow+1,
    316                 endColumn - startColumn+1);
    317         double[][] subMatrixData = subMatrix.getDataRef();
    318         for (int i = startRow; i <= endRow; i++) {
    319             for (int j = startColumn; j <= endColumn; j++) {
    320                     subMatrixData[i - startRow][j - startColumn] = data[i][j];
    321                 }
    322             }
    323         return subMatrix;
    324     }
    325
    326     /***
    327      * Gets a submatrix. Rows and columns are indicated
    328      * counting from 0 to n-1.
    329      *
    330      * @param selectedRows Array of row indices must be non-empty
    331      * @param selectedColumns Array of column indices must be non-empty
    332      * @return The subMatrix containing the data in the
    333      *     specified rows and columns
    334      * @exception MatrixIndexException  if supplied row or column index arrays
    335      *     are not valid
    336      */
    337     public RealMatrix getSubMatrix(int[] selectedRows, int[] selectedColumns)
    338     throws MatrixIndexException {
    339         if (selectedRows.length * selectedColumns.length == 0) {
    340             throw new MatrixIndexException(
    341                     "selected row and column index arrays must be non-empty");
    342         }
    343         RealMatrixImpl subMatrix = new RealMatrixImpl(selectedRows.length,
    344                 selectedColumns.length);
    345         double[][] subMatrixData = subMatrix.getDataRef();
    346         try  {
    347             for (int i = 0; i < selectedRows.length; i++) {
    348                 for (int j = 0; j < selectedColumns.length; j++) {
    349                     subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
    350                 }
    351             }
    352         }
    353         catch (ArrayIndexOutOfBoundsException e) {
    354             throw new MatrixIndexException("matrix dimension mismatch");
    355         }
    356         return subMatrix;
    357     }
    358
    359     /***
    360      * Replace the submatrix starting at <code>row, column</code> using data in
    361      * the input <code>subMatrix</code> array. Indexes are 0-based.
    362      * <p>
    363      * Example:<br>
    364      * Starting with <pre>
    365      * 1  2  3  4
    366      * 5  6  7  8
    367      * 9  0  1  2
    368      * </pre>
    369      * and <code>subMatrix = {{3, 4} {5,6}}</code>, invoking
    370      * <code>setSubMatrix(subMatrix,1,1))</code> will result in <pre>
    371      * 1  2  3  4
    372      * 5  3  4  8
    373      * 9  5  6  2
    374      * </pre>
    375      *
    376      * @param subMatrix  array containing the submatrix replacement data
    377      * @param row  row coordinate of the top, left element to be replaced
    378      * @param column  column coordinate of the top, left element to be replaced
    379      * @throws MatrixIndexException  if subMatrix does not fit into this
    380      *    matrix from element in (row, column)
    381      * @throws IllegalArgumentException if <code>subMatrix</code> is not rectangular
    382      *  (not all rows have the same length) or empty
    383      * @throws NullPointerException if <code>subMatrix</code> is null
    384      * @since 1.1
    385      */
    386     public void setSubMatrix(double[][] subMatrix, int row, int column)
    387         throws MatrixIndexException {
    388         if ((row < 0) || (column < 0)){
    389             throw new MatrixIndexException
    390                 ("invalid row or column index selection");
    391         }
    392         int nRows = subMatrix.length;
    393         if (nRows == 0) {
    394             throw new IllegalArgumentException(
    395             "Matrix must have at least one row.");
    396         }
    397         int nCols = subMatrix[0].length;
    398         if (nCols == 0) {
    399             throw new IllegalArgumentException(
    400             "Matrix must have at least one column.");
    401         }
    402         for (int r = 1; r < nRows; r++) {
    403             if (subMatrix[r].length != nCols) {
    404                 throw new IllegalArgumentException(
    405                 "All input rows must have the same length.");
    406             }
    407         }
    408         if (data == null) {
    409             if ((row > 0)||(column > 0)) throw new MatrixIndexException
    410                 ("matrix must be initialized to perfom this method");
    411             data = new double[nRows][nCols];
    412             System.arraycopy(subMatrix, 0, data, 0, subMatrix.length);
    413         }
    414         if (((nRows + row) > this.getRowDimension()) ||
    415             (nCols + column > this.getColumnDimension()))
    416             throw new MatrixIndexException(
    417                     "invalid row or column index selection");
    418         for (int i = 0; i < nRows; i++) {
    419             System.arraycopy(subMatrix[i], 0, data[row + i], column, nCols);
    420         }
    421         lu = null;
    422     }
    423
    424     /***
    425      * Returns the entries in row number <code>row</code> as a row matrix.
    426      * Row indices start at 0.
    427      *
    428      * @param row  the row to be fetched
    429      * @return row matrix
    430      * @throws MatrixIndexException if the specified row index is invalid
    431      */
    432     public RealMatrix getRowMatrix(int row) throws MatrixIndexException {
    433         if ( !isValidCoordinate( row, 0)) {
    434             throw new MatrixIndexException("illegal row argument");
    435         }
    436         int ncols = this.getColumnDimension();
    437         double[][] out = new double[1][ncols];
    438         System.arraycopy(data[row], 0, out[0], 0, ncols);
    439         return new RealMatrixImpl(out);
    440     }
    441
    442     /***
    443      * Returns the entries in column number <code>column</code>
    444      * as a column matrix.  Column indices start at 0.
    445      *
    446      * @param column the column to be fetched
    447      * @return column matrix
    448      * @throws MatrixIndexException if the specified column index is invalid
    449      */
    450     public RealMatrix getColumnMatrix(int column) throws MatrixIndexException {
    451         if ( !isValidCoordinate( 0, column)) {
    452             throw new MatrixIndexException("illegal column argument");
    453         }
    454         int nRows = this.getRowDimension();
    455         double[][] out = new double[nRows][1];
    456         for (int row = 0; row < nRows; row++) {
    457             out[row][0] = data[row][column];
    458         }
    459         return new RealMatrixImpl(out);
    460     }
    461
    462      /***
    463      * Returns the entries in row number <code>row</code> as an array.
    464      * <p>
    465      * Row indices start at 0.  A <code>MatrixIndexException</code> is thrown
    466      * unless <code>0 <= row < rowDimension.</code>
    467      *
    468      * @param row the row to be fetched
    469      * @return array of entries in the row
    470      * @throws MatrixIndexException if the specified row index is not valid
    471      */
    472     public double[] getRow(int row) throws MatrixIndexException {
    473         if ( !isValidCoordinate( row, 0 ) ) {
    474             throw new MatrixIndexException("illegal row argument");
    475         }
    476         int ncols = this.getColumnDimension();
    477         double[] out = new double[ncols];
    478         System.arraycopy(data[row], 0, out, 0, ncols);
    479         return out;
    480     }
    481
    482     /***
    483      * Returns the entries in column number <code>col</code> as an array.
    484      * <p>
    485      * Column indices start at 0.  A <code>MatrixIndexException</code> is thrown
    486      * unless <code>0 <= column < columnDimension.</code>
    487      *
    488      * @param col the column to be fetched
    489      * @return array of entries in the column
    490      * @throws MatrixIndexException if the specified column index is not valid
    491      */
    492     public double[] getColumn(int col) throws MatrixIndexException {
    493         if ( !isValidCoordinate(0, col) ) {
    494             throw new MatrixIndexException("illegal column argument");
    495         }
    496         int nRows = this.getRowDimension();
    497         double[] out = new double[nRows];
    498         for (int row = 0; row < nRows; row++) {
    499             out[row] = data[row][col];
    500         }
    501         return out;
    502     }
    503
    504     /***
    505      * Returns the entry in the specified row and column.
    506      * <p>
    507      * Row and column indices start at 0 and must satisfy
    508      * <ul>
    509      * <li><code>0 <= row < rowDimension</code></li>
    510      * <li><code> 0 <= column < columnDimension</code></li>
    511      * </ul>
    512      * otherwise a <code>MatrixIndexException</code> is thrown.
    513      *
    514      * @param row  row location of entry to be fetched
    515      * @param column  column location of entry to be fetched
    516      * @return matrix entry in row,column
    517      * @throws MatrixIndexException if the row or column index is not valid
    518      */
    519     public double getEntry(int row, int column)
    520         throws MatrixIndexException {
    521         if (!isValidCoordinate(row,column)) {
    522             throw new MatrixIndexException("matrix entry does not exist");
    523         }
    524         return data[row][column];
    525     }
    526
    527     /***
    528      * Returns the transpose matrix.
    529      *
    530      * @return transpose matrix
    531      */
    532     public RealMatrix transpose() {
    533         int nRows = this.getRowDimension();
    534         int nCols = this.getColumnDimension();
    535         RealMatrixImpl out = new RealMatrixImpl(nCols, nRows);
    536         double[][] outData = out.getDataRef();
    537         for (int row = 0; row < nRows; row++) {
    538             for (int col = 0; col < nCols; col++) {
    539                 outData[col][row] = data[row][col];
    540             }
    541         }
    542         return out;
    543     }
    544
    545     /***
    546      * Returns the inverse matrix if this matrix is invertible.
    547      *
    548      * @return inverse matrix
    549      * @throws InvalidMatrixException if this is not invertible
    550      */
    551     public RealMatrix inverse() throws InvalidMatrixException {
    552         return solve(MatrixUtils.createRealIdentityMatrix
    553                 (this.getRowDimension()));
    554     }
    555
    556     /***
    557      * @return determinant
    558      * @throws InvalidMatrixException if matrix is not square
    559      */
    560     public double getDeterminant() throws InvalidMatrixException {
    561         if (!isSquare()) {
    562             throw new InvalidMatrixException("matrix is not square");
    563         }
    564         if (isSingular()) {   // note: this has side effect of attempting LU decomp if lu == null
    565             return 0d;
    566         } else {
    567             double det = parity;
    568             for (int i = 0; i < this.getRowDimension(); i++) {
    569                 det *= lu[i][i];
    570             }
    571             return det;
    572         }
    573     }
    574
    575     /***
    576      * @return true if the matrix is square (rowDimension = columnDimension)
    577      */
    578     public boolean isSquare() {
    579         return (this.getColumnDimension() == this.getRowDimension());
    580     }
    581
    582     /***
    583      * @return true if the matrix is singular
    584      */
    585     public boolean isSingular() {
    586         if (lu == null) {
    587             try {
    588                 luDecompose();
    589                 return false;
    590             } catch (InvalidMatrixException ex) {
    591                 return true;
    592             }
    593         } else { // LU decomp must have been successfully performed
    594             return false; // so the matrix is not singular
    595         }
    596     }
    597
    598     /***
    599      * @return rowDimension
    600      */
    601     public int getRowDimension() {
    602         return data.length;
    603     }
    604
    605     /***
    606      * @return columnDimension
    607      */
    608     public int getColumnDimension() {
    609         return data[0].length;
    610     }
    611
    612     /***
    613      * @return trace
    614      * @throws IllegalArgumentException if the matrix is not square
    615      */
    616     public double getTrace() throws IllegalArgumentException {
    617         if (!isSquare()) {
    618             throw new IllegalArgumentException("matrix is not square");
    619         }
    620         double trace = data[0][0];
    621         for (int i = 1; i < this.getRowDimension(); i++) {
    622             trace += data[i][i];
    623         }
    624         return trace;
    625     }
    626
    627     /***
    628      * @param v vector to operate on
    629      * @throws IllegalArgumentException if columnDimension != v.length
    630      * @return resulting vector
    631      */
    632     public double[] operate(double[] v) throws IllegalArgumentException {
    633         if (v.length != this.getColumnDimension()) {
    634             throw new IllegalArgumentException("vector has wrong length");
    635         }
    636         int nRows = this.getRowDimension();
    637         int nCols = this.getColumnDimension();
    638         double[] out = new double[v.length];
    639         for (int row = 0; row < nRows; row++) {
    640             double sum = 0;
    641             for (int i = 0; i < nCols; i++) {
    642                 sum += data[row][i] * v[i];
    643             }
    644             out[row] = sum;
    645         }
    646         return out;
    647     }
    648
    649     /***
    650      * @param v vector to premultiply by
    651      * @throws IllegalArgumentException if rowDimension != v.length
    652      * @return resulting matrix
    653      */
    654     public double[] preMultiply(double[] v) throws IllegalArgumentException {
    655         int nRows = this.getRowDimension();
    656         if (v.length != nRows) {
    657             throw new IllegalArgumentException("vector has wrong length");
    658         }
    659         int nCols = this.getColumnDimension();
    660         double[] out = new double[nCols];
    661         for (int col = 0; col < nCols; col++) {
    662             double sum = 0;
    663             for (int i = 0; i < nRows; i++) {
    664                 sum += data[i][col] * v[i];
    665             }
    666             out[col] = sum;
    667         }
    668         return out;
    669     }
    670
    671     /***
    672      * Returns a matrix of (column) solution vectors for linear systems with
    673      * coefficient matrix = this and constant vectors = columns of
    674      * <code>b</code>.
    675      *
    676      * @param b  array of constant forming RHS of linear systems to
    677      * to solve
    678      * @return solution array
    679      * @throws IllegalArgumentException if this.rowDimension != row dimension
    680      * @throws InvalidMatrixException if this matrix is not square or is singular
    681      */
    682     public double[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
    683         int nRows = this.getRowDimension();
    684         if (b.length != nRows) {
    685             throw new IllegalArgumentException("constant vector has wrong length");
    686         }
    687         RealMatrix bMatrix = new RealMatrixImpl(b);
    688         double[][] solution = ((RealMatrixImpl) (solve(bMatrix))).getDataRef();
    689         double[] out = new double[nRows];
    690         for (int row = 0; row < nRows; row++) {
    691             out[row] = solution[row][0];
    692         }
    693         return out;
    694     }
    695
    696     /***
    697      * Returns a matrix of (column) solution vectors for linear systems with
    698      * coefficient matrix = this and constant vectors = columns of
    699      * <code>b</code>.
    700      *
    701      * @param b  matrix of constant vectors forming RHS of linear systems to
    702      * to solve
    703      * @return matrix of solution vectors
    704      * @throws IllegalArgumentException if this.rowDimension != row dimension
    705      * @throws InvalidMatrixException if this matrix is not square or is singular
    706      */
    707     public RealMatrix solve(RealMatrix b) throws IllegalArgumentException, InvalidMatrixException  {
    708         if (b.getRowDimension() != this.getRowDimension()) {
    709             throw new IllegalArgumentException("Incorrect row dimension");
    710         }
    711         if (!this.isSquare()) {
    712             throw new InvalidMatrixException("coefficient matrix is not square");
    713         }
    714         if (this.isSingular()) { // side effect: compute LU decomp
    715             throw new InvalidMatrixException("Matrix is singular.");
    716         }
    717
    718         int nCol = this.getColumnDimension();
    719         int nColB = b.getColumnDimension();
    720         int nRowB = b.getRowDimension();
    721
    722         // Apply permutations to b
    723         double[][] bp = new double[nRowB][nColB];
    724         for (int row = 0; row < nRowB; row++) {
    725             for (int col = 0; col < nColB; col++) {
    726                 bp[row][col] = b.getEntry(permutation[row], col);
    727             }
    728         }
    729
    730         // Solve LY = b
    731         for (int col = 0; col < nCol; col++) {
    732             for (int i = col + 1; i < nCol; i++) {
    733                 for (int j = 0; j < nColB; j++) {
    734                     bp[i][j] -= bp[col][j] * lu[i][col];
    735                 }
    736             }
    737         }
    738
    739         // Solve UX = Y
    740         for (int col = nCol - 1; col >= 0; col--) {
    741             for (int j = 0; j < nColB; j++) {
    742                 bp[col][j] /= lu[col][col];
    743             }
    744             for (int i = 0; i < col; i++) {
    745                 for (int j = 0; j < nColB; j++) {
    746                     bp[i][j] -= bp[col][j] * lu[i][col];
    747                 }
    748             }
    749         }
    750
    751         RealMatrixImpl outMat = new RealMatrixImpl(bp);
    752         return outMat;
    753     }
    754
    755     /***
    756      * Computes a new
    757      * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
    758      * LU decompostion</a> for this matrix, storing the result for use by other methods.
    759      * <p>
    760      * <strong>Implementation Note</strong>:<br>
    761      * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
    762      * Crout's algortithm</a>, with partial pivoting.
    763      * <p>
    764      * <strong>Usage Note</strong>:<br>
    765      * This method should rarely be invoked directly. Its only use is
    766      * to force recomputation of the LU decomposition when changes have been
    767      * made to the underlying data using direct array references. Changes
    768      * made using setXxx methods will trigger recomputation when needed
    769      * automatically.
    770      *
    771      * @throws InvalidMatrixException if the matrix is non-square or singular.
    772      */
    773     public void luDecompose() throws InvalidMatrixException {
    774
    775         int nRows = this.getRowDimension();
    776         int nCols = this.getColumnDimension();
    777         if (nRows != nCols) {
    778             throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
    779         }
    780         lu = this.getData();
    781
    782         // Initialize permutation array and parity
    783         permutation = new int[nRows];
    784         for (int row = 0; row < nRows; row++) {
    785             permutation[row] = row;
    786         }
    787         parity = 1;
    788
    789         // Loop over columns
    790         for (int col = 0; col < nCols; col++) {
    791
    792             double sum = 0;
    793
    794             // upper
    795             for (int row = 0; row < col; row++) {
    796                 sum = lu[row][col];
    797                 for (int i = 0; i < row; i++) {
    798                     sum -= lu[row][i] * lu[i][col];
    799                 }
    800                 lu[row][col] = sum;
    801             }
    802
    803             // lower
    804             int max = col; // permutation row
    805             double largest = 0d;
    806             for (int row = col; row < nRows; row++) {
    807                 sum = lu[row][col];
    808                 for (int i = 0; i < col; i++) {
    809                     sum -= lu[row][i] * lu[i][col];
    810                 }
    811                 lu[row][col] = sum;
    812
    813                 // maintain best permutation choice
    814                 if (Math.abs(sum) > largest) {
    815                     largest = Math.abs(sum);
    816                     max = row;
    817                 }
    818             }
    819
    820             // Singularity check
    821             if (Math.abs(lu[max][col]) < TOO_SMALL) {
    822                 lu = null;
    823                 throw new InvalidMatrixException("matrix is singular");
    824             }
    825
    826             // Pivot if necessary
    827             if (max != col) {
    828                 double tmp = 0;
    829                 for (int i = 0; i < nCols; i++) {
    830                     tmp = lu[max][i];
    831                     lu[max][i] = lu[col][i];
    832                     lu[col][i] = tmp;
    833                 }
    834                 int temp = permutation[max];
    835                 permutation[max] = permutation[col];
    836                 permutation[col] = temp;
    837                 parity = -parity;
    838             }
    839
    840             //Divide the lower elements by the "winning" diagonal elt.
    841             for (int row = col + 1; row < nRows; row++) {
    842                 lu[row][col] /= lu[col][col];
    843             }
    844         }
    845     }
    846
    847     /***
    848      *
    849      * @see java.lang.Object#toString()
    850      */
    851     public String toString() {
    852         StringBuffer res = new StringBuffer();
    853         res.append("RealMatrixImpl{");
    854         if (data != null) {
    855             for (int i = 0; i < data.length; i++) {
    856                 if (i > 0)
    857                     res.append(",");
    858                 res.append("{");
    859                 for (int j = 0; j < data[0].length; j++) {
    860                     if (j > 0)
    861                         res.append(",");
    862                     res.append(data[i][j]);
    863                 }
    864                 res.append("}");
    865             }
    866         }
    867         res.append("}");
    868         return res.toString();
    869     }
    870
    871     /***
    872      * Returns true iff <code>object</code> is a
    873      * <code>RealMatrixImpl</code> instance with the same dimensions as this
    874      * and all corresponding matrix entries are equal.  Corresponding entries
    875      * are compared using {@link java.lang.Double#doubleToLongBits(double)}
    876      *
    877      * @param object the object to test equality against.
    878      * @return true if object equals this
    879      */
    880     public boolean equals(Object object) {
    881         if (object == this ) {
    882             return true;
    883         }
    884         if (object instanceof RealMatrixImpl == false) {
    885             return false;
    886         }
    887         RealMatrix m = (RealMatrix) object;
    888         int nRows = getRowDimension();
    889         int nCols = getColumnDimension();
    890         if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
    891             return false;
    892         }
    893         for (int row = 0; row < nRows; row++) {
    894             for (int col = 0; col < nCols; col++) {
    895                 if (Double.doubleToLongBits(data[row][col]) !=
    896                     Double.doubleToLongBits(m.getEntry(row, col))) {
    897                     return false;
    898                 }
    899             }
    900         }
    901         return true;
    902     }
    903
    904     /***
    905      * Computes a hashcode for the matrix.
    906      *
    907      * @return hashcode for matrix
    908      */
    909     public int hashCode() {
    910         int ret = 7;
    911         int nRows = getRowDimension();
    912         int nCols = getColumnDimension();
    913         ret = ret * 31 + nRows;
    914         ret = ret * 31 + nCols;
    915         for (int row = 0; row < nRows; row++) {
    916            for (int col = 0; col < nCols; col++) {
    917                ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) *
    918                    MathUtils.hash(data[row][col]);
    919            }
    920         }
    921         return ret;
    922     }
    923
    924     //------------------------ Protected methods
    925
    926     /***
    927      * Returns <code>dimension x dimension</code> identity matrix.
    928      *
    929      * @param dimension dimension of identity matrix to generate
    930      * @return identity matrix
    931      * @throws IllegalArgumentException  if dimension is not positive
    932      * @deprecated use {@link MatrixUtils#createRealIdentityMatrix}
    933      */
    934     protected RealMatrix getIdentity(int dimension) {
    935         return MatrixUtils.createRealIdentityMatrix(dimension);
    936     }
    937
    938     /***
    939      *  Returns the LU decomposition as a RealMatrix.
    940      *  Returns a fresh copy of the cached LU matrix if this has been computed;
    941      *  otherwise the composition is computed and cached for use by other methods.
    942      *  Since a copy is returned in either case, changes to the returned matrix do not
    943      *  affect the LU decomposition property.
    944      * <p>
    945      * The matrix returned is a compact representation of the LU decomposition.
    946      * Elements below the main diagonal correspond to entries of the "L" matrix;
    947      * elements on and above the main diagonal correspond to entries of the "U"
    948      * matrix.
    949      * <p>
    950      * Example: <pre>
    951      *
    952      *     Returned matrix                L                  U
    953      *         2  3  1                   1  0  0            2  3  1
    954      *         5  4  6                   5  1  0            0  4  6
    955      *         1  7  8                   1  7  1            0  0  8
    956      * </pre>
    957      *
    958      * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
    959      *  where permuteRows reorders the rows of the matrix to follow the order determined
    960      *  by the <a href=#getPermutation()>permutation</a> property.
    961      *
    962      * @return LU decomposition matrix
    963      * @throws InvalidMatrixException if the matrix is non-square or singular.
    964      */
    965     protected RealMatrix getLUMatrix() throws InvalidMatrixException {
    966         if (lu == null) {
    967             luDecompose();
    968         }
    969         return new RealMatrixImpl(lu);
    970     }
    971
    972     /***
    973      * Returns the permutation associated with the lu decomposition.
    974      * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
    975      * <p>
    976      * Example:
    977      * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
    978      * and current first row is last.
    979      * <p>
    980      * Returns a fresh copy of the array.
    981      *
    982      * @return the permutation
    983      */
    984     protected int[] getPermutation() {
    985         int[] out = new int[permutation.length];
    986         System.arraycopy(permutation, 0, out, 0, permutation.length);
    987         return out;
    988     }
    989
    990     //------------------------ Private methods
    991
    992     /***
    993      * Returns a fresh copy of the underlying data array.
    994      *
    995      * @return a copy of the underlying data array.
    996      */
    997     private double[][] copyOut() {
    998         int nRows = this.getRowDimension();
    999         double[][] out = new double[nRows][this.getColumnDimension()];
    1000         // can't copy 2-d array in one shot, otherwise get row references
    1001         for (int i = 0; i < nRows; i++) {
    1002             System.arraycopy(data[i], 0, out[i], 0, data[i].length);
    1003         }
    1004         return out;
    1005     }
    1006
    1007     /***
    1008      * Replaces data with a fresh copy of the input array.
    1009      * <p>
    1010      * Verifies that the input array is rectangular and non-empty
    1011      *
    1012      * @param in data to copy in
    1013      * @throws IllegalArgumentException if input array is emtpy or not
    1014      *    rectangular
    1015      * @throws NullPointerException if input array is null
    1016      */
    1017     private void copyIn(double[][] in) {
    1018         setSubMatrix(in,0,0);
    1019     }
    1020
    1021     /***
    1022      * Tests a given coordinate as being valid or invalid
    1023      *
    1024      * @param row the row index.
    1025      * @param col the column index.
    1026      * @return true if the coordinate is with the current dimensions
    1027      */
    1028     private boolean isValidCoordinate(int row, int col) {
    1029         int nRows = this.getRowDimension();
    1030         int nCols = this.getColumnDimension();
    1031
    1032         return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols -1);
    1033     }
    1034
    1035 }