0001: /*
0002: * Copyright 2004 The Apache Software Foundation.
0003: *
0004: * Licensed under the Apache License, Version 2.0 (the "License");
0005: * you may not use this file except in compliance with the License.
0006: * You may obtain a copy of the License at
0007: *
0008: * http://www.apache.org/licenses/LICENSE-2.0
0009: *
0010: * Unless required by applicable law or agreed to in writing, software
0011: * distributed under the License is distributed on an "AS IS" BASIS,
0012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0013: * See the License for the specific language governing permissions and
0014: * limitations under the License.
0015: */
0016:
0017: package org.apache.commons.math.linear;
0018:
0019: import java.io.Serializable;
0020: import java.math.BigDecimal;
0021:
0022: /**
0023: * Implementation of {@link BigMatrix} using a BigDecimal[][] array to store entries
0024: * and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
0025: * LU decompostion</a> to support linear system
0026: * solution and inverse.
0027: * <p>
0028: * The LU decompostion is performed as needed, to support the following operations: <ul>
0029: * <li>solve</li>
0030: * <li>isSingular</li>
0031: * <li>getDeterminant</li>
0032: * <li>inverse</li> </ul>
0033: * <p>
0034: * <strong>Usage notes</strong>:<br>
0035: * <ul><li>
0036: * The LU decomposition is stored and reused on subsequent calls. If matrix
0037: * data are modified using any of the public setXxx methods, the saved
0038: * decomposition is discarded. If data are modified via references to the
0039: * underlying array obtained using <code>getDataRef()</code>, then the stored
0040: * LU decomposition will not be discarded. In this case, you need to
0041: * explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
0042: * before using any of the methods above.</li>
0043: * <li>
0044: * As specified in the {@link BigMatrix} interface, matrix element indexing
0045: * is 0-based -- e.g., <code>getEntry(0, 0)</code>
0046: * returns the element in the first row, first column of the matrix.</li></ul>
0047: * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
0048: */
0049: public class BigMatrixImpl implements BigMatrix, Serializable {
0050:
0051: /** Serialization id */
0052: private static final long serialVersionUID = -1011428905656140431L;
0053:
0054: /** Entries of the matrix */
0055: private BigDecimal data[][] = null;
0056:
0057: /** Entries of cached LU decomposition.
0058: * All updates to data (other than luDecompose()) *must* set this to null
0059: */
0060: private BigDecimal lu[][] = null;
0061:
0062: /** Permutation associated with LU decompostion */
0063: private int[] permutation = null;
0064:
0065: /** Parity of the permutation associated with the LU decomposition */
0066: private int parity = 1;
0067:
0068: /** Rounding mode for divisions **/
0069: private int roundingMode = BigDecimal.ROUND_HALF_UP;
0070:
0071: /*** BigDecimal scale ***/
0072: private int scale = 64;
0073:
0074: /** Bound to determine effective singularity in LU decomposition */
0075: protected static BigDecimal TOO_SMALL = new BigDecimal(10E-12);
0076:
0077: /** BigDecimal 0 */
0078: static final BigDecimal ZERO = new BigDecimal(0);
0079: /** BigDecimal 1 */
0080: static final BigDecimal ONE = new BigDecimal(1);
0081:
0082: /**
0083: * Creates a matrix with no data
0084: */
0085: public BigMatrixImpl() {
0086: }
0087:
0088: /**
0089: * Create a new BigMatrix with the supplied row and column dimensions.
0090: *
0091: * @param rowDimension the number of rows in the new matrix
0092: * @param columnDimension the number of columns in the new matrix
0093: * @throws IllegalArgumentException if row or column dimension is not
0094: * positive
0095: */
0096: public BigMatrixImpl(int rowDimension, int columnDimension) {
0097: if (rowDimension <= 0 || columnDimension <= 0) {
0098: throw new IllegalArgumentException(
0099: "row and column dimensions must be positive");
0100: }
0101: data = new BigDecimal[rowDimension][columnDimension];
0102: lu = null;
0103: }
0104:
0105: /**
0106: * Create a new BigMatrix using the <code>data</code> as the underlying
0107: * data array.
0108: * <p>
0109: * The input array is copied, not referenced.
0110: *
0111: * @param d data for new matrix
0112: * @throws IllegalArgumentException if <code>d</code> is not rectangular
0113: * (not all rows have the same length) or empty
0114: * @throws NullPointerException if <code>d</code> is null
0115: */
0116: public BigMatrixImpl(BigDecimal[][] d) {
0117: this .copyIn(d);
0118: lu = null;
0119: }
0120:
0121: /**
0122: * Create a new BigMatrix using the <code>data</code> as the underlying
0123: * data array.
0124: * <p>
0125: * The input array is copied, not referenced.
0126: *
0127: * @param d data for new matrix
0128: * @throws IllegalArgumentException if <code>d</code> is not rectangular
0129: * (not all rows have the same length) or empty
0130: * @throws NullPointerException if <code>d</code> is null
0131: */
0132: public BigMatrixImpl(double[][] d) {
0133: int nRows = d.length;
0134: if (nRows == 0) {
0135: throw new IllegalArgumentException(
0136: "Matrix must have at least one row.");
0137: }
0138: int nCols = d[0].length;
0139: if (nCols == 0) {
0140: throw new IllegalArgumentException(
0141: "Matrix must have at least one column.");
0142: }
0143: for (int row = 1; row < nRows; row++) {
0144: if (d[row].length != nCols) {
0145: throw new IllegalArgumentException(
0146: "All input rows must have the same length.");
0147: }
0148: }
0149: this .copyIn(d);
0150: lu = null;
0151: }
0152:
0153: /**
0154: * Create a new BigMatrix using the values represented by the strings in
0155: * <code>data</code> as the underlying data array.
0156: *
0157: * @param d data for new matrix
0158: * @throws IllegalArgumentException if <code>d</code> is not rectangular
0159: * (not all rows have the same length) or empty
0160: * @throws NullPointerException if <code>d</code> is null
0161: */
0162: public BigMatrixImpl(String[][] d) {
0163: int nRows = d.length;
0164: if (nRows == 0) {
0165: throw new IllegalArgumentException(
0166: "Matrix must have at least one row.");
0167: }
0168: int nCols = d[0].length;
0169: if (nCols == 0) {
0170: throw new IllegalArgumentException(
0171: "Matrix must have at least one column.");
0172: }
0173: for (int row = 1; row < nRows; row++) {
0174: if (d[row].length != nCols) {
0175: throw new IllegalArgumentException(
0176: "All input rows must have the same length.");
0177: }
0178: }
0179: this .copyIn(d);
0180: lu = null;
0181: }
0182:
0183: /**
0184: * Create a new (column) BigMatrix using <code>v</code> as the
0185: * data for the unique column of the <code>v.length x 1</code> matrix
0186: * created.
0187: * <p>
0188: * The input array is copied, not referenced.
0189: *
0190: * @param v column vector holding data for new matrix
0191: */
0192: public BigMatrixImpl(BigDecimal[] v) {
0193: int nRows = v.length;
0194: data = new BigDecimal[nRows][1];
0195: for (int row = 0; row < nRows; row++) {
0196: data[row][0] = v[row];
0197: }
0198: }
0199:
0200: /**
0201: * Create a new BigMatrix which is a copy of this.
0202: *
0203: * @return the cloned matrix
0204: */
0205: public BigMatrix copy() {
0206: return new BigMatrixImpl(this .copyOut());
0207: }
0208:
0209: /**
0210: * Compute the sum of this and <code>m</code>.
0211: *
0212: * @param m matrix to be added
0213: * @return this + m
0214: * @exception IllegalArgumentException if m is not the same size as this
0215: */
0216: public BigMatrix add(BigMatrix m) throws IllegalArgumentException {
0217: if (this .getColumnDimension() != m.getColumnDimension()
0218: || this .getRowDimension() != m.getRowDimension()) {
0219: throw new IllegalArgumentException(
0220: "matrix dimension mismatch");
0221: }
0222: int rowCount = this .getRowDimension();
0223: int columnCount = this .getColumnDimension();
0224: BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
0225: for (int row = 0; row < rowCount; row++) {
0226: for (int col = 0; col < columnCount; col++) {
0227: outData[row][col] = data[row][col].add(m.getEntry(row,
0228: col));
0229: }
0230: }
0231: return new BigMatrixImpl(outData);
0232: }
0233:
0234: /**
0235: * Compute this minus <code>m</code>.
0236: *
0237: * @param m matrix to be subtracted
0238: * @return this + m
0239: * @exception IllegalArgumentException if m is not the same size as *this
0240: */
0241: public BigMatrix subtract(BigMatrix m)
0242: throws IllegalArgumentException {
0243: if (this .getColumnDimension() != m.getColumnDimension()
0244: || this .getRowDimension() != m.getRowDimension()) {
0245: throw new IllegalArgumentException(
0246: "matrix dimension mismatch");
0247: }
0248: int rowCount = this .getRowDimension();
0249: int columnCount = this .getColumnDimension();
0250: BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
0251: for (int row = 0; row < rowCount; row++) {
0252: for (int col = 0; col < columnCount; col++) {
0253: outData[row][col] = data[row][col].subtract(m.getEntry(
0254: row, col));
0255: }
0256: }
0257: return new BigMatrixImpl(outData);
0258: }
0259:
0260: /**
0261: * Returns the result of adding d to each entry of this.
0262: *
0263: * @param d value to be added to each entry
0264: * @return d + this
0265: */
0266: public BigMatrix scalarAdd(BigDecimal d) {
0267: int rowCount = this .getRowDimension();
0268: int columnCount = this .getColumnDimension();
0269: BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
0270: for (int row = 0; row < rowCount; row++) {
0271: for (int col = 0; col < columnCount; col++) {
0272: outData[row][col] = data[row][col].add(d);
0273: }
0274: }
0275: return new BigMatrixImpl(outData);
0276: }
0277:
0278: /**
0279: * Returns the result multiplying each entry of this by <code>d</code>
0280: * @param d value to multiply all entries by
0281: * @return d * this
0282: */
0283: public BigMatrix scalarMultiply(BigDecimal d) {
0284: int rowCount = this .getRowDimension();
0285: int columnCount = this .getColumnDimension();
0286: BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
0287: for (int row = 0; row < rowCount; row++) {
0288: for (int col = 0; col < columnCount; col++) {
0289: outData[row][col] = data[row][col].multiply(d);
0290: }
0291: }
0292: return new BigMatrixImpl(outData);
0293: }
0294:
0295: /**
0296: * Returns the result of postmultiplying this by <code>m</code>.
0297: * @param m matrix to postmultiply by
0298: * @return this*m
0299: * @throws IllegalArgumentException
0300: * if columnDimension(this) != rowDimension(m)
0301: */
0302: public BigMatrix multiply(BigMatrix m)
0303: throws IllegalArgumentException {
0304: if (this .getColumnDimension() != m.getRowDimension()) {
0305: throw new IllegalArgumentException(
0306: "Matrices are not multiplication compatible.");
0307: }
0308: int nRows = this .getRowDimension();
0309: int nCols = m.getColumnDimension();
0310: int nSum = this .getColumnDimension();
0311: BigDecimal[][] outData = new BigDecimal[nRows][nCols];
0312: BigDecimal sum = ZERO;
0313: for (int row = 0; row < nRows; row++) {
0314: for (int col = 0; col < nCols; col++) {
0315: sum = ZERO;
0316: for (int i = 0; i < nSum; i++) {
0317: sum = sum.add(data[row][i].multiply(m.getEntry(i,
0318: col)));
0319: }
0320: outData[row][col] = sum;
0321: }
0322: }
0323: return new BigMatrixImpl(outData);
0324: }
0325:
0326: /**
0327: * Returns the result premultiplying this by <code>m</code>.
0328: * @param m matrix to premultiply by
0329: * @return m * this
0330: * @throws IllegalArgumentException
0331: * if rowDimension(this) != columnDimension(m)
0332: */
0333: public BigMatrix preMultiply(BigMatrix m)
0334: throws IllegalArgumentException {
0335: return m.multiply(this );
0336: }
0337:
0338: /**
0339: * Returns matrix entries as a two-dimensional array.
0340: * <p>
0341: * Makes a fresh copy of the underlying data.
0342: *
0343: * @return 2-dimensional array of entries
0344: */
0345: public BigDecimal[][] getData() {
0346: return copyOut();
0347: }
0348:
0349: /**
0350: * Returns matrix entries as a two-dimensional array.
0351: * <p>
0352: * Makes a fresh copy of the underlying data converted to
0353: * <code>double</code> values.
0354: *
0355: * @return 2-dimensional array of entries
0356: */
0357: public double[][] getDataAsDoubleArray() {
0358: int nRows = getRowDimension();
0359: int nCols = getColumnDimension();
0360: double d[][] = new double[nRows][nCols];
0361: for (int i = 0; i < nRows; i++) {
0362: for (int j = 0; j < nCols; j++) {
0363: d[i][j] = data[i][j].doubleValue();
0364: }
0365: }
0366: return d;
0367: }
0368:
0369: /**
0370: * Returns a reference to the underlying data array.
0371: * <p>
0372: * Does not make a fresh copy of the underlying data.
0373: *
0374: * @return 2-dimensional array of entries
0375: */
0376: public BigDecimal[][] getDataRef() {
0377: return data;
0378: }
0379:
0380: /***
0381: * Gets the rounding mode for division operations
0382: * The default is {@link java.math.BigDecimal#ROUND_HALF_UP}
0383: * @see BigDecimal
0384: * @return the rounding mode.
0385: */
0386: public int getRoundingMode() {
0387: return roundingMode;
0388: }
0389:
0390: /***
0391: * Sets the rounding mode for decimal divisions.
0392: * @see BigDecimal
0393: * @param roundingMode
0394: */
0395: public void setRoundingMode(int roundingMode) {
0396: this .roundingMode = roundingMode;
0397: }
0398:
0399: /***
0400: * Sets the scale for division operations.
0401: * The default is 64
0402: * @see BigDecimal
0403: * @return the scale
0404: */
0405: public int getScale() {
0406: return scale;
0407: }
0408:
0409: /***
0410: * Sets the scale for division operations.
0411: * @see BigDecimal
0412: * @param scale
0413: */
0414: public void setScale(int scale) {
0415: this .scale = scale;
0416: }
0417:
0418: /**
0419: * Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html">
0420: * maximum absolute row sum norm</a> of the matrix.
0421: *
0422: * @return norm
0423: */
0424: public BigDecimal getNorm() {
0425: BigDecimal maxColSum = ZERO;
0426: for (int col = 0; col < this .getColumnDimension(); col++) {
0427: BigDecimal sum = ZERO;
0428: for (int row = 0; row < this .getRowDimension(); row++) {
0429: sum = sum.add(data[row][col].abs());
0430: }
0431: maxColSum = maxColSum.max(sum);
0432: }
0433: return maxColSum;
0434: }
0435:
0436: /**
0437: * Gets a submatrix. Rows and columns are indicated
0438: * counting from 0 to n-1.
0439: *
0440: * @param startRow Initial row index
0441: * @param endRow Final row index
0442: * @param startColumn Initial column index
0443: * @param endColumn Final column index
0444: * @return The subMatrix containing the data of the
0445: * specified rows and columns
0446: * @exception MatrixIndexException if row or column selections are not valid
0447: */
0448: public BigMatrix getSubMatrix(int startRow, int endRow,
0449: int startColumn, int endColumn) throws MatrixIndexException {
0450: if (startRow < 0 || startRow > endRow || endRow > data.length
0451: || startColumn < 0 || startColumn > endColumn
0452: || endColumn > data[0].length) {
0453: throw new MatrixIndexException(
0454: "invalid row or column index selection");
0455: }
0456: BigMatrixImpl subMatrix = new BigMatrixImpl(endRow - startRow
0457: + 1, endColumn - startColumn + 1);
0458: BigDecimal[][] subMatrixData = subMatrix.getDataRef();
0459: for (int i = startRow; i <= endRow; i++) {
0460: for (int j = startColumn; j <= endColumn; j++) {
0461: subMatrixData[i - startRow][j - startColumn] = data[i][j];
0462: }
0463: }
0464: return subMatrix;
0465: }
0466:
0467: /**
0468: * Gets a submatrix. Rows and columns are indicated
0469: * counting from 0 to n-1.
0470: *
0471: * @param selectedRows Array of row indices must be non-empty
0472: * @param selectedColumns Array of column indices must be non-empty
0473: * @return The subMatrix containing the data in the
0474: * specified rows and columns
0475: * @exception MatrixIndexException if supplied row or column index arrays
0476: * are not valid
0477: */
0478: public BigMatrix getSubMatrix(int[] selectedRows,
0479: int[] selectedColumns) throws MatrixIndexException {
0480: if (selectedRows.length * selectedColumns.length == 0) {
0481: throw new MatrixIndexException(
0482: "selected row and column index arrays must be non-empty");
0483: }
0484: BigMatrixImpl subMatrix = new BigMatrixImpl(
0485: selectedRows.length, selectedColumns.length);
0486: BigDecimal[][] subMatrixData = subMatrix.getDataRef();
0487: try {
0488: for (int i = 0; i < selectedRows.length; i++) {
0489: for (int j = 0; j < selectedColumns.length; j++) {
0490: subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
0491: }
0492: }
0493: } catch (ArrayIndexOutOfBoundsException e) {
0494: throw new MatrixIndexException("matrix dimension mismatch");
0495: }
0496: return subMatrix;
0497: }
0498:
0499: /**
0500: * Replace the submatrix starting at <code>row, column</code> using data in
0501: * the input <code>subMatrix</code> array. Indexes are 0-based.
0502: * <p>
0503: * Example:<br>
0504: * Starting with <pre>
0505: * 1 2 3 4
0506: * 5 6 7 8
0507: * 9 0 1 2
0508: * </pre>
0509: * and <code>subMatrix = {{3, 4} {5,6}}</code>, invoking
0510: * <code>setSubMatrix(subMatrix,1,1))</code> will result in <pre>
0511: * 1 2 3 4
0512: * 5 3 4 8
0513: * 9 5 6 2
0514: * </pre>
0515: *
0516: * @param subMatrix array containing the submatrix replacement data
0517: * @param row row coordinate of the top, left element to be replaced
0518: * @param column column coordinate of the top, left element to be replaced
0519: * @throws MatrixIndexException if subMatrix does not fit into this
0520: * matrix from element in (row, column)
0521: * @throws IllegalArgumentException if <code>subMatrix</code> is not rectangular
0522: * (not all rows have the same length) or empty
0523: * @throws NullPointerException if <code>subMatrix</code> is null
0524: * @since 1.1
0525: */
0526: public void setSubMatrix(BigDecimal[][] subMatrix, int row,
0527: int column) throws MatrixIndexException {
0528: if ((row < 0) || (column < 0)) {
0529: throw new MatrixIndexException(
0530: "invalid row or column index selection");
0531: }
0532: int nRows = subMatrix.length;
0533: if (nRows == 0) {
0534: throw new IllegalArgumentException(
0535: "Matrix must have at least one row.");
0536: }
0537: int nCols = subMatrix[0].length;
0538: if (nCols == 0) {
0539: throw new IllegalArgumentException(
0540: "Matrix must have at least one column.");
0541: }
0542: for (int r = 1; r < nRows; r++) {
0543: if (subMatrix[r].length != nCols) {
0544: throw new IllegalArgumentException(
0545: "All input rows must have the same length.");
0546: }
0547: }
0548: if (data == null) {
0549: if ((row > 0) || (column > 0))
0550: throw new MatrixIndexException(
0551: "matrix must be initialized to perfom this method");
0552: data = new BigDecimal[nRows][nCols];
0553: System.arraycopy(subMatrix, 0, data, 0, subMatrix.length);
0554: }
0555: if (((nRows + row) > this .getRowDimension())
0556: || (nCols + column > this .getColumnDimension()))
0557: throw new MatrixIndexException(
0558: "invalid row or column index selection");
0559: for (int i = 0; i < nRows; i++) {
0560: System.arraycopy(subMatrix[i], 0, data[row + i], column,
0561: nCols);
0562: }
0563: lu = null;
0564: }
0565:
0566: /**
0567: * Returns the entries in row number <code>row</code>
0568: * as a row matrix. Row indices start at 0.
0569: *
0570: * @param row the row to be fetched
0571: * @return row matrix
0572: * @throws MatrixIndexException if the specified row index is invalid
0573: */
0574: public BigMatrix getRowMatrix(int row) throws MatrixIndexException {
0575: if (!isValidCoordinate(row, 0)) {
0576: throw new MatrixIndexException("illegal row argument");
0577: }
0578: int ncols = this .getColumnDimension();
0579: BigDecimal[][] out = new BigDecimal[1][ncols];
0580: System.arraycopy(data[row], 0, out[0], 0, ncols);
0581: return new BigMatrixImpl(out);
0582: }
0583:
0584: /**
0585: * Returns the entries in column number <code>column</code>
0586: * as a column matrix. Column indices start at 0.
0587: *
0588: * @param column the column to be fetched
0589: * @return column matrix
0590: * @throws MatrixIndexException if the specified column index is invalid
0591: */
0592: public BigMatrix getColumnMatrix(int column)
0593: throws MatrixIndexException {
0594: if (!isValidCoordinate(0, column)) {
0595: throw new MatrixIndexException("illegal column argument");
0596: }
0597: int nRows = this .getRowDimension();
0598: BigDecimal[][] out = new BigDecimal[nRows][1];
0599: for (int row = 0; row < nRows; row++) {
0600: out[row][0] = data[row][column];
0601: }
0602: return new BigMatrixImpl(out);
0603: }
0604:
0605: /**
0606: * Returns the entries in row number <code>row</code> as an array.
0607: * <p>
0608: * Row indices start at 0. A <code>MatrixIndexException</code> is thrown
0609: * unless <code>0 <= row < rowDimension.</code>
0610: *
0611: * @param row the row to be fetched
0612: * @return array of entries in the row
0613: * @throws MatrixIndexException if the specified row index is not valid
0614: */
0615: public BigDecimal[] getRow(int row) throws MatrixIndexException {
0616: if (!isValidCoordinate(row, 0)) {
0617: throw new MatrixIndexException("illegal row argument");
0618: }
0619: int ncols = this .getColumnDimension();
0620: BigDecimal[] out = new BigDecimal[ncols];
0621: System.arraycopy(data[row], 0, out, 0, ncols);
0622: return out;
0623: }
0624:
0625: /**
0626: * Returns the entries in row number <code>row</code> as an array
0627: * of double values.
0628: * <p>
0629: * Row indices start at 0. A <code>MatrixIndexException</code> is thrown
0630: * unless <code>0 <= row < rowDimension.</code>
0631: *
0632: * @param row the row to be fetched
0633: * @return array of entries in the row
0634: * @throws MatrixIndexException if the specified row index is not valid
0635: */
0636: public double[] getRowAsDoubleArray(int row)
0637: throws MatrixIndexException {
0638: if (!isValidCoordinate(row, 0)) {
0639: throw new MatrixIndexException("illegal row argument");
0640: }
0641: int ncols = this .getColumnDimension();
0642: double[] out = new double[ncols];
0643: for (int i = 0; i < ncols; i++) {
0644: out[i] = data[row][i].doubleValue();
0645: }
0646: return out;
0647: }
0648:
0649: /**
0650: * Returns the entries in column number <code>col</code> as an array.
0651: * <p>
0652: * Column indices start at 0. A <code>MatrixIndexException</code> is thrown
0653: * unless <code>0 <= column < columnDimension.</code>
0654: *
0655: * @param col the column to be fetched
0656: * @return array of entries in the column
0657: * @throws MatrixIndexException if the specified column index is not valid
0658: */
0659: public BigDecimal[] getColumn(int col) throws MatrixIndexException {
0660: if (!isValidCoordinate(0, col)) {
0661: throw new MatrixIndexException("illegal column argument");
0662: }
0663: int nRows = this .getRowDimension();
0664: BigDecimal[] out = new BigDecimal[nRows];
0665: for (int i = 0; i < nRows; i++) {
0666: out[i] = data[i][col];
0667: }
0668: return out;
0669: }
0670:
0671: /**
0672: * Returns the entries in column number <code>col</code> as an array
0673: * of double values.
0674: * <p>
0675: * Column indices start at 0. A <code>MatrixIndexException</code> is thrown
0676: * unless <code>0 <= column < columnDimension.</code>
0677: *
0678: * @param col the column to be fetched
0679: * @return array of entries in the column
0680: * @throws MatrixIndexException if the specified column index is not valid
0681: */
0682: public double[] getColumnAsDoubleArray(int col)
0683: throws MatrixIndexException {
0684: if (!isValidCoordinate(0, col)) {
0685: throw new MatrixIndexException("illegal column argument");
0686: }
0687: int nrows = this .getRowDimension();
0688: double[] out = new double[nrows];
0689: for (int i = 0; i < nrows; i++) {
0690: out[i] = data[i][col].doubleValue();
0691: }
0692: return out;
0693: }
0694:
0695: /**
0696: * Returns the entry in the specified row and column.
0697: * <p>
0698: * Row and column indices start at 0 and must satisfy
0699: * <ul>
0700: * <li><code>0 <= row < rowDimension</code></li>
0701: * <li><code> 0 <= column < columnDimension</code></li>
0702: * </ul>
0703: * otherwise a <code>MatrixIndexException</code> is thrown.
0704: *
0705: * @param row row location of entry to be fetched
0706: * @param column column location of entry to be fetched
0707: * @return matrix entry in row,column
0708: * @throws MatrixIndexException if the row or column index is not valid
0709: */
0710: public BigDecimal getEntry(int row, int column)
0711: throws MatrixIndexException {
0712: if (!isValidCoordinate(row, column)) {
0713: throw new MatrixIndexException(
0714: "matrix entry does not exist");
0715: }
0716: return data[row][column];
0717: }
0718:
0719: /**
0720: * Returns the entry in the specified row and column as a double.
0721: * <p>
0722: * Row and column indices start at 0 and must satisfy
0723: * <ul>
0724: * <li><code>0 <= row < rowDimension</code></li>
0725: * <li><code> 0 <= column < columnDimension</code></li>
0726: * </ul>
0727: * otherwise a <code>MatrixIndexException</code> is thrown.
0728: *
0729: * @param row row location of entry to be fetched
0730: * @param column column location of entry to be fetched
0731: * @return matrix entry in row,column
0732: * @throws MatrixIndexException if the row
0733: * or column index is not valid
0734: */
0735: public double getEntryAsDouble(int row, int column)
0736: throws MatrixIndexException {
0737: return getEntry(row, column).doubleValue();
0738: }
0739:
0740: /**
0741: * Returns the transpose matrix.
0742: *
0743: * @return transpose matrix
0744: */
0745: public BigMatrix transpose() {
0746: int nRows = this .getRowDimension();
0747: int nCols = this .getColumnDimension();
0748: BigMatrixImpl out = new BigMatrixImpl(nCols, nRows);
0749: BigDecimal[][] outData = out.getDataRef();
0750: for (int row = 0; row < nRows; row++) {
0751: for (int col = 0; col < nCols; col++) {
0752: outData[col][row] = data[row][col];
0753: }
0754: }
0755: return out;
0756: }
0757:
0758: /**
0759: * Returns the inverse matrix if this matrix is invertible.
0760: *
0761: * @return inverse matrix
0762: * @throws InvalidMatrixException if this is not invertible
0763: */
0764: public BigMatrix inverse() throws InvalidMatrixException {
0765: return solve(MatrixUtils.createBigIdentityMatrix(this
0766: .getRowDimension()));
0767: }
0768:
0769: /**
0770: * Returns the determinant of this matrix.
0771: *
0772: * @return determinant
0773: * @throws InvalidMatrixException if matrix is not square
0774: */
0775: public BigDecimal getDeterminant() throws InvalidMatrixException {
0776: if (!isSquare()) {
0777: throw new InvalidMatrixException("matrix is not square");
0778: }
0779: if (isSingular()) { // note: this has side effect of attempting LU decomp if lu == null
0780: return ZERO;
0781: } else {
0782: BigDecimal det = (parity == 1) ? ONE : ONE.negate();
0783: for (int i = 0; i < this .getRowDimension(); i++) {
0784: det = det.multiply(lu[i][i]);
0785: }
0786: return det;
0787: }
0788: }
0789:
0790: /**
0791: * Is this a square matrix?
0792: * @return true if the matrix is square (rowDimension = columnDimension)
0793: */
0794: public boolean isSquare() {
0795: return (this .getColumnDimension() == this .getRowDimension());
0796: }
0797:
0798: /**
0799: * Is this a singular matrix?
0800: * @return true if the matrix is singular
0801: */
0802: public boolean isSingular() {
0803: if (lu == null) {
0804: try {
0805: luDecompose();
0806: return false;
0807: } catch (InvalidMatrixException ex) {
0808: return true;
0809: }
0810: } else { // LU decomp must have been successfully performed
0811: return false; // so the matrix is not singular
0812: }
0813: }
0814:
0815: /**
0816: * Returns the number of rows in the matrix.
0817: *
0818: * @return rowDimension
0819: */
0820: public int getRowDimension() {
0821: return data.length;
0822: }
0823:
0824: /**
0825: * Returns the number of columns in the matrix.
0826: *
0827: * @return columnDimension
0828: */
0829: public int getColumnDimension() {
0830: return data[0].length;
0831: }
0832:
0833: /**
0834: * Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
0835: * trace</a> of the matrix (the sum of the elements on the main diagonal).
0836: *
0837: * @return trace
0838: *
0839: * @throws IllegalArgumentException if this matrix is not square.
0840: */
0841: public BigDecimal getTrace() throws IllegalArgumentException {
0842: if (!isSquare()) {
0843: throw new IllegalArgumentException("matrix is not square");
0844: }
0845: BigDecimal trace = data[0][0];
0846: for (int i = 1; i < this .getRowDimension(); i++) {
0847: trace = trace.add(data[i][i]);
0848: }
0849: return trace;
0850: }
0851:
0852: /**
0853: * Returns the result of multiplying this by the vector <code>v</code>.
0854: *
0855: * @param v the vector to operate on
0856: * @return this*v
0857: * @throws IllegalArgumentException if columnDimension != v.size()
0858: */
0859: public BigDecimal[] operate(BigDecimal[] v)
0860: throws IllegalArgumentException {
0861: if (v.length != this .getColumnDimension()) {
0862: throw new IllegalArgumentException(
0863: "vector has wrong length");
0864: }
0865: int nRows = this .getRowDimension();
0866: int nCols = this .getColumnDimension();
0867: BigDecimal[] out = new BigDecimal[v.length];
0868: for (int row = 0; row < nRows; row++) {
0869: BigDecimal sum = ZERO;
0870: for (int i = 0; i < nCols; i++) {
0871: sum = sum.add(data[row][i].multiply(v[i]));
0872: }
0873: out[row] = sum;
0874: }
0875: return out;
0876: }
0877:
0878: /**
0879: * Returns the result of multiplying this by the vector <code>v</code>.
0880: *
0881: * @param v the vector to operate on
0882: * @return this*v
0883: * @throws IllegalArgumentException if columnDimension != v.size()
0884: */
0885: public BigDecimal[] operate(double[] v)
0886: throws IllegalArgumentException {
0887: BigDecimal bd[] = new BigDecimal[v.length];
0888: for (int i = 0; i < bd.length; i++) {
0889: bd[i] = new BigDecimal(v[i]);
0890: }
0891: return operate(bd);
0892: }
0893:
0894: /**
0895: * Returns the (row) vector result of premultiplying this by the vector <code>v</code>.
0896: *
0897: * @param v the row vector to premultiply by
0898: * @return v*this
0899: * @throws IllegalArgumentException if rowDimension != v.size()
0900: */
0901: public BigDecimal[] preMultiply(BigDecimal[] v)
0902: throws IllegalArgumentException {
0903: int nRows = this .getRowDimension();
0904: if (v.length != nRows) {
0905: throw new IllegalArgumentException(
0906: "vector has wrong length");
0907: }
0908: int nCols = this .getColumnDimension();
0909: BigDecimal[] out = new BigDecimal[nCols];
0910: for (int col = 0; col < nCols; col++) {
0911: BigDecimal sum = ZERO;
0912: for (int i = 0; i < nRows; i++) {
0913: sum = sum.add(data[i][col].multiply(v[i]));
0914: }
0915: out[col] = sum;
0916: }
0917: return out;
0918: }
0919:
0920: /**
0921: * Returns a matrix of (column) solution vectors for linear systems with
0922: * coefficient matrix = this and constant vectors = columns of
0923: * <code>b</code>.
0924: *
0925: * @param b array of constants forming RHS of linear systems to
0926: * to solve
0927: * @return solution array
0928: * @throws IllegalArgumentException if this.rowDimension != row dimension
0929: * @throws InvalidMatrixException if this matrix is not square or is singular
0930: */
0931: public BigDecimal[] solve(BigDecimal[] b)
0932: throws IllegalArgumentException, InvalidMatrixException {
0933: int nRows = this .getRowDimension();
0934: if (b.length != nRows) {
0935: throw new IllegalArgumentException(
0936: "constant vector has wrong length");
0937: }
0938: BigMatrix bMatrix = new BigMatrixImpl(b);
0939: BigDecimal[][] solution = ((BigMatrixImpl) (solve(bMatrix)))
0940: .getDataRef();
0941: BigDecimal[] out = new BigDecimal[nRows];
0942: for (int row = 0; row < nRows; row++) {
0943: out[row] = solution[row][0];
0944: }
0945: return out;
0946: }
0947:
0948: /**
0949: * Returns a matrix of (column) solution vectors for linear systems with
0950: * coefficient matrix = this and constant vectors = columns of
0951: * <code>b</code>.
0952: *
0953: * @param b array of constants forming RHS of linear systems to
0954: * to solve
0955: * @return solution array
0956: * @throws IllegalArgumentException if this.rowDimension != row dimension
0957: * @throws InvalidMatrixException if this matrix is not square or is singular
0958: */
0959: public BigDecimal[] solve(double[] b)
0960: throws IllegalArgumentException, InvalidMatrixException {
0961: BigDecimal bd[] = new BigDecimal[b.length];
0962: for (int i = 0; i < bd.length; i++) {
0963: bd[i] = new BigDecimal(b[i]);
0964: }
0965: return solve(bd);
0966: }
0967:
0968: /**
0969: * Returns a matrix of (column) solution vectors for linear systems with
0970: * coefficient matrix = this and constant vectors = columns of
0971: * <code>b</code>.
0972: *
0973: * @param b matrix of constant vectors forming RHS of linear systems to
0974: * to solve
0975: * @return matrix of solution vectors
0976: * @throws IllegalArgumentException if this.rowDimension != row dimension
0977: * @throws InvalidMatrixException if this matrix is not square or is singular
0978: */
0979: public BigMatrix solve(BigMatrix b)
0980: throws IllegalArgumentException, InvalidMatrixException {
0981: if (b.getRowDimension() != this .getRowDimension()) {
0982: throw new IllegalArgumentException(
0983: "Incorrect row dimension");
0984: }
0985: if (!this .isSquare()) {
0986: throw new InvalidMatrixException(
0987: "coefficient matrix is not square");
0988: }
0989: if (this .isSingular()) { // side effect: compute LU decomp
0990: throw new InvalidMatrixException("Matrix is singular.");
0991: }
0992:
0993: int nCol = this .getColumnDimension();
0994: int nColB = b.getColumnDimension();
0995: int nRowB = b.getRowDimension();
0996:
0997: // Apply permutations to b
0998: BigDecimal[][] bp = new BigDecimal[nRowB][nColB];
0999: for (int row = 0; row < nRowB; row++) {
1000: for (int col = 0; col < nColB; col++) {
1001: bp[row][col] = b.getEntry(permutation[row], col);
1002: }
1003: }
1004:
1005: // Solve LY = b
1006: for (int col = 0; col < nCol; col++) {
1007: for (int i = col + 1; i < nCol; i++) {
1008: for (int j = 0; j < nColB; j++) {
1009: bp[i][j] = bp[i][j].subtract(bp[col][j]
1010: .multiply(lu[i][col]));
1011: }
1012: }
1013: }
1014:
1015: // Solve UX = Y
1016: for (int col = nCol - 1; col >= 0; col--) {
1017: for (int j = 0; j < nColB; j++) {
1018: bp[col][j] = bp[col][j].divide(lu[col][col], scale,
1019: roundingMode);
1020: }
1021: for (int i = 0; i < col; i++) {
1022: for (int j = 0; j < nColB; j++) {
1023: bp[i][j] = bp[i][j].subtract(bp[col][j]
1024: .multiply(lu[i][col]));
1025: }
1026: }
1027: }
1028:
1029: BigMatrixImpl outMat = new BigMatrixImpl(bp);
1030: return outMat;
1031: }
1032:
1033: /**
1034: * Computes a new
1035: * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
1036: * LU decompostion</a> for this matrix, storing the result for use by other methods.
1037: * <p>
1038: * <strong>Implementation Note</strong>:<br>
1039: * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
1040: * Crout's algortithm</a>, with partial pivoting.
1041: * <p>
1042: * <strong>Usage Note</strong>:<br>
1043: * This method should rarely be invoked directly. Its only use is
1044: * to force recomputation of the LU decomposition when changes have been
1045: * made to the underlying data using direct array references. Changes
1046: * made using setXxx methods will trigger recomputation when needed
1047: * automatically.
1048: *
1049: * @throws InvalidMatrixException if the matrix is non-square or singular.
1050: */
1051: public void luDecompose() throws InvalidMatrixException {
1052:
1053: int nRows = this .getRowDimension();
1054: int nCols = this .getColumnDimension();
1055: if (nRows != nCols) {
1056: throw new InvalidMatrixException(
1057: "LU decomposition requires that the matrix be square.");
1058: }
1059: lu = this .getData();
1060:
1061: // Initialize permutation array and parity
1062: permutation = new int[nRows];
1063: for (int row = 0; row < nRows; row++) {
1064: permutation[row] = row;
1065: }
1066: parity = 1;
1067:
1068: // Loop over columns
1069: for (int col = 0; col < nCols; col++) {
1070:
1071: BigDecimal sum = ZERO;
1072:
1073: // upper
1074: for (int row = 0; row < col; row++) {
1075: sum = lu[row][col];
1076: for (int i = 0; i < row; i++) {
1077: sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
1078: }
1079: lu[row][col] = sum;
1080: }
1081:
1082: // lower
1083: int max = col; // permutation row
1084: BigDecimal largest = ZERO;
1085: for (int row = col; row < nRows; row++) {
1086: sum = lu[row][col];
1087: for (int i = 0; i < col; i++) {
1088: sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
1089: }
1090: lu[row][col] = sum;
1091:
1092: // maintain best permutation choice
1093: if (sum.abs().compareTo(largest) == 1) {
1094: largest = sum.abs();
1095: max = row;
1096: }
1097: }
1098:
1099: // Singularity check
1100: if (lu[max][col].abs().compareTo(TOO_SMALL) <= 0) {
1101: lu = null;
1102: throw new InvalidMatrixException("matrix is singular");
1103: }
1104:
1105: // Pivot if necessary
1106: if (max != col) {
1107: BigDecimal tmp = ZERO;
1108: for (int i = 0; i < nCols; i++) {
1109: tmp = lu[max][i];
1110: lu[max][i] = lu[col][i];
1111: lu[col][i] = tmp;
1112: }
1113: int temp = permutation[max];
1114: permutation[max] = permutation[col];
1115: permutation[col] = temp;
1116: parity = -parity;
1117: }
1118:
1119: //Divide the lower elements by the "winning" diagonal elt.
1120: for (int row = col + 1; row < nRows; row++) {
1121: lu[row][col] = lu[row][col].divide(lu[col][col], scale,
1122: roundingMode);
1123: }
1124:
1125: }
1126:
1127: }
1128:
1129: /**
1130: *
1131: * @see Object#toString()
1132: */
1133: public String toString() {
1134: StringBuffer res = new StringBuffer();
1135: res.append("BigMatrixImpl{");
1136: if (data != null) {
1137: for (int i = 0; i < data.length; i++) {
1138: if (i > 0)
1139: res.append(",");
1140: res.append("{");
1141: for (int j = 0; j < data[0].length; j++) {
1142: if (j > 0)
1143: res.append(",");
1144: res.append(data[i][j]);
1145: }
1146: res.append("}");
1147: }
1148: }
1149: res.append("}");
1150: return res.toString();
1151: }
1152:
1153: /**
1154: * Returns true iff <code>object</code> is a
1155: * <code>BigMatrixImpl</code> instance with the same dimensions as this
1156: * and all corresponding matrix entries are equal. BigDecimal.equals
1157: * is used to compare corresponding entries.
1158: *
1159: * @param object the object to test equality against.
1160: * @return true if object equals this
1161: */
1162: public boolean equals(Object object) {
1163: if (object == this ) {
1164: return true;
1165: }
1166: if (object instanceof BigMatrixImpl == false) {
1167: return false;
1168: }
1169: BigMatrix m = (BigMatrix) object;
1170: int nRows = getRowDimension();
1171: int nCols = getColumnDimension();
1172: if (m.getColumnDimension() != nCols
1173: || m.getRowDimension() != nRows) {
1174: return false;
1175: }
1176: for (int row = 0; row < nRows; row++) {
1177: for (int col = 0; col < nCols; col++) {
1178: if (!data[row][col].equals(m.getEntry(row, col))) {
1179: return false;
1180: }
1181: }
1182: }
1183: return true;
1184: }
1185:
1186: /**
1187: * Computes a hashcode for the matrix.
1188: *
1189: * @return hashcode for matrix
1190: */
1191: public int hashCode() {
1192: int ret = 7;
1193: int nRows = getRowDimension();
1194: int nCols = getColumnDimension();
1195: ret = ret * 31 + nRows;
1196: ret = ret * 31 + nCols;
1197: for (int row = 0; row < nRows; row++) {
1198: for (int col = 0; col < nCols; col++) {
1199: ret = ret * 31 + (11 * (row + 1) + 17 * (col + 1))
1200: * data[row][col].hashCode();
1201: }
1202: }
1203: return ret;
1204: }
1205:
1206: //------------------------ Protected methods
1207:
1208: /**
1209: * Returns <code>dimension x dimension</code> identity matrix.
1210: *
1211: * @param dimension dimension of identity matrix to generate
1212: * @return identity matrix
1213: * @throws IllegalArgumentException if dimension is not positive
1214: * @deprecated use {@link MatrixUtils#createBigIdentityMatrix}
1215: */
1216: protected BigMatrix getIdentity(int dimension) {
1217: return MatrixUtils.createBigIdentityMatrix(dimension);
1218: }
1219:
1220: /**
1221: * Returns the LU decomposition as a BigMatrix.
1222: * Returns a fresh copy of the cached LU matrix if this has been computed;
1223: * otherwise the composition is computed and cached for use by other methods.
1224: * Since a copy is returned in either case, changes to the returned matrix do not
1225: * affect the LU decomposition property.
1226: * <p>
1227: * The matrix returned is a compact representation of the LU decomposition.
1228: * Elements below the main diagonal correspond to entries of the "L" matrix;
1229: * elements on and above the main diagonal correspond to entries of the "U"
1230: * matrix.
1231: * <p>
1232: * Example: <pre>
1233: *
1234: * Returned matrix L U
1235: * 2 3 1 1 0 0 2 3 1
1236: * 5 4 6 5 1 0 0 4 6
1237: * 1 7 8 1 7 1 0 0 8
1238: * </pre>
1239: *
1240: * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
1241: * where permuteRows reorders the rows of the matrix to follow the order determined
1242: * by the <a href=#getPermutation()>permutation</a> property.
1243: *
1244: * @return LU decomposition matrix
1245: * @throws InvalidMatrixException if the matrix is non-square or singular.
1246: */
1247: protected BigMatrix getLUMatrix() throws InvalidMatrixException {
1248: if (lu == null) {
1249: luDecompose();
1250: }
1251: return new BigMatrixImpl(lu);
1252: }
1253:
1254: /**
1255: * Returns the permutation associated with the lu decomposition.
1256: * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
1257: * <p>
1258: * Example:
1259: * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
1260: * and current first row is last.
1261: * <p>
1262: * Returns a fresh copy of the array.
1263: *
1264: * @return the permutation
1265: */
1266: protected int[] getPermutation() {
1267: int[] out = new int[permutation.length];
1268: System.arraycopy(permutation, 0, out, 0, permutation.length);
1269: return out;
1270: }
1271:
1272: //------------------------ Private methods
1273:
1274: /**
1275: * Returns a fresh copy of the underlying data array.
1276: *
1277: * @return a copy of the underlying data array.
1278: */
1279: private BigDecimal[][] copyOut() {
1280: int nRows = this .getRowDimension();
1281: BigDecimal[][] out = new BigDecimal[nRows][this
1282: .getColumnDimension()];
1283: // can't copy 2-d array in one shot, otherwise get row references
1284: for (int i = 0; i < nRows; i++) {
1285: System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1286: }
1287: return out;
1288: }
1289:
1290: /**
1291: * Replaces data with a fresh copy of the input array.
1292: * <p>
1293: * Verifies that the input array is rectangular and non-empty.
1294: *
1295: * @param in data to copy in
1296: * @throws IllegalArgumentException if input array is emtpy or not
1297: * rectangular
1298: * @throws NullPointerException if input array is null
1299: */
1300: private void copyIn(BigDecimal[][] in) {
1301: setSubMatrix(in, 0, 0);
1302: }
1303:
1304: /**
1305: * Replaces data with a fresh copy of the input array.
1306: *
1307: * @param in data to copy in
1308: */
1309: private void copyIn(double[][] in) {
1310: int nRows = in.length;
1311: int nCols = in[0].length;
1312: data = new BigDecimal[nRows][nCols];
1313: for (int i = 0; i < nRows; i++) {
1314: for (int j = 0; j < nCols; j++) {
1315: data[i][j] = new BigDecimal(in[i][j]);
1316: }
1317: }
1318: lu = null;
1319: }
1320:
1321: /**
1322: * Replaces data with BigDecimals represented by the strings in the input
1323: * array.
1324: *
1325: * @param in data to copy in
1326: */
1327: private void copyIn(String[][] in) {
1328: int nRows = in.length;
1329: int nCols = in[0].length;
1330: data = new BigDecimal[nRows][nCols];
1331: for (int i = 0; i < nRows; i++) {
1332: for (int j = 0; j < nCols; j++) {
1333: data[i][j] = new BigDecimal(in[i][j]);
1334: }
1335: }
1336: lu = null;
1337: }
1338:
1339: /**
1340: * Tests a given coordinate as being valid or invalid
1341: *
1342: * @param row the row index.
1343: * @param col the column index.
1344: * @return true if the coordinate is with the current dimensions
1345: */
1346: private boolean isValidCoordinate(int row, int col) {
1347: int nRows = this .getRowDimension();
1348: int nCols = this .getColumnDimension();
1349:
1350: return !(row < 0 || row >= nRows || col < 0 || col >= nCols);
1351: }
1352:
1353: }
|