View Javadoc

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