Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
RealMatrixImpl |
|
| 4.0;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 | } |