Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
BigMatrixImpl |
|
| 3.7586206896551726;3.759 |
1 | /* |
|
2 | * Copyright 2004 The Apache Software Foundation. |
|
3 | * |
|
4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
|
5 | * you may not use this file except in compliance with the License. |
|
6 | * You may obtain a copy of the License at |
|
7 | * |
|
8 | * http://www.apache.org/licenses/LICENSE-2.0 |
|
9 | * |
|
10 | * Unless required by applicable law or agreed to in writing, software |
|
11 | * distributed under the License is distributed on an "AS IS" BASIS, |
|
12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
13 | * See the License for the specific language governing permissions and |
|
14 | * limitations under the License. |
|
15 | */ |
|
16 | ||
17 | package org.apache.commons.math.linear; |
|
18 | import java.io.Serializable; |
|
19 | import java.math.BigDecimal; |
|
20 | ||
21 | /** |
|
22 | * Implementation for {@link BigMatrix} using a BigDecimal[][] array to store entries |
|
23 | * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf"> |
|
24 | * LU decompostion</a> to support linear system |
|
25 | * solution and inverse. |
|
26 | * <p> |
|
27 | * The LU decompostion is performed as needed, to support the following operations: <ul> |
|
28 | * <li>solve</li> |
|
29 | * <li>isSingular</li> |
|
30 | * <li>getDeterminant</li> |
|
31 | * <li>inverse</li> </ul> |
|
32 | * <p> |
|
33 | * <strong>Usage notes</strong>:<br> |
|
34 | * <ul><li> |
|
35 | * The LU decomposition is stored and reused on subsequent calls. If matrix |
|
36 | * data are modified using any of the public setXxx methods, the saved |
|
37 | * decomposition is discarded. If data are modified via references to the |
|
38 | * underlying array obtained using <code>getDataRef()</code>, then the stored |
|
39 | * LU decomposition will not be discarded. In this case, you need to |
|
40 | * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition |
|
41 | * before using any of the methods above.</li> |
|
42 | * <li> |
|
43 | * As specified in the {@link BigMatrix} interface, matrix element indexing |
|
44 | * is 0-based -- e.g., <code>getEntry(0, 0)</code> |
|
45 | * returns the element in the first row, first column of the matrix.</li></ul> |
|
46 | * @version $Revision$ $Date: 2005-07-04 15:16:48 -0700 (Mon, 04 Jul 2005) $ |
|
47 | */ |
|
48 | public class BigMatrixImpl implements BigMatrix, Serializable { |
|
49 | ||
50 | /** Serialization id */ |
|
51 | static final long serialVersionUID = -1011428905656140431L; |
|
52 | ||
53 | /** Entries of the matrix */ |
|
54 | 452 | private BigDecimal data[][] = null; |
55 | ||
56 | /** Entries of cached LU decomposition. |
|
57 | * All updates to data (other than luDecompose()) *must* set this to null |
|
58 | */ |
|
59 | 452 | private BigDecimal lu[][] = null; |
60 | ||
61 | /** Permutation associated with LU decompostion */ |
|
62 | 452 | private int[] permutation = null; |
63 | ||
64 | /** Parity of the permutation associated with the LU decomposition */ |
|
65 | 452 | private int parity = 1; |
66 | ||
67 | /** Rounding mode for divisions **/ |
|
68 | 452 | private int roundingMode = BigDecimal.ROUND_HALF_UP; |
69 | ||
70 | /*** BigDecimal scale ***/ |
|
71 | 452 | private int scale = 64; |
72 | ||
73 | /** Bound to determine effective singularity in LU decomposition */ |
|
74 | 4 | protected static BigDecimal TOO_SMALL = new BigDecimal(10E-12); |
75 | ||
76 | /** BigDecimal 0 */ |
|
77 | 4 | static final BigDecimal ZERO = new BigDecimal(0); |
78 | /** BigDecimal 1 */ |
|
79 | 4 | static final BigDecimal ONE = new BigDecimal(1); |
80 | ||
81 | /** |
|
82 | * Creates a matrix with no data |
|
83 | */ |
|
84 | 2 | public BigMatrixImpl() { |
85 | 2 | } |
86 | ||
87 | /** |
|
88 | * Create a new BigMatrix with the supplied row and column dimensions. |
|
89 | * |
|
90 | * @param rowDimension the number of rows in the new matrix |
|
91 | * @param columnDimension the number of columns in the new matrix |
|
92 | * @throws IllegalArgumentException if row or column dimension is not |
|
93 | * positive |
|
94 | */ |
|
95 | 52 | public BigMatrixImpl(int rowDimension, int columnDimension) { |
96 | 52 | if (rowDimension <=0 || columnDimension <=0) { |
97 | 4 | throw new IllegalArgumentException |
98 | ("row and column dimensions must be positive"); |
|
99 | } |
|
100 | 48 | data = new BigDecimal[rowDimension][columnDimension]; |
101 | 48 | lu = null; |
102 | 48 | } |
103 | ||
104 | /** |
|
105 | * Create a new BigMatrix using the <code>data</code> as the underlying |
|
106 | * data array. |
|
107 | * <p> |
|
108 | * The input array is copied, not referenced. |
|
109 | * |
|
110 | * @param d data for new matrix |
|
111 | * @throws IllegalArgumentException if <code>d</code> is not rectangular |
|
112 | * (not all rows have the same length) or empty |
|
113 | * @throws NullPointerException if <code>d</code> is null |
|
114 | */ |
|
115 | 164 | public BigMatrixImpl(BigDecimal[][] d) { |
116 | 164 | this.copyIn(d); |
117 | 164 | lu = null; |
118 | 164 | } |
119 | ||
120 | /** |
|
121 | * Create a new BigMatrix using the <code>data</code> as the underlying |
|
122 | * data array. |
|
123 | * <p> |
|
124 | * The input array is copied, not referenced. |
|
125 | * |
|
126 | * @param d data for new matrix |
|
127 | * @throws IllegalArgumentException if <code>d</code> is not rectangular |
|
128 | * (not all rows have the same length) or empty |
|
129 | * @throws NullPointerException if <code>d</code> is null |
|
130 | */ |
|
131 | 206 | public BigMatrixImpl(double[][] d) { |
132 | 206 | int nRows = d.length; |
133 | 204 | if (nRows == 0) { |
134 | 2 | throw new IllegalArgumentException( |
135 | "Matrix must have at least one row."); |
|
136 | } |
|
137 | 202 | int nCols = d[0].length; |
138 | 202 | if (nCols == 0) { |
139 | 4 | throw new IllegalArgumentException( |
140 | "Matrix must have at least one column."); |
|
141 | } |
|
142 | 546 | for (int row = 1; row < nRows; row++) { |
143 | 350 | if (d[row].length != nCols) { |
144 | 2 | throw new IllegalArgumentException( |
145 | "All input rows must have the same length."); |
|
146 | } |
|
147 | } |
|
148 | 196 | this.copyIn(d); |
149 | 196 | lu = null; |
150 | 196 | } |
151 | ||
152 | /** |
|
153 | * Create a new BigMatrix using the values represented by the strings in |
|
154 | * <code>data</code> as the underlying data array. |
|
155 | * |
|
156 | * @param d data for new matrix |
|
157 | * @throws IllegalArgumentException if <code>d</code> is not rectangular |
|
158 | * (not all rows have the same length) or empty |
|
159 | * @throws NullPointerException if <code>d</code> is null |
|
160 | */ |
|
161 | 26 | public BigMatrixImpl(String[][] d) { |
162 | 26 | int nRows = d.length; |
163 | 26 | if (nRows == 0) { |
164 | 2 | throw new IllegalArgumentException( |
165 | "Matrix must have at least one row."); |
|
166 | } |
|
167 | 24 | int nCols = d[0].length; |
168 | 24 | if (nCols == 0) { |
169 | 2 | throw new IllegalArgumentException( |
170 | "Matrix must have at least one column."); |
|
171 | } |
|
172 | 42 | for (int row = 1; row < nRows; row++) { |
173 | 22 | if (d[row].length != nCols) { |
174 | 2 | throw new IllegalArgumentException( |
175 | "All input rows must have the same length."); |
|
176 | } |
|
177 | } |
|
178 | 20 | this.copyIn(d); |
179 | 18 | lu = null; |
180 | 18 | } |
181 | ||
182 | /** |
|
183 | * Create a new (column) BigMatrix using <code>v</code> as the |
|
184 | * data for the unique column of the <code>v.length x 1</code> matrix |
|
185 | * created. |
|
186 | * <p> |
|
187 | * The input array is copied, not referenced. |
|
188 | * |
|
189 | * @param v column vector holding data for new matrix |
|
190 | */ |
|
191 | 2 | public BigMatrixImpl(BigDecimal[] v) { |
192 | 2 | int nRows = v.length; |
193 | 2 | data = new BigDecimal[nRows][1]; |
194 | 8 | for (int row = 0; row < nRows; row++) { |
195 | 6 | data[row][0] = v[row]; |
196 | } |
|
197 | 2 | } |
198 | ||
199 | /** |
|
200 | * Create a new BigMatrix which is a copy of this. |
|
201 | * |
|
202 | * @return the cloned matrix |
|
203 | */ |
|
204 | public BigMatrix copy() { |
|
205 | 2 | return new BigMatrixImpl(this.copyOut()); |
206 | } |
|
207 | ||
208 | /** |
|
209 | * Compute the sum of this and <code>m</code>. |
|
210 | * |
|
211 | * @param m matrix to be added |
|
212 | * @return this + m |
|
213 | * @exception IllegalArgumentException if m is not the same size as this |
|
214 | */ |
|
215 | public BigMatrix add(BigMatrix m) throws IllegalArgumentException { |
|
216 | 6 | if (this.getColumnDimension() != m.getColumnDimension() || |
217 | this.getRowDimension() != m.getRowDimension()) { |
|
218 | 2 | throw new IllegalArgumentException("matrix dimension mismatch"); |
219 | } |
|
220 | 4 | int rowCount = this.getRowDimension(); |
221 | 4 | int columnCount = this.getColumnDimension(); |
222 | 4 | BigDecimal[][] outData = new BigDecimal[rowCount][columnCount]; |
223 | 16 | for (int row = 0; row < rowCount; row++) { |
224 | 48 | for (int col = 0; col < columnCount; col++) { |
225 | 36 | outData[row][col] = data[row][col].add(m.getEntry(row, col)); |
226 | } |
|
227 | } |
|
228 | 4 | return new BigMatrixImpl(outData); |
229 | } |
|
230 | ||
231 | /** |
|
232 | * Compute this minus <code>m</code>. |
|
233 | * |
|
234 | * @param m matrix to be subtracted |
|
235 | * @return this + m |
|
236 | * @exception IllegalArgumentException if m is not the same size as *this |
|
237 | */ |
|
238 | public BigMatrix subtract(BigMatrix m) throws IllegalArgumentException { |
|
239 | 56 | if (this.getColumnDimension() != m.getColumnDimension() || |
240 | this.getRowDimension() != m.getRowDimension()) { |
|
241 | 2 | throw new IllegalArgumentException("matrix dimension mismatch"); |
242 | } |
|
243 | 54 | int rowCount = this.getRowDimension(); |
244 | 54 | int columnCount = this.getColumnDimension(); |
245 | 54 | BigDecimal[][] outData = new BigDecimal[rowCount][columnCount]; |
246 | 210 | for (int row = 0; row < rowCount; row++) { |
247 | 602 | for (int col = 0; col < columnCount; col++) { |
248 | 446 | outData[row][col] = data[row][col].subtract(m.getEntry(row, col)); |
249 | } |
|
250 | } |
|
251 | 54 | return new BigMatrixImpl(outData); |
252 | } |
|
253 | ||
254 | /** |
|
255 | * Returns the result of adding d to each entry of this. |
|
256 | * |
|
257 | * @param d value to be added to each entry |
|
258 | * @return d + this |
|
259 | */ |
|
260 | public BigMatrix scalarAdd(BigDecimal d) { |
|
261 | 2 | int rowCount = this.getRowDimension(); |
262 | 2 | int columnCount = this.getColumnDimension(); |
263 | 2 | BigDecimal[][] outData = new BigDecimal[rowCount][columnCount]; |
264 | 8 | for (int row = 0; row < rowCount; row++) { |
265 | 24 | for (int col = 0; col < columnCount; col++) { |
266 | 18 | outData[row][col] = data[row][col].add(d); |
267 | } |
|
268 | } |
|
269 | 2 | return new BigMatrixImpl(outData); |
270 | } |
|
271 | ||
272 | /** |
|
273 | * Returns the result multiplying each entry of this by <code>d</code> |
|
274 | * @param d value to multiply all entries by |
|
275 | * @return d * this |
|
276 | */ |
|
277 | public BigMatrix scalarMultiply(BigDecimal d) { |
|
278 | 2 | int rowCount = this.getRowDimension(); |
279 | 2 | int columnCount = this.getColumnDimension(); |
280 | 2 | BigDecimal[][] outData = new BigDecimal[rowCount][columnCount]; |
281 | 8 | for (int row = 0; row < rowCount; row++) { |
282 | 24 | for (int col = 0; col < columnCount; col++) { |
283 | 18 | outData[row][col] = data[row][col].multiply(d); |
284 | } |
|
285 | } |
|
286 | 2 | return new BigMatrixImpl(outData); |
287 | } |
|
288 | ||
289 | /** |
|
290 | * Returns the result of postmultiplying this by <code>m</code>. |
|
291 | * @param m matrix to postmultiply by |
|
292 | * @return this*m |
|
293 | * @throws IllegalArgumentException |
|
294 | * if columnDimension(this) != rowDimension(m) |
|
295 | */ |
|
296 | public BigMatrix multiply(BigMatrix m) throws IllegalArgumentException { |
|
297 | 34 | if (this.getColumnDimension() != m.getRowDimension()) { |
298 | 4 | throw new IllegalArgumentException("Matrices are not multiplication compatible."); |
299 | } |
|
300 | 30 | int nRows = this.getRowDimension(); |
301 | 30 | int nCols = m.getColumnDimension(); |
302 | 30 | int nSum = this.getColumnDimension(); |
303 | 30 | BigDecimal[][] outData = new BigDecimal[nRows][nCols]; |
304 | 30 | BigDecimal sum = ZERO; |
305 | 114 | for (int row = 0; row < nRows; row++) { |
306 | 320 | for (int col = 0; col < nCols; col++) { |
307 | 236 | sum = ZERO; |
308 | 952 | for (int i = 0; i < nSum; i++) { |
309 | 716 | sum = sum.add(data[row][i].multiply(m.getEntry(i, col))); |
310 | } |
|
311 | 236 | outData[row][col] = sum; |
312 | } |
|
313 | } |
|
314 | 30 | return new BigMatrixImpl(outData); |
315 | } |
|
316 | ||
317 | /** |
|
318 | * Returns the result premultiplying this by <code>m</code>. |
|
319 | * @param m matrix to premultiply by |
|
320 | * @return m * this |
|
321 | * @throws IllegalArgumentException |
|
322 | * if rowDimension(this) != columnDimension(m) |
|
323 | */ |
|
324 | public BigMatrix preMultiply(BigMatrix m) throws IllegalArgumentException { |
|
325 | 12 | return m.multiply(this); |
326 | } |
|
327 | ||
328 | /** |
|
329 | * Returns matrix entries as a two-dimensional array. |
|
330 | * <p> |
|
331 | * Makes a fresh copy of the underlying data. |
|
332 | * |
|
333 | * @return 2-dimensional array of entries |
|
334 | */ |
|
335 | public BigDecimal[][] getData() { |
|
336 | 50 | return copyOut(); |
337 | } |
|
338 | ||
339 | /** |
|
340 | * Returns matrix entries as a two-dimensional array. |
|
341 | * <p> |
|
342 | * Makes a fresh copy of the underlying data converted to |
|
343 | * <code>double</code> values. |
|
344 | * |
|
345 | * @return 2-dimensional array of entries |
|
346 | */ |
|
347 | public double[][] getDataAsDoubleArray() { |
|
348 | 0 | int nRows = getRowDimension(); |
349 | 0 | int nCols = getColumnDimension(); |
350 | 0 | double d[][] = new double[nRows][nCols]; |
351 | 0 | for (int i = 0; i < nRows; i++) { |
352 | 0 | for (int j=0; j<nCols;j++) { |
353 | 0 | d[i][j] = data[i][j].doubleValue(); |
354 | } |
|
355 | } |
|
356 | 0 | return d; |
357 | } |
|
358 | ||
359 | /** |
|
360 | * Returns a reference to the underlying data array. |
|
361 | * <p> |
|
362 | * Does not make a fresh copy of the underlying data. |
|
363 | * |
|
364 | * @return 2-dimensional array of entries |
|
365 | */ |
|
366 | public BigDecimal[][] getDataRef() { |
|
367 | 50 | return data; |
368 | } |
|
369 | ||
370 | /*** |
|
371 | * Gets the rounding mode for division operations |
|
372 | * The default is {@link java.math.BigDecimal#ROUND_HALF_UP} |
|
373 | * @see BigDecimal |
|
374 | * @return the rounding mode. |
|
375 | */ |
|
376 | public int getRoundingMode() { |
|
377 | 0 | return roundingMode; |
378 | } |
|
379 | ||
380 | /*** |
|
381 | * Sets the rounding mode for decimal divisions. |
|
382 | * @see BigDecimal |
|
383 | * @param roundingMode |
|
384 | */ |
|
385 | public void setRoundingMode(int roundingMode) { |
|
386 | 0 | this.roundingMode = roundingMode; |
387 | 0 | } |
388 | ||
389 | /*** |
|
390 | * Sets the scale for division operations. |
|
391 | * The default is 64 |
|
392 | * @see BigDecimal |
|
393 | * @return the scale |
|
394 | */ |
|
395 | public int getScale() { |
|
396 | 0 | return scale; |
397 | } |
|
398 | ||
399 | /*** |
|
400 | * Sets the scale for division operations. |
|
401 | * @see BigDecimal |
|
402 | * @param scale |
|
403 | */ |
|
404 | public void setScale(int scale) { |
|
405 | 0 | this.scale = scale; |
406 | 0 | } |
407 | ||
408 | /** |
|
409 | * Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html"> |
|
410 | * maximum absolute row sum norm</a> of the matrix. |
|
411 | * |
|
412 | * @return norm |
|
413 | */ |
|
414 | public BigDecimal getNorm() { |
|
415 | 56 | BigDecimal maxColSum = ZERO; |
416 | 214 | for (int col = 0; col < this.getColumnDimension(); col++) { |
417 | 158 | BigDecimal sum = ZERO; |
418 | 616 | for (int row = 0; row < this.getRowDimension(); row++) { |
419 | 458 | sum = sum.add(data[row][col].abs()); |
420 | } |
|
421 | 158 | maxColSum = maxColSum.max(sum); |
422 | } |
|
423 | 56 | return maxColSum; |
424 | } |
|
425 | ||
426 | /** |
|
427 | * Gets a submatrix. Rows and columns are indicated |
|
428 | * counting from 0 to n-1. |
|
429 | * |
|
430 | * @param startRow Initial row index |
|
431 | * @param endRow Final row index |
|
432 | * @param startColumn Initial column index |
|
433 | * @param endColumn Final column index |
|
434 | * @return The subMatrix containing the data of the |
|
435 | * specified rows and columns |
|
436 | * @exception MatrixIndexException if row or column selections are not valid |
|
437 | */ |
|
438 | public BigMatrix getSubMatrix(int startRow, int endRow, int startColumn, |
|
439 | int endColumn) throws MatrixIndexException { |
|
440 | 14 | if (startRow < 0 || startRow > endRow || endRow > data.length || |
441 | startColumn < 0 || startColumn > endColumn || |
|
442 | endColumn > data[0].length ) { |
|
443 | 8 | throw new MatrixIndexException( |
444 | "invalid row or column index selection"); |
|
445 | } |
|
446 | 6 | BigMatrixImpl subMatrix = new BigMatrixImpl(endRow - startRow+1, |
447 | endColumn - startColumn+1); |
|
448 | 6 | BigDecimal[][] subMatrixData = subMatrix.getDataRef(); |
449 | 16 | for (int i = startRow; i <= endRow; i++) { |
450 | 24 | for (int j = startColumn; j <= endColumn; j++) { |
451 | 14 | subMatrixData[i - startRow][j - startColumn] = data[i][j]; |
452 | } |
|
453 | } |
|
454 | 6 | return subMatrix; |
455 | } |
|
456 | ||
457 | /** |
|
458 | * Gets a submatrix. Rows and columns are indicated |
|
459 | * counting from 0 to n-1. |
|
460 | * |
|
461 | * @param selectedRows Array of row indices must be non-empty |
|
462 | * @param selectedColumns Array of column indices must be non-empty |
|
463 | * @return The subMatrix containing the data in the |
|
464 | * specified rows and columns |
|
465 | * @exception MatrixIndexException if supplied row or column index arrays |
|
466 | * are not valid |
|
467 | */ |
|
468 | public BigMatrix getSubMatrix(int[] selectedRows, int[] selectedColumns) |
|
469 | throws MatrixIndexException { |
|
470 | 16 | if (selectedRows.length * selectedColumns.length == 0) { |
471 | 2 | throw new MatrixIndexException( |
472 | "selected row and column index arrays must be non-empty"); |
|
473 | } |
|
474 | 14 | BigMatrixImpl subMatrix = new BigMatrixImpl(selectedRows.length, |
475 | selectedColumns.length); |
|
476 | 14 | BigDecimal[][] subMatrixData = subMatrix.getDataRef(); |
477 | try { |
|
478 | 38 | for (int i = 0; i < selectedRows.length; i++) { |
479 | 82 | for (int j = 0; j < selectedColumns.length; j++) { |
480 | 58 | subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]]; |
481 | } |
|
482 | } |
|
483 | } |
|
484 | 2 | catch (ArrayIndexOutOfBoundsException e) { |
485 | 2 | throw new MatrixIndexException("matrix dimension mismatch"); |
486 | 12 | } |
487 | 12 | return subMatrix; |
488 | } |
|
489 | ||
490 | /** |
|
491 | * Replace the submatrix starting at <code>row, column</code> using data in |
|
492 | * the input <code>subMatrix</code> array. Indexes are 0-based. |
|
493 | * <p> |
|
494 | * Example:<br> |
|
495 | * Starting with <pre> |
|
496 | * 1 2 3 4 |
|
497 | * 5 6 7 8 |
|
498 | * 9 0 1 2 |
|
499 | * </pre> |
|
500 | * and <code>subMatrix = {{3, 4} {5,6}}</code>, invoking |
|
501 | * <code>setSubMatrix(subMatrix,1,1))</code> will result in <pre> |
|
502 | * 1 2 3 4 |
|
503 | * 5 3 4 8 |
|
504 | * 9 5 6 2 |
|
505 | * </pre> |
|
506 | * |
|
507 | * @param subMatrix array containing the submatrix replacement data |
|
508 | * @param row row coordinate of the top, left element to be replaced |
|
509 | * @param column column coordinate of the top, left element to be replaced |
|
510 | * @throws MatrixIndexException if subMatrix does not fit into this |
|
511 | * matrix from element in (row, column) |
|
512 | * @throws IllegalArgumentException if <code>subMatrix</code> is not rectangular |
|
513 | * (not all rows have the same length) or empty |
|
514 | * @throws NullPointerException if <code>subMatrix</code> is null |
|
515 | * @since 1.1 |
|
516 | */ |
|
517 | public void setSubMatrix(BigDecimal[][] subMatrix, int row, int column) |
|
518 | throws MatrixIndexException { |
|
519 | 180 | if ((row < 0) || (column < 0)){ |
520 | 0 | throw new MatrixIndexException |
521 | ("invalid row or column index selection"); |
|
522 | } |
|
523 | 180 | int nRows = subMatrix.length; |
524 | 178 | if (nRows == 0) { |
525 | 0 | throw new IllegalArgumentException( |
526 | "Matrix must have at least one row."); |
|
527 | } |
|
528 | 178 | int nCols = subMatrix[0].length; |
529 | 178 | if (nCols == 0) { |
530 | 2 | throw new IllegalArgumentException( |
531 | "Matrix must have at least one column."); |
|
532 | } |
|
533 | 494 | for (int r = 1; r < nRows; r++) { |
534 | 320 | if (subMatrix[r].length != nCols) { |
535 | 2 | throw new IllegalArgumentException( |
536 | "All input rows must have the same length."); |
|
537 | } |
|
538 | } |
|
539 | 174 | if (data == null) { |
540 | 164 | if ((row > 0)||(column > 0)) throw new MatrixIndexException |
541 | ("matrix must be initialized to perfom this method"); |
|
542 | 164 | data = new BigDecimal[nRows][nCols]; |
543 | 164 | System.arraycopy(subMatrix, 0, data, 0, subMatrix.length); |
544 | } |
|
545 | 174 | if (((nRows + row) > this.getRowDimension()) || |
546 | (nCols + column > this.getColumnDimension())) |
|
547 | 2 | throw new MatrixIndexException( |
548 | "invalid row or column index selection"); |
|
549 | 658 | for (int i = 0; i < nRows; i++) { |
550 | 486 | System.arraycopy(subMatrix[i], 0, data[row + i], column, nCols); |
551 | } |
|
552 | 172 | lu = null; |
553 | 172 | } |
554 | ||
555 | /** |
|
556 | * Returns the entries in row number <code>row</code> |
|
557 | * as a row matrix. Row indices start at 0. |
|
558 | * |
|
559 | * @param row the row to be fetched |
|
560 | * @return row matrix |
|
561 | * @throws MatrixIndexException if the specified row index is invalid |
|
562 | */ |
|
563 | public BigMatrix getRowMatrix(int row) throws MatrixIndexException { |
|
564 | 8 | if ( !isValidCoordinate( row, 0)) { |
565 | 4 | throw new MatrixIndexException("illegal row argument"); |
566 | } |
|
567 | 4 | int ncols = this.getColumnDimension(); |
568 | 4 | BigDecimal[][] out = new BigDecimal[1][ncols]; |
569 | 4 | System.arraycopy(data[row], 0, out[0], 0, ncols); |
570 | 4 | return new BigMatrixImpl(out); |
571 | } |
|
572 | ||
573 | /** |
|
574 | * Returns the entries in column number <code>column</code> |
|
575 | * as a column matrix. Column indices start at 0. |
|
576 | * |
|
577 | * @param column the column to be fetched |
|
578 | * @return column matrix |
|
579 | * @throws MatrixIndexException if the specified column index is invalid |
|
580 | */ |
|
581 | public BigMatrix getColumnMatrix(int column) throws MatrixIndexException { |
|
582 | 8 | if ( !isValidCoordinate( 0, column)) { |
583 | 4 | throw new MatrixIndexException("illegal column argument"); |
584 | } |
|
585 | 4 | int nRows = this.getRowDimension(); |
586 | 4 | BigDecimal[][] out = new BigDecimal[nRows][1]; |
587 | 20 | for (int row = 0; row < nRows; row++) { |
588 | 16 | out[row][0] = data[row][column]; |
589 | } |
|
590 | 4 | return new BigMatrixImpl(out); |
591 | } |
|
592 | ||
593 | /** |
|
594 | * Returns the entries in row number <code>row</code> as an array. |
|
595 | * <p> |
|
596 | * Row indices start at 0. A <code>MatrixIndexException</code> is thrown |
|
597 | * unless <code>0 <= row < rowDimension.</code> |
|
598 | * |
|
599 | * @param row the row to be fetched |
|
600 | * @return array of entries in the row |
|
601 | * @throws MatrixIndexException if the specified row index is not valid |
|
602 | */ |
|
603 | public BigDecimal[] getRow(int row) throws MatrixIndexException { |
|
604 | 0 | if ( !isValidCoordinate( row, 0 ) ) { |
605 | 0 | throw new MatrixIndexException("illegal row argument"); |
606 | } |
|
607 | 0 | int ncols = this.getColumnDimension(); |
608 | 0 | BigDecimal[] out = new BigDecimal[ncols]; |
609 | 0 | System.arraycopy(data[row], 0, out, 0, ncols); |
610 | 0 | return out; |
611 | } |
|
612 | ||
613 | /** |
|
614 | * Returns the entries in row number <code>row</code> as an array |
|
615 | * of double values. |
|
616 | * <p> |
|
617 | * Row indices start at 0. A <code>MatrixIndexException</code> is thrown |
|
618 | * unless <code>0 <= row < rowDimension.</code> |
|
619 | * |
|
620 | * @param row the row to be fetched |
|
621 | * @return array of entries in the row |
|
622 | * @throws MatrixIndexException if the specified row index is not valid |
|
623 | */ |
|
624 | public double[] getRowAsDoubleArray(int row) throws MatrixIndexException { |
|
625 | 4 | if ( !isValidCoordinate( row, 0 ) ) { |
626 | 2 | throw new MatrixIndexException("illegal row argument"); |
627 | } |
|
628 | 2 | int ncols = this.getColumnDimension(); |
629 | 2 | double[] out = new double[ncols]; |
630 | 8 | for (int i=0;i<ncols;i++) { |
631 | 6 | out[i] = data[row][i].doubleValue(); |
632 | } |
|
633 | 2 | return out; |
634 | } |
|
635 | ||
636 | /** |
|
637 | * Returns the entries in column number <code>col</code> as an array. |
|
638 | * <p> |
|
639 | * Column indices start at 0. A <code>MatrixIndexException</code> is thrown |
|
640 | * unless <code>0 <= column < columnDimension.</code> |
|
641 | * |
|
642 | * @param col the column to be fetched |
|
643 | * @return array of entries in the column |
|
644 | * @throws MatrixIndexException if the specified column index is not valid |
|
645 | */ |
|
646 | public BigDecimal[] getColumn(int col) throws MatrixIndexException { |
|
647 | 0 | if ( !isValidCoordinate(0, col) ) { |
648 | 0 | throw new MatrixIndexException("illegal column argument"); |
649 | } |
|
650 | 0 | int nRows = this.getRowDimension(); |
651 | 0 | BigDecimal[] out = new BigDecimal[nRows]; |
652 | 0 | for (int i = 0; i < nRows; i++) { |
653 | 0 | out[i] = data[i][col]; |
654 | } |
|
655 | 0 | return out; |
656 | } |
|
657 | ||
658 | /** |
|
659 | * Returns the entries in column number <code>col</code> as an array |
|
660 | * of double values. |
|
661 | * <p> |
|
662 | * Column indices start at 0. A <code>MatrixIndexException</code> is thrown |
|
663 | * unless <code>0 <= column < columnDimension.</code> |
|
664 | * |
|
665 | * @param col the column to be fetched |
|
666 | * @return array of entries in the column |
|
667 | * @throws MatrixIndexException if the specified column index is not valid |
|
668 | */ |
|
669 | public double[] getColumnAsDoubleArray(int col) throws MatrixIndexException { |
|
670 | 4 | if ( !isValidCoordinate( 0, col ) ) { |
671 | 2 | throw new MatrixIndexException("illegal column argument"); |
672 | } |
|
673 | 2 | int nrows = this.getRowDimension(); |
674 | 2 | double[] out = new double[nrows]; |
675 | 8 | for (int i=0;i<nrows;i++) { |
676 | 6 | out[i] = data[i][col].doubleValue(); |
677 | } |
|
678 | 2 | return out; |
679 | } |
|
680 | ||
681 | /** |
|
682 | * Returns the entry in the specified row and column. |
|
683 | * <p> |
|
684 | * Row and column indices start at 0 and must satisfy |
|
685 | * <ul> |
|
686 | * <li><code>0 <= row < rowDimension</code></li> |
|
687 | * <li><code> 0 <= column < columnDimension</code></li> |
|
688 | * </ul> |
|
689 | * otherwise a <code>MatrixIndexException</code> is thrown. |
|
690 | * |
|
691 | * @param row row location of entry to be fetched |
|
692 | * @param column column location of entry to be fetched |
|
693 | * @return matrix entry in row,column |
|
694 | * @throws MatrixIndexException if the row or column index is not valid |
|
695 | */ |
|
696 | public BigDecimal getEntry(int row, int column) |
|
697 | throws MatrixIndexException { |
|
698 | 1756 | if (!isValidCoordinate(row,column)) { |
699 | 0 | throw new MatrixIndexException("matrix entry does not exist"); |
700 | } |
|
701 | 1756 | return data[row][column]; |
702 | } |
|
703 | ||
704 | /** |
|
705 | * Returns the entry in the specified row and column as a double. |
|
706 | * <p> |
|
707 | * Row and column indices start at 0 and must satisfy |
|
708 | * <ul> |
|
709 | * <li><code>0 <= row < rowDimension</code></li> |
|
710 | * <li><code> 0 <= column < columnDimension</code></li> |
|
711 | * </ul> |
|
712 | * otherwise a <code>MatrixIndexException</code> is thrown. |
|
713 | * |
|
714 | * @param row row location of entry to be fetched |
|
715 | * @param column column location of entry to be fetched |
|
716 | * @return matrix entry in row,column |
|
717 | * @throws MatrixIndexException if the row |
|
718 | * or column index is not valid |
|
719 | */ |
|
720 | public double getEntryAsDouble(int row, int column) throws MatrixIndexException { |
|
721 | 0 | return getEntry(row,column).doubleValue(); |
722 | } |
|
723 | ||
724 | /** |
|
725 | * Returns the transpose matrix. |
|
726 | * |
|
727 | * @return transpose matrix |
|
728 | */ |
|
729 | public BigMatrix transpose() { |
|
730 | 8 | int nRows = this.getRowDimension(); |
731 | 8 | int nCols = this.getColumnDimension(); |
732 | 8 | BigMatrixImpl out = new BigMatrixImpl(nCols, nRows); |
733 | 8 | BigDecimal[][] outData = out.getDataRef(); |
734 | 30 | for (int row = 0; row < nRows; row++) { |
735 | 88 | for (int col = 0; col < nCols; col++) { |
736 | 66 | outData[col][row] = data[row][col]; |
737 | } |
|
738 | } |
|
739 | 8 | return out; |
740 | } |
|
741 | ||
742 | /** |
|
743 | * Returns the inverse matrix if this matrix is invertible. |
|
744 | * |
|
745 | * @return inverse matrix |
|
746 | * @throws InvalidMatrixException if this is not invertible |
|
747 | */ |
|
748 | public BigMatrix inverse() throws InvalidMatrixException { |
|
749 | 14 | return solve(MatrixUtils.createBigIdentityMatrix |
750 | (this.getRowDimension())); |
|
751 | } |
|
752 | ||
753 | /** |
|
754 | * Returns the determinant of this matrix. |
|
755 | * |
|
756 | * @return determinant |
|
757 | * @throws InvalidMatrixException if matrix is not square |
|
758 | */ |
|
759 | public BigDecimal getDeterminant() throws InvalidMatrixException { |
|
760 | 10 | if (!isSquare()) { |
761 | 2 | throw new InvalidMatrixException("matrix is not square"); |
762 | } |
|
763 | 8 | if (isSingular()) { // note: this has side effect of attempting LU decomp if lu == null |
764 | 2 | return ZERO; |
765 | } else { |
|
766 | 6 | BigDecimal det = (parity == 1) ? ONE : ONE.negate(); |
767 | 22 | for (int i = 0; i < this.getRowDimension(); i++) { |
768 | 16 | det = det.multiply(lu[i][i]); |
769 | } |
|
770 | 6 | return det; |
771 | } |
|
772 | } |
|
773 | ||
774 | /** |
|
775 | * Is this a square matrix? |
|
776 | * @return true if the matrix is square (rowDimension = columnDimension) |
|
777 | */ |
|
778 | public boolean isSquare() { |
|
779 | 52 | return (this.getColumnDimension() == this.getRowDimension()); |
780 | } |
|
781 | ||
782 | /** |
|
783 | * Is this a singular matrix? |
|
784 | * @return true if the matrix is singular |
|
785 | */ |
|
786 | public boolean isSingular() { |
|
787 | 32 | if (lu == null) { |
788 | try { |
|
789 | 30 | luDecompose(); |
790 | 20 | return false; |
791 | 10 | } catch (InvalidMatrixException ex) { |
792 | 10 | return true; |
793 | } |
|
794 | } else { // LU decomp must have been successfully performed |
|
795 | 2 | return false; // so the matrix is not singular |
796 | } |
|
797 | } |
|
798 | ||
799 | /** |
|
800 | * Returns the number of rows in the matrix. |
|
801 | * |
|
802 | * @return rowDimension |
|
803 | */ |
|
804 | public int getRowDimension() { |
|
805 | 3296 | return data.length; |
806 | } |
|
807 | ||
808 | /** |
|
809 | * Returns the number of columns in the matrix. |
|
810 | * |
|
811 | * @return columnDimension |
|
812 | */ |
|
813 | public int getColumnDimension() { |
|
814 | 2856 | return data[0].length; |
815 | } |
|
816 | ||
817 | /** |
|
818 | * Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html"> |
|
819 | * trace</a> of the matrix (the sum of the elements on the main diagonal). |
|
820 | * |
|
821 | * @return trace |
|
822 | * |
|
823 | * @throws IllegalArgumentException if this matrix is not square. |
|
824 | */ |
|
825 | public BigDecimal getTrace() throws IllegalArgumentException { |
|
826 | 4 | if (!isSquare()) { |
827 | 2 | throw new IllegalArgumentException("matrix is not square"); |
828 | } |
|
829 | 2 | BigDecimal trace = data[0][0]; |
830 | 6 | for (int i = 1; i < this.getRowDimension(); i++) { |
831 | 4 | trace = trace.add(data[i][i]); |
832 | } |
|
833 | 2 | return trace; |
834 | } |
|
835 | ||
836 | /** |
|
837 | * Returns the result of multiplying this by the vector <code>v</code>. |
|
838 | * |
|
839 | * @param v the vector to operate on |
|
840 | * @return this*v |
|
841 | * @throws IllegalArgumentException if columnDimension != v.size() |
|
842 | */ |
|
843 | public BigDecimal[] operate(BigDecimal[] v) throws IllegalArgumentException { |
|
844 | 6 | if (v.length != this.getColumnDimension()) { |
845 | 2 | throw new IllegalArgumentException("vector has wrong length"); |
846 | } |
|
847 | 4 | int nRows = this.getRowDimension(); |
848 | 4 | int nCols = this.getColumnDimension(); |
849 | 4 | BigDecimal[] out = new BigDecimal[v.length]; |
850 | 16 | for (int row = 0; row < nRows; row++) { |
851 | 12 | BigDecimal sum = ZERO; |
852 | 48 | for (int i = 0; i < nCols; i++) { |
853 | 36 | sum = sum.add(data[row][i].multiply(v[i])); |
854 | } |
|
855 | 12 | out[row] = sum; |
856 | } |
|
857 | 4 | return out; |
858 | } |
|
859 | ||
860 | /** |
|
861 | * Returns the result of multiplying this by the vector <code>v</code>. |
|
862 | * |
|
863 | * @param v the vector to operate on |
|
864 | * @return this*v |
|
865 | * @throws IllegalArgumentException if columnDimension != v.size() |
|
866 | */ |
|
867 | public BigDecimal[] operate(double[] v) throws IllegalArgumentException { |
|
868 | 0 | BigDecimal bd[] = new BigDecimal[v.length]; |
869 | 0 | for (int i=0;i<bd.length;i++) { |
870 | 0 | bd[i] = new BigDecimal(v[i]); |
871 | } |
|
872 | 0 | return operate(bd); |
873 | } |
|
874 | ||
875 | /** |
|
876 | * Returns the (row) vector result of premultiplying this by the vector <code>v</code>. |
|
877 | * |
|
878 | * @param v the row vector to premultiply by |
|
879 | * @return v*this |
|
880 | * @throws IllegalArgumentException if rowDimension != v.size() |
|
881 | */ |
|
882 | public BigDecimal[] preMultiply(BigDecimal[] v) throws IllegalArgumentException { |
|
883 | 4 | int nRows = this.getRowDimension(); |
884 | 4 | if (v.length != nRows) { |
885 | 2 | throw new IllegalArgumentException("vector has wrong length"); |
886 | } |
|
887 | 2 | int nCols = this.getColumnDimension(); |
888 | 2 | BigDecimal[] out = new BigDecimal[nCols]; |
889 | 8 | for (int col = 0; col < nCols; col++) { |
890 | 6 | BigDecimal sum = ZERO; |
891 | 24 | for (int i = 0; i < nRows; i++) { |
892 | 18 | sum = sum.add(data[i][col].multiply(v[i])); |
893 | } |
|
894 | 6 | out[col] = sum; |
895 | } |
|
896 | 2 | return out; |
897 | } |
|
898 | ||
899 | /** |
|
900 | * Returns a matrix of (column) solution vectors for linear systems with |
|
901 | * coefficient matrix = this and constant vectors = columns of |
|
902 | * <code>b</code>. |
|
903 | * |
|
904 | * @param b array of constants forming RHS of linear systems to |
|
905 | * to solve |
|
906 | * @return solution array |
|
907 | * @throws IllegalArgumentException if this.rowDimension != row dimension |
|
908 | * @throws InvalidMatrixException if this matrix is not square or is singular |
|
909 | */ |
|
910 | public BigDecimal[] solve(BigDecimal[] b) throws IllegalArgumentException, InvalidMatrixException { |
|
911 | 4 | int nRows = this.getRowDimension(); |
912 | 4 | if (b.length != nRows) { |
913 | 2 | throw new IllegalArgumentException("constant vector has wrong length"); |
914 | } |
|
915 | 2 | BigMatrix bMatrix = new BigMatrixImpl(b); |
916 | 2 | BigDecimal[][] solution = ((BigMatrixImpl) (solve(bMatrix))).getDataRef(); |
917 | 2 | BigDecimal[] out = new BigDecimal[nRows]; |
918 | 8 | for (int row = 0; row < nRows; row++) { |
919 | 6 | out[row] = solution[row][0]; |
920 | } |
|
921 | 2 | return out; |
922 | } |
|
923 | ||
924 | /** |
|
925 | * Returns a matrix of (column) solution vectors for linear systems with |
|
926 | * coefficient matrix = this and constant vectors = columns of |
|
927 | * <code>b</code>. |
|
928 | * |
|
929 | * @param b array of constants forming RHS of linear systems to |
|
930 | * to solve |
|
931 | * @return solution array |
|
932 | * @throws IllegalArgumentException if this.rowDimension != row dimension |
|
933 | * @throws InvalidMatrixException if this matrix is not square or is singular |
|
934 | */ |
|
935 | public BigDecimal[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException { |
|
936 | 0 | BigDecimal bd[] = new BigDecimal[b.length]; |
937 | 0 | for (int i=0;i<bd.length;i++) { |
938 | 0 | bd[i] = new BigDecimal(b[i]); |
939 | } |
|
940 | 0 | return solve(bd); |
941 | } |
|
942 | ||
943 | /** |
|
944 | * Returns a matrix of (column) solution vectors for linear systems with |
|
945 | * coefficient matrix = this and constant vectors = columns of |
|
946 | * <code>b</code>. |
|
947 | * |
|
948 | * @param b matrix of constant vectors forming RHS of linear systems to |
|
949 | * to solve |
|
950 | * @return matrix of solution vectors |
|
951 | * @throws IllegalArgumentException if this.rowDimension != row dimension |
|
952 | * @throws InvalidMatrixException if this matrix is not square or is singular |
|
953 | */ |
|
954 | public BigMatrix solve(BigMatrix b) throws IllegalArgumentException, InvalidMatrixException { |
|
955 | 22 | if (b.getRowDimension() != this.getRowDimension()) { |
956 | 4 | throw new IllegalArgumentException("Incorrect row dimension"); |
957 | } |
|
958 | 18 | if (!this.isSquare()) { |
959 | 2 | throw new InvalidMatrixException("coefficient matrix is not square"); |
960 | } |
|
961 | 16 | if (this.isSingular()) { // side effect: compute LU decomp |
962 | 4 | throw new InvalidMatrixException("Matrix is singular."); |
963 | } |
|
964 | ||
965 | 12 | int nCol = this.getColumnDimension(); |
966 | 12 | int nColB = b.getColumnDimension(); |
967 | 12 | int nRowB = b.getRowDimension(); |
968 | ||
969 | // Apply permutations to b |
|
970 | 12 | BigDecimal[][] bp = new BigDecimal[nRowB][nColB]; |
971 | 48 | for (int row = 0; row < nRowB; row++) { |
972 | 132 | for (int col = 0; col < nColB; col++) { |
973 | 96 | bp[row][col] = b.getEntry(permutation[row], col); |
974 | } |
|
975 | } |
|
976 | ||
977 | // Solve LY = b |
|
978 | 48 | for (int col = 0; col < nCol; col++) { |
979 | 72 | for (int i = col + 1; i < nCol; i++) { |
980 | 132 | for (int j = 0; j < nColB; j++) { |
981 | 96 | bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col])); |
982 | } |
|
983 | } |
|
984 | } |
|
985 | ||
986 | // Solve UX = Y |
|
987 | 48 | for (int col = nCol - 1; col >= 0; col--) { |
988 | 132 | for (int j = 0; j < nColB; j++) { |
989 | 96 | bp[col][j] = bp[col][j].divide(lu[col][col], scale, roundingMode); |
990 | } |
|
991 | 72 | for (int i = 0; i < col; i++) { |
992 | 132 | for (int j = 0; j < nColB; j++) { |
993 | 96 | bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col])); |
994 | } |
|
995 | } |
|
996 | } |
|
997 | ||
998 | 12 | BigMatrixImpl outMat = new BigMatrixImpl(bp); |
999 | 12 | return outMat; |
1000 | } |
|
1001 | ||
1002 | /** |
|
1003 | * Computes a new |
|
1004 | * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf"> |
|
1005 | * LU decompostion</a> for this matrix, storing the result for use by other methods. |
|
1006 | * <p> |
|
1007 | * <strong>Implementation Note</strong>:<br> |
|
1008 | * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm"> |
|
1009 | * Crout's algortithm</a>, with partial pivoting. |
|
1010 | * <p> |
|
1011 | * <strong>Usage Note</strong>:<br> |
|
1012 | * This method should rarely be invoked directly. Its only use is |
|
1013 | * to force recomputation of the LU decomposition when changes have been |
|
1014 | * made to the underlying data using direct array references. Changes |
|
1015 | * made using setXxx methods will trigger recomputation when needed |
|
1016 | * automatically. |
|
1017 | * |
|
1018 | * @throws InvalidMatrixException if the matrix is non-square or singular. |
|
1019 | */ |
|
1020 | public void luDecompose() throws InvalidMatrixException { |
|
1021 | ||
1022 | 44 | int nRows = this.getRowDimension(); |
1023 | 44 | int nCols = this.getColumnDimension(); |
1024 | 44 | if (nRows != nCols) { |
1025 | 4 | throw new InvalidMatrixException("LU decomposition requires that the matrix be square."); |
1026 | } |
|
1027 | 40 | lu = this.getData(); |
1028 | ||
1029 | // Initialize permutation array and parity |
|
1030 | 40 | permutation = new int[nRows]; |
1031 | 162 | for (int row = 0; row < nRows; row++) { |
1032 | 122 | permutation[row] = row; |
1033 | } |
|
1034 | 40 | parity = 1; |
1035 | ||
1036 | // Loop over columns |
|
1037 | 150 | for (int col = 0; col < nCols; col++) { |
1038 | ||
1039 | 122 | BigDecimal sum = ZERO; |
1040 | ||
1041 | // upper |
|
1042 | 254 | for (int row = 0; row < col; row++) { |
1043 | 132 | sum = lu[row][col]; |
1044 | 190 | for (int i = 0; i < row; i++) { |
1045 | 58 | sum = sum.subtract(lu[row][i].multiply(lu[i][col])); |
1046 | } |
|
1047 | 132 | lu[row][col] = sum; |
1048 | } |
|
1049 | ||
1050 | // lower |
|
1051 | 122 | int max = col; // permutation row |
1052 | 122 | BigDecimal largest = ZERO; |
1053 | 376 | for (int row = col; row < nRows; row++) { |
1054 | 254 | sum = lu[row][col]; |
1055 | 444 | for (int i = 0; i < col; i++) { |
1056 | 190 | sum = sum.subtract(lu[row][i].multiply(lu[i][col])); |
1057 | } |
|
1058 | 254 | lu[row][col] = sum; |
1059 | ||
1060 | // maintain best permutation choice |
|
1061 | 254 | if (sum.abs().compareTo(largest) == 1) { |
1062 | 176 | largest = sum.abs(); |
1063 | 176 | max = row; |
1064 | } |
|
1065 | } |
|
1066 | ||
1067 | // Singularity check |
|
1068 | 122 | if (lu[max][col].abs().compareTo(TOO_SMALL) <= 0) { |
1069 | 12 | lu = null; |
1070 | 12 | throw new InvalidMatrixException("matrix is singular"); |
1071 | } |
|
1072 | ||
1073 | // Pivot if necessary |
|
1074 | 110 | if (max != col) { |
1075 | 54 | BigDecimal tmp = ZERO; |
1076 | 230 | for (int i = 0; i < nCols; i++) { |
1077 | 176 | tmp = lu[max][i]; |
1078 | 176 | lu[max][i] = lu[col][i]; |
1079 | 176 | lu[col][i] = tmp; |
1080 | } |
|
1081 | 54 | int temp = permutation[max]; |
1082 | 54 | permutation[max] = permutation[col]; |
1083 | 54 | permutation[col] = temp; |
1084 | 54 | parity = -parity; |
1085 | } |
|
1086 | ||
1087 | //Divide the lower elements by the "winning" diagonal elt. |
|
1088 | 242 | for (int row = col + 1; row < nRows; row++) { |
1089 | 132 | lu[row][col] = lu[row][col].divide(lu[col][col], scale, roundingMode); |
1090 | } |
|
1091 | ||
1092 | } |
|
1093 | ||
1094 | 28 | } |
1095 | ||
1096 | /** |
|
1097 | * |
|
1098 | * @see Object#toString() |
|
1099 | */ |
|
1100 | public String toString() { |
|
1101 | 4 | StringBuffer res = new StringBuffer(); |
1102 | 4 | res.append("BigMatrixImpl{"); |
1103 | 4 | if (data != null) { |
1104 | 8 | for (int i = 0; i < data.length; i++) { |
1105 | 6 | if (i > 0) |
1106 | 4 | res.append(","); |
1107 | 6 | res.append("{"); |
1108 | 24 | for (int j = 0; j < data[0].length; j++) { |
1109 | 18 | if (j > 0) |
1110 | 12 | res.append(","); |
1111 | 18 | res.append(data[i][j]); |
1112 | } |
|
1113 | 6 | res.append("}"); |
1114 | } |
|
1115 | } |
|
1116 | 4 | res.append("}"); |
1117 | 4 | return res.toString(); |
1118 | } |
|
1119 | ||
1120 | /** |
|
1121 | * Returns true iff <code>object</code> is a |
|
1122 | * <code>BigMatrixImpl</code> instance with the same dimensions as this |
|
1123 | * and all corresponding matrix entries are equal. BigDecimal.equals |
|
1124 | * is used to compare corresponding entries. |
|
1125 | * |
|
1126 | * @param object the object to test equality against. |
|
1127 | * @return true if object equals this |
|
1128 | */ |
|
1129 | public boolean equals(Object object) { |
|
1130 | 66 | if (object == this ) { |
1131 | 2 | return true; |
1132 | } |
|
1133 | 64 | if (object instanceof BigMatrixImpl == false) { |
1134 | 2 | return false; |
1135 | } |
|
1136 | 62 | BigMatrix m = (BigMatrix) object; |
1137 | 62 | int nRows = getRowDimension(); |
1138 | 62 | int nCols = getColumnDimension(); |
1139 | 62 | if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) { |
1140 | 2 | return false; |
1141 | } |
|
1142 | 192 | for (int row = 0; row < nRows; row++) { |
1143 | 422 | for (int col = 0; col < nCols; col++) { |
1144 | 290 | if (!data[row][col].equals(m.getEntry(row, col))) { |
1145 | 4 | return false; |
1146 | } |
|
1147 | } |
|
1148 | } |
|
1149 | 56 | return true; |
1150 | } |
|
1151 | ||
1152 | /** |
|
1153 | * Computes a hashcode for the matrix. |
|
1154 | * |
|
1155 | * @return hashcode for matrix |
|
1156 | */ |
|
1157 | public int hashCode() { |
|
1158 | 12 | int ret = 7; |
1159 | 12 | int nRows = getRowDimension(); |
1160 | 12 | int nCols = getColumnDimension(); |
1161 | 12 | ret = ret * 31 + nRows; |
1162 | 12 | ret = ret * 31 + nCols; |
1163 | 40 | for (int row = 0; row < nRows; row++) { |
1164 | 104 | for (int col = 0; col < nCols; col++) { |
1165 | 76 | ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) * |
1166 | data[row][col].hashCode(); |
|
1167 | } |
|
1168 | } |
|
1169 | 12 | return ret; |
1170 | } |
|
1171 | ||
1172 | //------------------------ Protected methods |
|
1173 | ||
1174 | /** |
|
1175 | * Returns <code>dimension x dimension</code> identity matrix. |
|
1176 | * |
|
1177 | * @param dimension dimension of identity matrix to generate |
|
1178 | * @return identity matrix |
|
1179 | * @throws IllegalArgumentException if dimension is not positive |
|
1180 | * @deprecated use {@link MatrixUtils#createBigIdentityMatrix} |
|
1181 | */ |
|
1182 | protected BigMatrix getIdentity(int dimension) { |
|
1183 | 0 | return MatrixUtils.createBigIdentityMatrix(dimension); |
1184 | } |
|
1185 | ||
1186 | /** |
|
1187 | * Returns the LU decomposition as a BigMatrix. |
|
1188 | * Returns a fresh copy of the cached LU matrix if this has been computed; |
|
1189 | * otherwise the composition is computed and cached for use by other methods. |
|
1190 | * Since a copy is returned in either case, changes to the returned matrix do not |
|
1191 | * affect the LU decomposition property. |
|
1192 | * <p> |
|
1193 | * The matrix returned is a compact representation of the LU decomposition. |
|
1194 | * Elements below the main diagonal correspond to entries of the "L" matrix; |
|
1195 | * elements on and above the main diagonal correspond to entries of the "U" |
|
1196 | * matrix. |
|
1197 | * <p> |
|
1198 | * Example: <pre> |
|
1199 | * |
|
1200 | * Returned matrix L U |
|
1201 | * 2 3 1 1 0 0 2 3 1 |
|
1202 | * 5 4 6 5 1 0 0 4 6 |
|
1203 | * 1 7 8 1 7 1 0 0 8 |
|
1204 | * </pre> |
|
1205 | * |
|
1206 | * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br> |
|
1207 | * where permuteRows reorders the rows of the matrix to follow the order determined |
|
1208 | * by the <a href=#getPermutation()>permutation</a> property. |
|
1209 | * |
|
1210 | * @return LU decomposition matrix |
|
1211 | * @throws InvalidMatrixException if the matrix is non-square or singular. |
|
1212 | */ |
|
1213 | protected BigMatrix getLUMatrix() throws InvalidMatrixException { |
|
1214 | 12 | if (lu == null) { |
1215 | 12 | luDecompose(); |
1216 | } |
|
1217 | 8 | return new BigMatrixImpl(lu); |
1218 | } |
|
1219 | ||
1220 | /** |
|
1221 | * Returns the permutation associated with the lu decomposition. |
|
1222 | * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1. |
|
1223 | * <p> |
|
1224 | * Example: |
|
1225 | * permutation = [1, 2, 0] means current 2nd row is first, current third row is second |
|
1226 | * and current first row is last. |
|
1227 | * <p> |
|
1228 | * Returns a fresh copy of the array. |
|
1229 | * |
|
1230 | * @return the permutation |
|
1231 | */ |
|
1232 | protected int[] getPermutation() { |
|
1233 | 8 | int[] out = new int[permutation.length]; |
1234 | 8 | System.arraycopy(permutation, 0, out, 0, permutation.length); |
1235 | 8 | return out; |
1236 | } |
|
1237 | ||
1238 | //------------------------ Private methods |
|
1239 | ||
1240 | /** |
|
1241 | * Returns a fresh copy of the underlying data array. |
|
1242 | * |
|
1243 | * @return a copy of the underlying data array. |
|
1244 | */ |
|
1245 | private BigDecimal[][] copyOut() { |
|
1246 | 52 | int nRows = this.getRowDimension(); |
1247 | 52 | BigDecimal[][] out = new BigDecimal[nRows][this.getColumnDimension()]; |
1248 | // can't copy 2-d array in one shot, otherwise get row references |
|
1249 | 208 | for (int i = 0; i < nRows; i++) { |
1250 | 156 | System.arraycopy(data[i], 0, out[i], 0, data[i].length); |
1251 | } |
|
1252 | 52 | return out; |
1253 | } |
|
1254 | ||
1255 | /** |
|
1256 | * Replaces data with a fresh copy of the input array. |
|
1257 | * <p> |
|
1258 | * Verifies that the input array is rectangular and non-empty. |
|
1259 | * |
|
1260 | * @param in data to copy in |
|
1261 | * @throws IllegalArgumentException if input array is emtpy or not |
|
1262 | * rectangular |
|
1263 | * @throws NullPointerException if input array is null |
|
1264 | */ |
|
1265 | private void copyIn(BigDecimal[][] in) { |
|
1266 | 164 | setSubMatrix(in,0,0); |
1267 | 164 | } |
1268 | ||
1269 | /** |
|
1270 | * Replaces data with a fresh copy of the input array. |
|
1271 | * |
|
1272 | * @param in data to copy in |
|
1273 | */ |
|
1274 | private void copyIn(double[][] in) { |
|
1275 | 196 | int nRows = in.length; |
1276 | 196 | int nCols = in[0].length; |
1277 | 196 | data = new BigDecimal[nRows][nCols]; |
1278 | 740 | for (int i = 0; i < nRows; i++) { |
1279 | 2136 | for (int j=0; j < nCols; j++) { |
1280 | 1592 | data[i][j] = new BigDecimal(in[i][j]); |
1281 | } |
|
1282 | } |
|
1283 | 196 | lu = null; |
1284 | 196 | } |
1285 | ||
1286 | /** |
|
1287 | * Replaces data with BigDecimals represented by the strings in the input |
|
1288 | * array. |
|
1289 | * |
|
1290 | * @param in data to copy in |
|
1291 | */ |
|
1292 | private void copyIn(String[][] in) { |
|
1293 | 20 | int nRows = in.length; |
1294 | 20 | int nCols = in[0].length; |
1295 | 20 | data = new BigDecimal[nRows][nCols]; |
1296 | 58 | for (int i = 0; i < nRows; i++) { |
1297 | 100 | for (int j=0; j < nCols; j++) { |
1298 | 62 | data[i][j] = new BigDecimal(in[i][j]); |
1299 | } |
|
1300 | } |
|
1301 | 18 | lu = null; |
1302 | 18 | } |
1303 | ||
1304 | /** |
|
1305 | * Tests a given coordinate as being valid or invalid |
|
1306 | * |
|
1307 | * @param row the row index. |
|
1308 | * @param col the column index. |
|
1309 | * @return true if the coordinate is with the current dimensions |
|
1310 | */ |
|
1311 | private boolean isValidCoordinate(int row, int col) { |
|
1312 | 1780 | int nRows = this.getRowDimension(); |
1313 | 1780 | int nCols = this.getColumnDimension(); |
1314 | ||
1315 | 1780 | return !(row < 0 || row >= nRows || col < 0 || col >= nCols); |
1316 | } |
|
1317 | ||
1318 | } |