Java > Open Source Codes > org > apache > commons > math > linear > RealMatrixImpl


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     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.
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         if (object == this ) {
881             return true;
882         }
883         if (object instanceof RealMatrixImpl == false) {
884             return false;
885         }
886         RealMatrix m = (RealMatrix) object;
887         int nRows = getRowDimension();
888         int nCols = getColumnDimension();
889         if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
890             return false;
891         }
892         for (int row = 0; row < nRows; row++) {
893             for (int col = 0; col < nCols; col++) {
894                 if (data[row][col] != m.getEntry(row, col)) {
895                     return false;
896                 }
897             }
898         }
899         return true;
900     }
901     
902     /**
903      * Computes a hashcode for the matrix.
904      *
905      * @return hashcode for matrix
906      */
907     public int hashCode() {
908         int ret = 7;
909         int nRows = getRowDimension();
910         int nCols = getColumnDimension();
911         ret = ret * 31 + nRows;
912         ret = ret * 31 + nCols;
913         for (int row = 0; row < nRows; row++) {
914            for (int col = 0; col < nCols; col++) {
915                ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) *
916                    MathUtils.hash(data[row][col]);
917            }
918         }
919         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         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         if (lu == null) {
965             luDecompose();
966         }
967         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         int[] out = new int[permutation.length];
984         System.arraycopy(permutation, 0, out, 0, permutation.length);
985         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         int nRows = this.getRowDimension();
997         double[][] out = new double[nRows][this.getColumnDimension()];
998         // can't copy 2-d array in one shot, otherwise get row references
999 for (int i = 0; i < nRows; i++) {
1000            System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1001        }
1002        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        setSubMatrix(in,0,0);
1017    }
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        int nRows = this.getRowDimension();
1028        int nCols = this.getColumnDimension();
1029
1030        return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols -1);
1031    }
1032
1033}
1034

A to Z: JavaDoc & Examples

---

Daily Java News & Articles

---

Open Source Projects

---

Open Source Codes

---

Free Computer Books

---

Remove Frame

---

Popular Tags