001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.linear; 019 import java.io.Serializable; 020 import java.math.BigDecimal; 021 022 import org.apache.commons.math.MathRuntimeException; 023 024 /** 025 * Implementation of {@link BigMatrix} using a BigDecimal[][] array to store entries 026 * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf"> 027 * LU decompostion</a> to support linear system 028 * solution and inverse. 029 * <p> 030 * The LU decompostion is performed as needed, to support the following operations: <ul> 031 * <li>solve</li> 032 * <li>isSingular</li> 033 * <li>getDeterminant</li> 034 * <li>inverse</li> </ul></p> 035 * <p> 036 * <strong>Usage notes</strong>:<br> 037 * <ul><li> 038 * The LU decomposition is stored and reused on subsequent calls. If matrix 039 * data are modified using any of the public setXxx methods, the saved 040 * decomposition is discarded. If data are modified via references to the 041 * underlying array obtained using <code>getDataRef()</code>, then the stored 042 * LU decomposition will not be discarded. In this case, you need to 043 * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition 044 * before using any of the methods above.</li> 045 * <li> 046 * As specified in the {@link BigMatrix} interface, matrix element indexing 047 * is 0-based -- e.g., <code>getEntry(0, 0)</code> 048 * returns the element in the first row, first column of the matrix.</li></ul></p> 049 * 050 * @deprecated as of 2.0, replaced by {@link Array2DRowFieldMatrix} with a {@link 051 * org.apache.commons.math.util.BigReal} parameter 052 * @version $Revision: 811833 $ $Date: 2009-09-06 12:27:50 -0400 (Sun, 06 Sep 2009) $ 053 */ 054 @Deprecated 055 public class BigMatrixImpl implements BigMatrix, Serializable { 056 057 /** BigDecimal 0 */ 058 static final BigDecimal ZERO = new BigDecimal(0); 059 060 /** BigDecimal 1 */ 061 static final BigDecimal ONE = new BigDecimal(1); 062 063 /** Bound to determine effective singularity in LU decomposition */ 064 private static final BigDecimal TOO_SMALL = new BigDecimal(10E-12); 065 066 /** Serialization id */ 067 private static final long serialVersionUID = -1011428905656140431L; 068 069 /** Entries of the matrix */ 070 protected BigDecimal data[][] = null; 071 072 /** Entries of cached LU decomposition. 073 * All updates to data (other than luDecompose()) *must* set this to null 074 */ 075 protected BigDecimal lu[][] = null; 076 077 /** Permutation associated with LU decomposition */ 078 protected int[] permutation = null; 079 080 /** Parity of the permutation associated with the LU decomposition */ 081 protected int parity = 1; 082 083 /** Rounding mode for divisions **/ 084 private int roundingMode = BigDecimal.ROUND_HALF_UP; 085 086 /*** BigDecimal scale ***/ 087 private int scale = 64; 088 089 /** 090 * Creates a matrix with no data 091 */ 092 public BigMatrixImpl() { 093 } 094 095 /** 096 * Create a new BigMatrix with the supplied row and column dimensions. 097 * 098 * @param rowDimension the number of rows in the new matrix 099 * @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 }