Coverage Report - org.apache.commons.math.linear.RealMatrixImpl

Classes in this File Line Coverage Branch Coverage Complexity
RealMatrixImpl
100% 
100% 
4

 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 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$ $Date: 2005-07-04 15:16:48 -0700 (Mon, 04 Jul 2005) $
 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  436
     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  436
     private double lu[][] = null;
 62  
 
 63  
     /** Permutation associated with LU decompostion */
 64  436
     private int[] permutation = null;
 65  
 
 66  
     /** Parity of the permutation associated with the LU decomposition */
 67  436
     private int parity = 1;
 68  
 
 69  
     /** Bound to determine effective singularity in LU decomposition */
 70  4
     protected static double TOO_SMALL = 10E-12;
 71  
 
 72  
     /**
 73  
      * Creates a matrix with no data
 74  
      */
 75  4
     public RealMatrixImpl() {
 76  4
     }
 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  54
     public RealMatrixImpl(int rowDimension, int columnDimension) {
 87  54
         if (rowDimension <= 0 || columnDimension <= 0) {
 88  4
             throw new IllegalArgumentException(
 89  
                     "row and column dimensions must be postive");
 90  
         }
 91  50
         data = new double[rowDimension][columnDimension];
 92  50
         lu = null;
 93  50
     }
 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  374
     public RealMatrixImpl(double[][] d) {
 107  374
         this.copyIn(d);
 108  364
         lu = null;
 109  364
     }
 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  4
     public RealMatrixImpl(double[] v) {
 121  4
         int nRows = v.length;
 122  4
         data = new double[nRows][1];
 123  16
         for (int row = 0; row < nRows; row++) {
 124  12
             data[row][0] = v[row];
 125  
         }
 126  4
     }
 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  2
         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  6
         if (this.getColumnDimension() != m.getColumnDimension() ||
 146  
                 this.getRowDimension() != m.getRowDimension()) {
 147  2
             throw new IllegalArgumentException("matrix dimension mismatch");
 148  
         }
 149  4
         int rowCount = this.getRowDimension();
 150  4
         int columnCount = this.getColumnDimension();
 151  4
         double[][] outData = new double[rowCount][columnCount];
 152  16
         for (int row = 0; row < rowCount; row++) {
 153  48
             for (int col = 0; col < columnCount; col++) {
 154  36
                 outData[row][col] = data[row][col] + m.getEntry(row, col);
 155  
             }  
 156  
         }
 157  4
         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  54
         if (this.getColumnDimension() != m.getColumnDimension() ||
 169  
                 this.getRowDimension() != m.getRowDimension()) {
 170  2
             throw new IllegalArgumentException("matrix dimension mismatch");
 171  
         }
 172  52
         int rowCount = this.getRowDimension();
 173  52
         int columnCount = this.getColumnDimension();
 174  52
         double[][] outData = new double[rowCount][columnCount];
 175  202
         for (int row = 0; row < rowCount; row++) {
 176  578
             for (int col = 0; col < columnCount; col++) {
 177  428
                 outData[row][col] = data[row][col] - m.getEntry(row, col);
 178  
             }
 179  
         }
 180  52
         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  2
         int rowCount = this.getRowDimension();
 191  2
         int columnCount = this.getColumnDimension();
 192  2
         double[][] outData = new double[rowCount][columnCount];
 193  8
         for (int row = 0; row < rowCount; row++) {
 194  24
             for (int col = 0; col < columnCount; col++) {
 195  18
                 outData[row][col] = data[row][col] + d;
 196  
             }
 197  
         }
 198  2
         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  2
         int rowCount = this.getRowDimension();
 208  2
         int columnCount = this.getColumnDimension();
 209  2
         double[][] outData = new double[rowCount][columnCount];
 210  8
         for (int row = 0; row < rowCount; row++) {
 211  24
             for (int col = 0; col < columnCount; col++) {
 212  18
                 outData[row][col] = data[row][col] * d;
 213  
             }
 214  
         }
 215  2
         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  38
         if (this.getColumnDimension() != m.getRowDimension()) {
 227  4
             throw new IllegalArgumentException("Matrices are not multiplication compatible.");
 228  
         }
 229  34
         int nRows = this.getRowDimension();
 230  34
         int nCols = m.getColumnDimension();
 231  34
         int nSum = this.getColumnDimension();
 232  34
         double[][] outData = new double[nRows][nCols];
 233  34
         double sum = 0;
 234  128
         for (int row = 0; row < nRows; row++) {
 235  356
             for (int col = 0; col < nCols; col++) {
 236  262
                 sum = 0;
 237  1056
                 for (int i = 0; i < nSum; i++) {
 238  794
                     sum += data[row][i] * m.getEntry(i, col);
 239  
                 }
 240  262
                 outData[row][col] = sum;
 241  
             }
 242  
         }
 243  34
         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  12
         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  48
         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  54
         return data;
 277  
     }
 278  
 
 279  
     /**
 280  
      *
 281  
      * @return norm
 282  
      */
 283  
     public double getNorm() {
 284  54
         double maxColSum = 0;
 285  206
         for (int col = 0; col < this.getColumnDimension(); col++) {
 286  152
             double sum = 0;
 287  592
             for (int row = 0; row < this.getRowDimension(); row++) {
 288  440
                 sum += Math.abs(data[row][col]);
 289  
             }
 290  152
             maxColSum = Math.max(maxColSum, sum);
 291  
         }
 292  54
         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  14
         if (startRow < 0 || startRow > endRow || endRow > data.length ||
 310  
              startColumn < 0 || startColumn > endColumn ||
 311  
              endColumn > data[0].length ) {
 312  8
             throw new MatrixIndexException(
 313  
                     "invalid row or column index selection");
 314  
         }
 315  6
         RealMatrixImpl subMatrix = new RealMatrixImpl(endRow - startRow+1,
 316  
                 endColumn - startColumn+1);
 317  6
         double[][] subMatrixData = subMatrix.getDataRef();
 318  16
         for (int i = startRow; i <= endRow; i++) {
 319  24
             for (int j = startColumn; j <= endColumn; j++) {
 320  14
                     subMatrixData[i - startRow][j - startColumn] = data[i][j];
 321  
                 }
 322  
             }
 323  6
         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  16
         if (selectedRows.length * selectedColumns.length == 0) {
 340  2
             throw new MatrixIndexException(
 341  
                     "selected row and column index arrays must be non-empty");
 342  
         }
 343  14
         RealMatrixImpl subMatrix = new RealMatrixImpl(selectedRows.length,
 344  
                 selectedColumns.length);
 345  14
         double[][] subMatrixData = subMatrix.getDataRef();
 346  
         try  {
 347  38
             for (int i = 0; i < selectedRows.length; i++) {
 348  82
                 for (int j = 0; j < selectedColumns.length; j++) {
 349  58
                     subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
 350  
                 }
 351  
             }
 352  
         }
 353  2
         catch (ArrayIndexOutOfBoundsException e) {
 354  2
             throw new MatrixIndexException("matrix dimension mismatch");
 355  12
         }
 356  12
         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  398
         if ((row < 0) || (column < 0)){
 389  4
             throw new MatrixIndexException
 390  
                 ("invalid row or column index selection");          
 391  
         }
 392  394
         int nRows = subMatrix.length;
 393  390
         if (nRows == 0) {
 394  2
             throw new IllegalArgumentException(
 395  
             "Matrix must have at least one row."); 
 396  
         }
 397  388
         int nCols = subMatrix[0].length;
 398  388
         if (nCols == 0) {
 399  6
             throw new IllegalArgumentException(
 400  
             "Matrix must have at least one column."); 
 401  
         }
 402  1068
         for (int r = 1; r < nRows; r++) {
 403  690
             if (subMatrix[r].length != nCols) {
 404  4
                 throw new IllegalArgumentException(
 405  
                 "All input rows must have the same length.");
 406  
             }
 407  
         }       
 408  378
         if (data == null) {
 409  368
             if ((row > 0)||(column > 0)) throw new MatrixIndexException
 410  
                 ("matrix must be initialized to perfom this method");
 411  364
             data = new double[nRows][nCols];
 412  364
             System.arraycopy(subMatrix, 0, data, 0, subMatrix.length);          
 413  
         }   
 414  374
         if (((nRows + row) > this.getRowDimension()) ||
 415  
             (nCols + column > this.getColumnDimension()))
 416  2
             throw new MatrixIndexException(
 417  
                     "invalid row or column index selection");                   
 418  1418
         for (int i = 0; i < nRows; i++) {
 419  1046
             System.arraycopy(subMatrix[i], 0, data[row + i], column, nCols);
 420  
         } 
 421  372
         lu = null;
 422  372
     }
 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  8
         if ( !isValidCoordinate( row, 0)) {
 434  4
             throw new MatrixIndexException("illegal row argument");
 435  
         }
 436  4
         int ncols = this.getColumnDimension();
 437  4
         double[][] out = new double[1][ncols]; 
 438  4
         System.arraycopy(data[row], 0, out[0], 0, ncols);
 439  4
         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  8
         if ( !isValidCoordinate( 0, column)) {
 452  4
             throw new MatrixIndexException("illegal column argument");
 453  
         }
 454  4
         int nRows = this.getRowDimension();
 455  4
         double[][] out = new double[nRows][1]; 
 456  20
         for (int row = 0; row < nRows; row++) {
 457  16
             out[row][0] = data[row][column];
 458  
         }
 459  4
         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  4
         if ( !isValidCoordinate( row, 0 ) ) {
 474  2
             throw new MatrixIndexException("illegal row argument");
 475  
         }
 476  2
         int ncols = this.getColumnDimension();
 477  2
         double[] out = new double[ncols];
 478  2
         System.arraycopy(data[row], 0, out, 0, ncols);
 479  2
         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  4
         if ( !isValidCoordinate(0, col) ) {
 494  2
             throw new MatrixIndexException("illegal column argument");
 495  
         }
 496  2
         int nRows = this.getRowDimension();
 497  2
         double[] out = new double[nRows];
 498  8
         for (int row = 0; row < nRows; row++) {
 499  6
             out[row] = data[row][col];
 500  
         }
 501  2
         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  1832
         if (!isValidCoordinate(row,column)) {
 522  2
             throw new MatrixIndexException("matrix entry does not exist");
 523  
         }
 524  1830
         return data[row][column];
 525  
     }
 526  
 
 527  
     /**
 528  
      * Returns the transpose matrix.
 529  
      *
 530  
      * @return transpose matrix
 531  
      */
 532  
     public RealMatrix transpose() {
 533  8
         int nRows = this.getRowDimension();
 534  8
         int nCols = this.getColumnDimension();
 535  8
         RealMatrixImpl out = new RealMatrixImpl(nCols, nRows);
 536  8
         double[][] outData = out.getDataRef();
 537  30
         for (int row = 0; row < nRows; row++) {
 538  88
             for (int col = 0; col < nCols; col++) {
 539  66
                 outData[col][row] = data[row][col];
 540  
             }
 541  
         }
 542  8
         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  16
         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  10
         if (!isSquare()) {
 562  2
             throw new InvalidMatrixException("matrix is not square");
 563  
         }
 564  8
         if (isSingular()) {   // note: this has side effect of attempting LU decomp if lu == null
 565  2
             return 0d;
 566  
         } else {
 567  6
             double det = parity;
 568  22
             for (int i = 0; i < this.getRowDimension(); i++) {
 569  16
                 det *= lu[i][i];
 570  
             }
 571  6
             return det;
 572  
         }
 573  
     }
 574  
 
 575  
     /**
 576  
      * @return true if the matrix is square (rowDimension = columnDimension)
 577  
      */
 578  
     public boolean isSquare() {
 579  60
         return (this.getColumnDimension() == this.getRowDimension());
 580  
     }
 581  
 
 582  
     /**
 583  
      * @return true if the matrix is singular
 584  
      */
 585  
     public boolean isSingular() {
 586  36
         if (lu == null) {
 587  
             try {
 588  34
                 luDecompose();
 589  24
                 return false;
 590  10
             } catch (InvalidMatrixException ex) {
 591  10
                 return true;
 592  
             }
 593  
         } else { // LU decomp must have been successfully performed
 594  2
             return false; // so the matrix is not singular
 595  
         }
 596  
     }
 597  
 
 598  
     /**
 599  
      * @return rowDimension
 600  
      */
 601  
     public int getRowDimension() {
 602  3558
         return data.length;
 603  
     }
 604  
 
 605  
     /**
 606  
      * @return columnDimension
 607  
      */
 608  
     public int getColumnDimension() {
 609  3122
         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  4
         if (!isSquare()) {
 618  2
             throw new IllegalArgumentException("matrix is not square");
 619  
         }
 620  2
         double trace = data[0][0];
 621  6
         for (int i = 1; i < this.getRowDimension(); i++) {
 622  4
             trace += data[i][i];
 623  
         }
 624  2
         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  6
         if (v.length != this.getColumnDimension()) {
 634  2
             throw new IllegalArgumentException("vector has wrong length");
 635  
         }
 636  4
         int nRows = this.getRowDimension();
 637  4
         int nCols = this.getColumnDimension();
 638  4
         double[] out = new double[v.length];
 639  16
         for (int row = 0; row < nRows; row++) {
 640  12
             double sum = 0;
 641  48
             for (int i = 0; i < nCols; i++) {
 642  36
                 sum += data[row][i] * v[i];
 643  
             }
 644  12
             out[row] = sum;
 645  
         }
 646  4
         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  4
         int nRows = this.getRowDimension();
 656  4
         if (v.length != nRows) {
 657  2
             throw new IllegalArgumentException("vector has wrong length");
 658  
         }
 659  2
         int nCols = this.getColumnDimension();
 660  2
         double[] out = new double[nCols];
 661  8
         for (int col = 0; col < nCols; col++) {
 662  6
             double sum = 0;
 663  24
             for (int i = 0; i < nRows; i++) {
 664  18
                 sum += data[i][col] * v[i];
 665  
             }
 666  6
             out[col] = sum;
 667  
         }
 668  2
         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  6
         int nRows = this.getRowDimension();
 684  6
         if (b.length != nRows) {
 685  2
             throw new IllegalArgumentException("constant vector has wrong length");
 686  
         }
 687  4
         RealMatrix bMatrix = new RealMatrixImpl(b);
 688  4
         double[][] solution = ((RealMatrixImpl) (solve(bMatrix))).getDataRef();
 689  4
         double[] out = new double[nRows];
 690  16
         for (int row = 0; row < nRows; row++) {
 691  12
             out[row] = solution[row][0];
 692  
         }
 693  4
         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  26
         if (b.getRowDimension() != this.getRowDimension()) {
 709  4
             throw new IllegalArgumentException("Incorrect row dimension");
 710  
         }
 711  22
         if (!this.isSquare()) {
 712  2
             throw new InvalidMatrixException("coefficient matrix is not square");
 713  
         }
 714  20
         if (this.isSingular()) { // side effect: compute LU decomp
 715  4
             throw new InvalidMatrixException("Matrix is singular.");
 716  
         }
 717  
 
 718  16
         int nCol = this.getColumnDimension();
 719  16
         int nColB = b.getColumnDimension();
 720  16
         int nRowB = b.getRowDimension();
 721  
 
 722  
         // Apply permutations to b
 723  16
         double[][] bp = new double[nRowB][nColB];
 724  62
         for (int row = 0; row < nRowB; row++) {
 725  156
             for (int col = 0; col < nColB; col++) {
 726  110
                 bp[row][col] = b.getEntry(permutation[row], col);
 727  
             }
 728  
         }
 729  
 
 730  
         // Solve LY = b
 731  62
         for (int col = 0; col < nCol; col++) {
 732  90
             for (int i = col + 1; i < nCol; i++) {
 733  150
                 for (int j = 0; j < nColB; j++) {
 734  106
                     bp[i][j] -= bp[col][j] * lu[i][col];
 735  
                 }
 736  
             }
 737  
         }
 738  
 
 739  
         // Solve UX = Y
 740  62
         for (int col = nCol - 1; col >= 0; col--) {
 741  156
             for (int j = 0; j < nColB; j++) {
 742  110
                 bp[col][j] /= lu[col][col];
 743  
             }
 744  90
             for (int i = 0; i < col; i++) {
 745  150
                 for (int j = 0; j < nColB; j++) {
 746  106
                     bp[i][j] -= bp[col][j] * lu[i][col];
 747  
                 }
 748  
             }
 749  
         }
 750  
 
 751  16
         RealMatrixImpl outMat = new RealMatrixImpl(bp);
 752  16
         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  48
         int nRows = this.getRowDimension();
 776  48
         int nCols = this.getColumnDimension();
 777  48
         if (nRows != nCols) {
 778  4
             throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
 779  
         }
 780  44
         lu = this.getData();
 781  
 
 782  
         // Initialize permutation array and parity
 783  44
         permutation = new int[nRows];
 784  176
         for (int row = 0; row < nRows; row++) {
 785  132
             permutation[row] = row;
 786  
         }
 787  44
         parity = 1;
 788  
 
 789  
         // Loop over columns
 790  164
         for (int col = 0; col < nCols; col++) {
 791  
 
 792  132
             double sum = 0;
 793  
 
 794  
             // upper
 795  272
             for (int row = 0; row < col; row++) {
 796  140
                 sum = lu[row][col];
 797  200
                 for (int i = 0; i < row; i++) {
 798  60
                     sum -= lu[row][i] * lu[i][col];
 799  
                 }
 800  140
                 lu[row][col] = sum;
 801  
             }
 802  
 
 803  
             // lower
 804  132
             int max = col; // permutation row
 805  132
             double largest = 0d;
 806  404
             for (int row = col; row < nRows; row++) {
 807  272
                 sum = lu[row][col];
 808  472
                 for (int i = 0; i < col; i++) {
 809  200
                     sum -= lu[row][i] * lu[i][col];
 810  
                 }
 811  272
                 lu[row][col] = sum;
 812  
 
 813  
                 // maintain best permutation choice
 814  272
                 if (Math.abs(sum) > largest) {
 815  206
                     largest = Math.abs(sum);
 816  206
                     max = row;
 817  
                 }
 818  
             }
 819  
 
 820  
             // Singularity check
 821  132
             if (Math.abs(lu[max][col]) < TOO_SMALL) {
 822  12
                 lu = null;
 823  12
                 throw new InvalidMatrixException("matrix is singular");
 824  
             }
 825  
 
 826  
             // Pivot if necessary
 827  120
             if (max != col) {
 828  66
                 double tmp = 0;
 829  284
                 for (int i = 0; i < nCols; i++) {
 830  218
                     tmp = lu[max][i];
 831  218
                     lu[max][i] = lu[col][i];
 832  218
                     lu[col][i] = tmp;
 833  
                 }
 834  66
                 int temp = permutation[max];
 835  66
                 permutation[max] = permutation[col];
 836  66
                 permutation[col] = temp;
 837  66
                 parity = -parity;
 838  
             }
 839  
 
 840  
             //Divide the lower elements by the "winning" diagonal elt.
 841  260
             for (int row = col + 1; row < nRows; row++) {
 842  140
                 lu[row][col] /= lu[col][col];
 843  
             }
 844  
         }
 845  32
     }
 846  
 
 847  
     /**
 848  
      *
 849  
      * @see java.lang.Object#toString()
 850  
      */
 851  
     public String toString() {
 852  4
         StringBuffer res = new StringBuffer();
 853  4
         res.append("RealMatrixImpl{");
 854  4
         if (data != null) {
 855  8
             for (int i = 0; i < data.length; i++) {
 856  6
                 if (i > 0)
 857  4
                     res.append(",");
 858  6
                 res.append("{");
 859  24
                 for (int j = 0; j < data[0].length; j++) {
 860  18
                     if (j > 0)
 861  12
                         res.append(",");
 862  18
                     res.append(data[i][j]);
 863  
                 } 
 864  6
                 res.append("}");
 865  
             } 
 866  
         }
 867  4
         res.append("}");
 868  4
         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.
 875  
      * 
 876  
      * @param object the object to test equality against.
 877  
      * @return true if object equals this
 878  
      */
 879  
     public boolean equals(Object object) {
 880  52
         if (object == this ) {
 881  2
             return true;
 882  
         }
 883  50
         if (object instanceof RealMatrixImpl == false) {
 884  2
             return false;
 885  
         }
 886  48
         RealMatrix m = (RealMatrix) object;
 887  48
         int nRows = getRowDimension();
 888  48
         int nCols = getColumnDimension();
 889  48
         if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
 890  2
             return false;
 891  
         }
 892  150
         for (int row = 0; row < nRows; row++) {
 893  356
             for (int col = 0; col < nCols; col++) {
 894  252
                 if (data[row][col] != m.getEntry(row, col)) {
 895  2
                     return false;
 896  
                 }
 897  
             }
 898  
         }
 899  44
         return true;
 900  
     }
 901  
     
 902  
     /**
 903  
      * Computes a hashcode for the matrix.
 904  
      * 
 905  
      * @return hashcode for matrix
 906  
      */
 907  
     public int hashCode() {
 908  8
         int ret = 7;
 909  8
         int nRows = getRowDimension();
 910  8
         int nCols = getColumnDimension();
 911  8
         ret = ret * 31 + nRows;
 912  8
         ret = ret * 31 + nCols;
 913  32
         for (int row = 0; row < nRows; row++) {
 914  96
            for (int col = 0; col < nCols; col++) {
 915  72
                ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) * 
 916  
                    MathUtils.hash(data[row][col]);
 917  
            }
 918  
         }
 919  8
         return ret;
 920  
     }
 921  
 
 922  
     //------------------------ Protected methods
 923  
 
 924  
     /**
 925  
      * Returns <code>dimension x dimension</code> identity matrix.
 926  
      *
 927  
      * @param dimension dimension of identity matrix to generate
 928  
      * @return identity matrix
 929  
      * @throws IllegalArgumentException  if dimension is not positive
 930  
      * @deprecated use {@link MatrixUtils#createRealIdentityMatrix}
 931  
      */
 932  
     protected RealMatrix getIdentity(int dimension) {
 933  0
         return MatrixUtils.createRealIdentityMatrix(dimension);
 934  
     }
 935  
 
 936  
     /**
 937  
      *  Returns the LU decomposition as a RealMatrix.
 938  
      *  Returns a fresh copy of the cached LU matrix if this has been computed;
 939  
      *  otherwise the composition is computed and cached for use by other methods.
 940  
      *  Since a copy is returned in either case, changes to the returned matrix do not
 941  
      *  affect the LU decomposition property.
 942  
      * <p>
 943  
      * The matrix returned is a compact representation of the LU decomposition.
 944  
      * Elements below the main diagonal correspond to entries of the "L" matrix;
 945  
      * elements on and above the main diagonal correspond to entries of the "U"
 946  
      * matrix.
 947  
      * <p>
 948  
      * Example: <pre>
 949  
      *
 950  
      *     Returned matrix                L                  U
 951  
      *         2  3  1                   1  0  0            2  3  1
 952  
      *         5  4  6                   5  1  0            0  4  6
 953  
      *         1  7  8                   1  7  1            0  0  8
 954  
      * </pre>
 955  
      *
 956  
      * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
 957  
      *  where permuteRows reorders the rows of the matrix to follow the order determined
 958  
      *  by the <a href=#getPermutation()>permutation</a> property.
 959  
      *
 960  
      * @return LU decomposition matrix
 961  
      * @throws InvalidMatrixException if the matrix is non-square or singular.
 962  
      */
 963  
     protected RealMatrix getLUMatrix() throws InvalidMatrixException {
 964  14
         if (lu == null) {
 965  12
             luDecompose();
 966  
         }
 967  10
         return new RealMatrixImpl(lu);
 968  
     }
 969  
 
 970  
     /**
 971  
      * Returns the permutation associated with the lu decomposition.
 972  
      * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
 973  
      * <p>
 974  
      * Example:
 975  
      * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
 976  
      * and current first row is last.
 977  
      * <p>
 978  
      * Returns a fresh copy of the array.
 979  
      *
 980  
      * @return the permutation
 981  
      */
 982  
     protected int[] getPermutation() {
 983  10
         int[] out = new int[permutation.length];
 984  10
         System.arraycopy(permutation, 0, out, 0, permutation.length);
 985  10
         return out;
 986  
     }
 987  
 
 988  
     //------------------------ Private methods
 989  
 
 990  
     /**
 991  
      * Returns a fresh copy of the underlying data array.
 992  
      *
 993  
      * @return a copy of the underlying data array.
 994  
      */
 995  
     private double[][] copyOut() {
 996  50
         int nRows = this.getRowDimension();
 997  50
         double[][] out = new double[nRows][this.getColumnDimension()];
 998  
         // can't copy 2-d array in one shot, otherwise get row references
 999  200
         for (int i = 0; i < nRows; i++) {
 1000  150
             System.arraycopy(data[i], 0, out[i], 0, data[i].length);
 1001  
         }
 1002  50
         return out;
 1003  
     }
 1004  
 
 1005  
     /**
 1006  
      * Replaces data with a fresh copy of the input array.
 1007  
      * <p>
 1008  
      * Verifies that the input array is rectangular and non-empty
 1009  
      *
 1010  
      * @param in data to copy in
 1011  
      * @throws IllegalArgumentException if input array is emtpy or not
 1012  
      *    rectangular
 1013  
      * @throws NullPointerException if input array is null
 1014  
      */
 1015  
     private void copyIn(double[][] in) {
 1016  374
         setSubMatrix(in,0,0);
 1017  364
     }
 1018  
 
 1019  
     /**
 1020  
      * Tests a given coordinate as being valid or invalid
 1021  
      *
 1022  
      * @param row the row index.
 1023  
      * @param col the column index.
 1024  
      * @return true if the coordinate is with the current dimensions
 1025  
      */
 1026  
     private boolean isValidCoordinate(int row, int col) {
 1027  1856
         int nRows = this.getRowDimension();
 1028  1856
         int nCols = this.getColumnDimension();
 1029  
 
 1030  1856
         return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols -1);
 1031  
     }
 1032  
 
 1033  
 }