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