0001: /*
0002: * Copyright 2003-2005 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 org.apache.commons.math.util.MathUtils;
0021:
0022: /**
0023: * Implementation of RealMatrix using a double[][] array to store entries and
0024: * <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 cached and reused on subsequent calls.
0037: * If data are modified via references to the underlying array obtained using
0038: * <code>getDataRef()</code>, then the stored LU decomposition will not be
0039: * discarded. In this case, you need to explicitly invoke
0040: * <code>LUDecompose()</code> to recompute the decomposition
0041: * before using any of the methods above.</li>
0042: * <li>
0043: * As specified in the {@link RealMatrix} interface, matrix element indexing
0044: * is 0-based -- e.g., <code>getEntry(0, 0)</code>
0045: * returns the element in the first row, first column of the matrix.</li></ul>
0046: *
0047: * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
0048: */
0049: public class RealMatrixImpl implements RealMatrix, Serializable {
0050:
0051: /** Serializable version identifier */
0052: private static final long serialVersionUID = 4237564493130426188L;
0053:
0054: /** Entries of the matrix */
0055: private double 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 double 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: /** Bound to determine effective singularity in LU decomposition */
0069: protected static double TOO_SMALL = 10E-12;
0070:
0071: /**
0072: * Creates a matrix with no data
0073: */
0074: public RealMatrixImpl() {
0075: }
0076:
0077: /**
0078: * Create a new RealMatrix with the supplied row and column dimensions.
0079: *
0080: * @param rowDimension the number of rows in the new matrix
0081: * @param columnDimension the number of columns in the new matrix
0082: * @throws IllegalArgumentException if row or column dimension is not
0083: * positive
0084: */
0085: public RealMatrixImpl(int rowDimension, int columnDimension) {
0086: if (rowDimension <= 0 || columnDimension <= 0) {
0087: throw new IllegalArgumentException(
0088: "row and column dimensions must be postive");
0089: }
0090: data = new double[rowDimension][columnDimension];
0091: lu = null;
0092: }
0093:
0094: /**
0095: * Create a new RealMatrix using the input array as the underlying
0096: * data array.
0097: * <p>
0098: * The input array is copied, not referenced.
0099: *
0100: * @param d data for new matrix
0101: * @throws IllegalArgumentException if <code>data</code> is not rectangular
0102: * (not all rows have the same length) or empty
0103: * @throws NullPointerException if <code>data</code> is null
0104: */
0105: public RealMatrixImpl(double[][] d) {
0106: this .copyIn(d);
0107: lu = null;
0108: }
0109:
0110: /**
0111: * Create a new (column) RealMatrix using <code>v</code> as the
0112: * data for the unique column of the <code>v.length x 1</code> matrix
0113: * created.
0114: * <p>
0115: * The input array is copied, not referenced.
0116: *
0117: * @param v column vector holding data for new matrix
0118: */
0119: public RealMatrixImpl(double[] v) {
0120: int nRows = v.length;
0121: data = new double[nRows][1];
0122: for (int row = 0; row < nRows; row++) {
0123: data[row][0] = v[row];
0124: }
0125: }
0126:
0127: /**
0128: * Create a new RealMatrix which is a copy of this.
0129: *
0130: * @return the cloned matrix
0131: */
0132: public RealMatrix copy() {
0133: return new RealMatrixImpl(this .copyOut());
0134: }
0135:
0136: /**
0137: * Compute the sum of this and <code>m</code>.
0138: *
0139: * @param m matrix to be added
0140: * @return this + m
0141: * @throws IllegalArgumentException if m is not the same size as this
0142: */
0143: public RealMatrix add(RealMatrix m) throws IllegalArgumentException {
0144: if (this .getColumnDimension() != m.getColumnDimension()
0145: || this .getRowDimension() != m.getRowDimension()) {
0146: throw new IllegalArgumentException(
0147: "matrix dimension mismatch");
0148: }
0149: int rowCount = this .getRowDimension();
0150: int columnCount = this .getColumnDimension();
0151: double[][] outData = new double[rowCount][columnCount];
0152: for (int row = 0; row < rowCount; row++) {
0153: for (int col = 0; col < columnCount; col++) {
0154: outData[row][col] = data[row][col]
0155: + m.getEntry(row, col);
0156: }
0157: }
0158: return new RealMatrixImpl(outData);
0159: }
0160:
0161: /**
0162: * Compute this minus <code>m</code>.
0163: *
0164: * @param m matrix to be subtracted
0165: * @return this + m
0166: * @throws IllegalArgumentException if m is not the same size as *this
0167: */
0168: public RealMatrix subtract(RealMatrix m)
0169: throws IllegalArgumentException {
0170: if (this .getColumnDimension() != m.getColumnDimension()
0171: || this .getRowDimension() != m.getRowDimension()) {
0172: throw new IllegalArgumentException(
0173: "matrix dimension mismatch");
0174: }
0175: int rowCount = this .getRowDimension();
0176: int columnCount = this .getColumnDimension();
0177: double[][] outData = new double[rowCount][columnCount];
0178: for (int row = 0; row < rowCount; row++) {
0179: for (int col = 0; col < columnCount; col++) {
0180: outData[row][col] = data[row][col]
0181: - m.getEntry(row, col);
0182: }
0183: }
0184: return new RealMatrixImpl(outData);
0185: }
0186:
0187: /**
0188: * Returns the result of adding d to each entry of this.
0189: *
0190: * @param d value to be added to each entry
0191: * @return d + this
0192: */
0193: public RealMatrix scalarAdd(double d) {
0194: int rowCount = this .getRowDimension();
0195: int columnCount = this .getColumnDimension();
0196: double[][] outData = new double[rowCount][columnCount];
0197: for (int row = 0; row < rowCount; row++) {
0198: for (int col = 0; col < columnCount; col++) {
0199: outData[row][col] = data[row][col] + d;
0200: }
0201: }
0202: return new RealMatrixImpl(outData);
0203: }
0204:
0205: /**
0206: * Returns the result multiplying each entry of this by <code>d</code>
0207: * @param d value to multiply all entries by
0208: * @return d * this
0209: */
0210: public RealMatrix scalarMultiply(double d) {
0211: int rowCount = this .getRowDimension();
0212: int columnCount = this .getColumnDimension();
0213: double[][] outData = new double[rowCount][columnCount];
0214: for (int row = 0; row < rowCount; row++) {
0215: for (int col = 0; col < columnCount; col++) {
0216: outData[row][col] = data[row][col] * d;
0217: }
0218: }
0219: return new RealMatrixImpl(outData);
0220: }
0221:
0222: /**
0223: * Returns the result of postmultiplying this by <code>m</code>.
0224: * @param m matrix to postmultiply by
0225: * @return this*m
0226: * @throws IllegalArgumentException
0227: * if columnDimension(this) != rowDimension(m)
0228: */
0229: public RealMatrix multiply(RealMatrix m)
0230: throws IllegalArgumentException {
0231: if (this .getColumnDimension() != m.getRowDimension()) {
0232: throw new IllegalArgumentException(
0233: "Matrices are not multiplication compatible.");
0234: }
0235: int nRows = this .getRowDimension();
0236: int nCols = m.getColumnDimension();
0237: int nSum = this .getColumnDimension();
0238: double[][] outData = new double[nRows][nCols];
0239: double sum = 0;
0240: for (int row = 0; row < nRows; row++) {
0241: for (int col = 0; col < nCols; col++) {
0242: sum = 0;
0243: for (int i = 0; i < nSum; i++) {
0244: sum += data[row][i] * m.getEntry(i, col);
0245: }
0246: outData[row][col] = sum;
0247: }
0248: }
0249: return new RealMatrixImpl(outData);
0250: }
0251:
0252: /**
0253: * Returns the result premultiplying this by <code>m</code>.
0254: * @param m matrix to premultiply by
0255: * @return m * this
0256: * @throws IllegalArgumentException
0257: * if rowDimension(this) != columnDimension(m)
0258: */
0259: public RealMatrix preMultiply(RealMatrix m)
0260: throws IllegalArgumentException {
0261: return m.multiply(this );
0262: }
0263:
0264: /**
0265: * Returns matrix entries as a two-dimensional array.
0266: * <p>
0267: * Makes a fresh copy of the underlying data.
0268: *
0269: * @return 2-dimensional array of entries
0270: */
0271: public double[][] getData() {
0272: return copyOut();
0273: }
0274:
0275: /**
0276: * Returns a reference to the underlying data array.
0277: * <p>
0278: * Does not make a fresh copy of the underlying data.
0279: *
0280: * @return 2-dimensional array of entries
0281: */
0282: public double[][] getDataRef() {
0283: return data;
0284: }
0285:
0286: /**
0287: *
0288: * @return norm
0289: */
0290: public double getNorm() {
0291: double maxColSum = 0;
0292: for (int col = 0; col < this .getColumnDimension(); col++) {
0293: double sum = 0;
0294: for (int row = 0; row < this .getRowDimension(); row++) {
0295: sum += Math.abs(data[row][col]);
0296: }
0297: maxColSum = Math.max(maxColSum, sum);
0298: }
0299: return maxColSum;
0300: }
0301:
0302: /**
0303: * Gets a submatrix. Rows and columns are indicated
0304: * counting from 0 to n-1.
0305: *
0306: * @param startRow Initial row index
0307: * @param endRow Final row index
0308: * @param startColumn Initial column index
0309: * @param endColumn Final column index
0310: * @return The subMatrix containing the data of the
0311: * specified rows and columns
0312: * @exception MatrixIndexException if row or column selections are not valid
0313: */
0314: public RealMatrix getSubMatrix(int startRow, int endRow,
0315: int startColumn, int endColumn) throws MatrixIndexException {
0316: if (startRow < 0 || startRow > endRow || endRow > data.length
0317: || startColumn < 0 || startColumn > endColumn
0318: || endColumn > data[0].length) {
0319: throw new MatrixIndexException(
0320: "invalid row or column index selection");
0321: }
0322: RealMatrixImpl subMatrix = new RealMatrixImpl(endRow - startRow
0323: + 1, endColumn - startColumn + 1);
0324: double[][] subMatrixData = subMatrix.getDataRef();
0325: for (int i = startRow; i <= endRow; i++) {
0326: for (int j = startColumn; j <= endColumn; j++) {
0327: subMatrixData[i - startRow][j - startColumn] = data[i][j];
0328: }
0329: }
0330: return subMatrix;
0331: }
0332:
0333: /**
0334: * Gets a submatrix. Rows and columns are indicated
0335: * counting from 0 to n-1.
0336: *
0337: * @param selectedRows Array of row indices must be non-empty
0338: * @param selectedColumns Array of column indices must be non-empty
0339: * @return The subMatrix containing the data in the
0340: * specified rows and columns
0341: * @exception MatrixIndexException if supplied row or column index arrays
0342: * are not valid
0343: */
0344: public RealMatrix getSubMatrix(int[] selectedRows,
0345: int[] selectedColumns) throws MatrixIndexException {
0346: if (selectedRows.length * selectedColumns.length == 0) {
0347: throw new MatrixIndexException(
0348: "selected row and column index arrays must be non-empty");
0349: }
0350: RealMatrixImpl subMatrix = new RealMatrixImpl(
0351: selectedRows.length, selectedColumns.length);
0352: double[][] subMatrixData = subMatrix.getDataRef();
0353: try {
0354: for (int i = 0; i < selectedRows.length; i++) {
0355: for (int j = 0; j < selectedColumns.length; j++) {
0356: subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
0357: }
0358: }
0359: } catch (ArrayIndexOutOfBoundsException e) {
0360: throw new MatrixIndexException("matrix dimension mismatch");
0361: }
0362: return subMatrix;
0363: }
0364:
0365: /**
0366: * Replace the submatrix starting at <code>row, column</code> using data in
0367: * the input <code>subMatrix</code> array. Indexes are 0-based.
0368: * <p>
0369: * Example:<br>
0370: * Starting with <pre>
0371: * 1 2 3 4
0372: * 5 6 7 8
0373: * 9 0 1 2
0374: * </pre>
0375: * and <code>subMatrix = {{3, 4} {5,6}}</code>, invoking
0376: * <code>setSubMatrix(subMatrix,1,1))</code> will result in <pre>
0377: * 1 2 3 4
0378: * 5 3 4 8
0379: * 9 5 6 2
0380: * </pre>
0381: *
0382: * @param subMatrix array containing the submatrix replacement data
0383: * @param row row coordinate of the top, left element to be replaced
0384: * @param column column coordinate of the top, left element to be replaced
0385: * @throws MatrixIndexException if subMatrix does not fit into this
0386: * matrix from element in (row, column)
0387: * @throws IllegalArgumentException if <code>subMatrix</code> is not rectangular
0388: * (not all rows have the same length) or empty
0389: * @throws NullPointerException if <code>subMatrix</code> is null
0390: * @since 1.1
0391: */
0392: public void setSubMatrix(double[][] subMatrix, int row, int column)
0393: throws MatrixIndexException {
0394: if ((row < 0) || (column < 0)) {
0395: throw new MatrixIndexException(
0396: "invalid row or column index selection");
0397: }
0398: int nRows = subMatrix.length;
0399: if (nRows == 0) {
0400: throw new IllegalArgumentException(
0401: "Matrix must have at least one row.");
0402: }
0403: int nCols = subMatrix[0].length;
0404: if (nCols == 0) {
0405: throw new IllegalArgumentException(
0406: "Matrix must have at least one column.");
0407: }
0408: for (int r = 1; r < nRows; r++) {
0409: if (subMatrix[r].length != nCols) {
0410: throw new IllegalArgumentException(
0411: "All input rows must have the same length.");
0412: }
0413: }
0414: if (data == null) {
0415: if ((row > 0) || (column > 0))
0416: throw new MatrixIndexException(
0417: "matrix must be initialized to perfom this method");
0418: data = new double[nRows][nCols];
0419: System.arraycopy(subMatrix, 0, data, 0, subMatrix.length);
0420: }
0421: if (((nRows + row) > this .getRowDimension())
0422: || (nCols + column > this .getColumnDimension()))
0423: throw new MatrixIndexException(
0424: "invalid row or column index selection");
0425: for (int i = 0; i < nRows; i++) {
0426: System.arraycopy(subMatrix[i], 0, data[row + i], column,
0427: nCols);
0428: }
0429: lu = null;
0430: }
0431:
0432: /**
0433: * Returns the entries in row number <code>row</code> as a row matrix.
0434: * Row indices start at 0.
0435: *
0436: * @param row the row to be fetched
0437: * @return row matrix
0438: * @throws MatrixIndexException if the specified row index is invalid
0439: */
0440: public RealMatrix getRowMatrix(int row) throws MatrixIndexException {
0441: if (!isValidCoordinate(row, 0)) {
0442: throw new MatrixIndexException("illegal row argument");
0443: }
0444: int ncols = this .getColumnDimension();
0445: double[][] out = new double[1][ncols];
0446: System.arraycopy(data[row], 0, out[0], 0, ncols);
0447: return new RealMatrixImpl(out);
0448: }
0449:
0450: /**
0451: * Returns the entries in column number <code>column</code>
0452: * as a column matrix. Column indices start at 0.
0453: *
0454: * @param column the column to be fetched
0455: * @return column matrix
0456: * @throws MatrixIndexException if the specified column index is invalid
0457: */
0458: public RealMatrix getColumnMatrix(int column)
0459: throws MatrixIndexException {
0460: if (!isValidCoordinate(0, column)) {
0461: throw new MatrixIndexException("illegal column argument");
0462: }
0463: int nRows = this .getRowDimension();
0464: double[][] out = new double[nRows][1];
0465: for (int row = 0; row < nRows; row++) {
0466: out[row][0] = data[row][column];
0467: }
0468: return new RealMatrixImpl(out);
0469: }
0470:
0471: /**
0472: * Returns the entries in row number <code>row</code> as an array.
0473: * <p>
0474: * Row indices start at 0. A <code>MatrixIndexException</code> is thrown
0475: * unless <code>0 <= row < rowDimension.</code>
0476: *
0477: * @param row the row to be fetched
0478: * @return array of entries in the row
0479: * @throws MatrixIndexException if the specified row index is not valid
0480: */
0481: public double[] getRow(int row) throws MatrixIndexException {
0482: if (!isValidCoordinate(row, 0)) {
0483: throw new MatrixIndexException("illegal row argument");
0484: }
0485: int ncols = this .getColumnDimension();
0486: double[] out = new double[ncols];
0487: System.arraycopy(data[row], 0, out, 0, ncols);
0488: return out;
0489: }
0490:
0491: /**
0492: * Returns the entries in column number <code>col</code> as an array.
0493: * <p>
0494: * Column indices start at 0. A <code>MatrixIndexException</code> is thrown
0495: * unless <code>0 <= column < columnDimension.</code>
0496: *
0497: * @param col the column to be fetched
0498: * @return array of entries in the column
0499: * @throws MatrixIndexException if the specified column index is not valid
0500: */
0501: public double[] getColumn(int col) throws MatrixIndexException {
0502: if (!isValidCoordinate(0, col)) {
0503: throw new MatrixIndexException("illegal column argument");
0504: }
0505: int nRows = this .getRowDimension();
0506: double[] out = new double[nRows];
0507: for (int row = 0; row < nRows; row++) {
0508: out[row] = data[row][col];
0509: }
0510: return out;
0511: }
0512:
0513: /**
0514: * Returns the entry in the specified row and column.
0515: * <p>
0516: * Row and column indices start at 0 and must satisfy
0517: * <ul>
0518: * <li><code>0 <= row < rowDimension</code></li>
0519: * <li><code> 0 <= column < columnDimension</code></li>
0520: * </ul>
0521: * otherwise a <code>MatrixIndexException</code> is thrown.
0522: *
0523: * @param row row location of entry to be fetched
0524: * @param column column location of entry to be fetched
0525: * @return matrix entry in row,column
0526: * @throws MatrixIndexException if the row or column index is not valid
0527: */
0528: public double getEntry(int row, int column)
0529: throws MatrixIndexException {
0530: if (!isValidCoordinate(row, column)) {
0531: throw new MatrixIndexException(
0532: "matrix entry does not exist");
0533: }
0534: return data[row][column];
0535: }
0536:
0537: /**
0538: * Returns the transpose matrix.
0539: *
0540: * @return transpose matrix
0541: */
0542: public RealMatrix transpose() {
0543: int nRows = this .getRowDimension();
0544: int nCols = this .getColumnDimension();
0545: RealMatrixImpl out = new RealMatrixImpl(nCols, nRows);
0546: double[][] outData = out.getDataRef();
0547: for (int row = 0; row < nRows; row++) {
0548: for (int col = 0; col < nCols; col++) {
0549: outData[col][row] = data[row][col];
0550: }
0551: }
0552: return out;
0553: }
0554:
0555: /**
0556: * Returns the inverse matrix if this matrix is invertible.
0557: *
0558: * @return inverse matrix
0559: * @throws InvalidMatrixException if this is not invertible
0560: */
0561: public RealMatrix inverse() throws InvalidMatrixException {
0562: return solve(MatrixUtils.createRealIdentityMatrix(this
0563: .getRowDimension()));
0564: }
0565:
0566: /**
0567: * @return determinant
0568: * @throws InvalidMatrixException if matrix is not square
0569: */
0570: public double getDeterminant() throws InvalidMatrixException {
0571: if (!isSquare()) {
0572: throw new InvalidMatrixException("matrix is not square");
0573: }
0574: if (isSingular()) { // note: this has side effect of attempting LU decomp if lu == null
0575: return 0d;
0576: } else {
0577: double det = parity;
0578: for (int i = 0; i < this .getRowDimension(); i++) {
0579: det *= lu[i][i];
0580: }
0581: return det;
0582: }
0583: }
0584:
0585: /**
0586: * @return true if the matrix is square (rowDimension = columnDimension)
0587: */
0588: public boolean isSquare() {
0589: return (this .getColumnDimension() == this .getRowDimension());
0590: }
0591:
0592: /**
0593: * @return true if the matrix is singular
0594: */
0595: public boolean isSingular() {
0596: if (lu == null) {
0597: try {
0598: luDecompose();
0599: return false;
0600: } catch (InvalidMatrixException ex) {
0601: return true;
0602: }
0603: } else { // LU decomp must have been successfully performed
0604: return false; // so the matrix is not singular
0605: }
0606: }
0607:
0608: /**
0609: * @return rowDimension
0610: */
0611: public int getRowDimension() {
0612: return data.length;
0613: }
0614:
0615: /**
0616: * @return columnDimension
0617: */
0618: public int getColumnDimension() {
0619: return data[0].length;
0620: }
0621:
0622: /**
0623: * @return trace
0624: * @throws IllegalArgumentException if the matrix is not square
0625: */
0626: public double getTrace() throws IllegalArgumentException {
0627: if (!isSquare()) {
0628: throw new IllegalArgumentException("matrix is not square");
0629: }
0630: double trace = data[0][0];
0631: for (int i = 1; i < this .getRowDimension(); i++) {
0632: trace += data[i][i];
0633: }
0634: return trace;
0635: }
0636:
0637: /**
0638: * @param v vector to operate on
0639: * @throws IllegalArgumentException if columnDimension != v.length
0640: * @return resulting vector
0641: */
0642: public double[] operate(double[] v) throws IllegalArgumentException {
0643: if (v.length != this .getColumnDimension()) {
0644: throw new IllegalArgumentException(
0645: "vector has wrong length");
0646: }
0647: int nRows = this .getRowDimension();
0648: int nCols = this .getColumnDimension();
0649: double[] out = new double[v.length];
0650: for (int row = 0; row < nRows; row++) {
0651: double sum = 0;
0652: for (int i = 0; i < nCols; i++) {
0653: sum += data[row][i] * v[i];
0654: }
0655: out[row] = sum;
0656: }
0657: return out;
0658: }
0659:
0660: /**
0661: * @param v vector to premultiply by
0662: * @throws IllegalArgumentException if rowDimension != v.length
0663: * @return resulting matrix
0664: */
0665: public double[] preMultiply(double[] v)
0666: throws IllegalArgumentException {
0667: int nRows = this .getRowDimension();
0668: if (v.length != nRows) {
0669: throw new IllegalArgumentException(
0670: "vector has wrong length");
0671: }
0672: int nCols = this .getColumnDimension();
0673: double[] out = new double[nCols];
0674: for (int col = 0; col < nCols; col++) {
0675: double sum = 0;
0676: for (int i = 0; i < nRows; i++) {
0677: sum += data[i][col] * v[i];
0678: }
0679: out[col] = sum;
0680: }
0681: return out;
0682: }
0683:
0684: /**
0685: * Returns a matrix of (column) solution vectors for linear systems with
0686: * coefficient matrix = this and constant vectors = columns of
0687: * <code>b</code>.
0688: *
0689: * @param b array of constant forming RHS of linear systems to
0690: * to solve
0691: * @return solution array
0692: * @throws IllegalArgumentException if this.rowDimension != row dimension
0693: * @throws InvalidMatrixException if this matrix is not square or is singular
0694: */
0695: public double[] solve(double[] b) throws IllegalArgumentException,
0696: InvalidMatrixException {
0697: int nRows = this .getRowDimension();
0698: if (b.length != nRows) {
0699: throw new IllegalArgumentException(
0700: "constant vector has wrong length");
0701: }
0702: RealMatrix bMatrix = new RealMatrixImpl(b);
0703: double[][] solution = ((RealMatrixImpl) (solve(bMatrix)))
0704: .getDataRef();
0705: double[] out = new double[nRows];
0706: for (int row = 0; row < nRows; row++) {
0707: out[row] = solution[row][0];
0708: }
0709: return out;
0710: }
0711:
0712: /**
0713: * Returns a matrix of (column) solution vectors for linear systems with
0714: * coefficient matrix = this and constant vectors = columns of
0715: * <code>b</code>.
0716: *
0717: * @param b matrix of constant vectors forming RHS of linear systems to
0718: * to solve
0719: * @return matrix of solution vectors
0720: * @throws IllegalArgumentException if this.rowDimension != row dimension
0721: * @throws InvalidMatrixException if this matrix is not square or is singular
0722: */
0723: public RealMatrix solve(RealMatrix b)
0724: throws IllegalArgumentException, InvalidMatrixException {
0725: if (b.getRowDimension() != this .getRowDimension()) {
0726: throw new IllegalArgumentException(
0727: "Incorrect row dimension");
0728: }
0729: if (!this .isSquare()) {
0730: throw new InvalidMatrixException(
0731: "coefficient matrix is not square");
0732: }
0733: if (this .isSingular()) { // side effect: compute LU decomp
0734: throw new InvalidMatrixException("Matrix is singular.");
0735: }
0736:
0737: int nCol = this .getColumnDimension();
0738: int nColB = b.getColumnDimension();
0739: int nRowB = b.getRowDimension();
0740:
0741: // Apply permutations to b
0742: double[][] bp = new double[nRowB][nColB];
0743: for (int row = 0; row < nRowB; row++) {
0744: for (int col = 0; col < nColB; col++) {
0745: bp[row][col] = b.getEntry(permutation[row], col);
0746: }
0747: }
0748:
0749: // Solve LY = b
0750: for (int col = 0; col < nCol; col++) {
0751: for (int i = col + 1; i < nCol; i++) {
0752: for (int j = 0; j < nColB; j++) {
0753: bp[i][j] -= bp[col][j] * lu[i][col];
0754: }
0755: }
0756: }
0757:
0758: // Solve UX = Y
0759: for (int col = nCol - 1; col >= 0; col--) {
0760: for (int j = 0; j < nColB; j++) {
0761: bp[col][j] /= lu[col][col];
0762: }
0763: for (int i = 0; i < col; i++) {
0764: for (int j = 0; j < nColB; j++) {
0765: bp[i][j] -= bp[col][j] * lu[i][col];
0766: }
0767: }
0768: }
0769:
0770: RealMatrixImpl outMat = new RealMatrixImpl(bp);
0771: return outMat;
0772: }
0773:
0774: /**
0775: * Computes a new
0776: * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
0777: * LU decompostion</a> for this matrix, storing the result for use by other methods.
0778: * <p>
0779: * <strong>Implementation Note</strong>:<br>
0780: * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
0781: * Crout's algortithm</a>, with partial pivoting.
0782: * <p>
0783: * <strong>Usage Note</strong>:<br>
0784: * This method should rarely be invoked directly. Its only use is
0785: * to force recomputation of the LU decomposition when changes have been
0786: * made to the underlying data using direct array references. Changes
0787: * made using setXxx methods will trigger recomputation when needed
0788: * automatically.
0789: *
0790: * @throws InvalidMatrixException if the matrix is non-square or singular.
0791: */
0792: public void luDecompose() throws InvalidMatrixException {
0793:
0794: int nRows = this .getRowDimension();
0795: int nCols = this .getColumnDimension();
0796: if (nRows != nCols) {
0797: throw new InvalidMatrixException(
0798: "LU decomposition requires that the matrix be square.");
0799: }
0800: lu = this .getData();
0801:
0802: // Initialize permutation array and parity
0803: permutation = new int[nRows];
0804: for (int row = 0; row < nRows; row++) {
0805: permutation[row] = row;
0806: }
0807: parity = 1;
0808:
0809: // Loop over columns
0810: for (int col = 0; col < nCols; col++) {
0811:
0812: double sum = 0;
0813:
0814: // upper
0815: for (int row = 0; row < col; row++) {
0816: sum = lu[row][col];
0817: for (int i = 0; i < row; i++) {
0818: sum -= lu[row][i] * lu[i][col];
0819: }
0820: lu[row][col] = sum;
0821: }
0822:
0823: // lower
0824: int max = col; // permutation row
0825: double largest = 0d;
0826: for (int row = col; row < nRows; row++) {
0827: sum = lu[row][col];
0828: for (int i = 0; i < col; i++) {
0829: sum -= lu[row][i] * lu[i][col];
0830: }
0831: lu[row][col] = sum;
0832:
0833: // maintain best permutation choice
0834: if (Math.abs(sum) > largest) {
0835: largest = Math.abs(sum);
0836: max = row;
0837: }
0838: }
0839:
0840: // Singularity check
0841: if (Math.abs(lu[max][col]) < TOO_SMALL) {
0842: lu = null;
0843: throw new InvalidMatrixException("matrix is singular");
0844: }
0845:
0846: // Pivot if necessary
0847: if (max != col) {
0848: double tmp = 0;
0849: for (int i = 0; i < nCols; i++) {
0850: tmp = lu[max][i];
0851: lu[max][i] = lu[col][i];
0852: lu[col][i] = tmp;
0853: }
0854: int temp = permutation[max];
0855: permutation[max] = permutation[col];
0856: permutation[col] = temp;
0857: parity = -parity;
0858: }
0859:
0860: //Divide the lower elements by the "winning" diagonal elt.
0861: for (int row = col + 1; row < nRows; row++) {
0862: lu[row][col] /= lu[col][col];
0863: }
0864: }
0865: }
0866:
0867: /**
0868: *
0869: * @see java.lang.Object#toString()
0870: */
0871: public String toString() {
0872: StringBuffer res = new StringBuffer();
0873: res.append("RealMatrixImpl{");
0874: if (data != null) {
0875: for (int i = 0; i < data.length; i++) {
0876: if (i > 0)
0877: res.append(",");
0878: res.append("{");
0879: for (int j = 0; j < data[0].length; j++) {
0880: if (j > 0)
0881: res.append(",");
0882: res.append(data[i][j]);
0883: }
0884: res.append("}");
0885: }
0886: }
0887: res.append("}");
0888: return res.toString();
0889: }
0890:
0891: /**
0892: * Returns true iff <code>object</code> is a
0893: * <code>RealMatrixImpl</code> instance with the same dimensions as this
0894: * and all corresponding matrix entries are equal. Corresponding entries
0895: * are compared using {@link java.lang.Double#doubleToLongBits(double)}
0896: *
0897: * @param object the object to test equality against.
0898: * @return true if object equals this
0899: */
0900: public boolean equals(Object object) {
0901: if (object == this ) {
0902: return true;
0903: }
0904: if (object instanceof RealMatrixImpl == false) {
0905: return false;
0906: }
0907: RealMatrix m = (RealMatrix) object;
0908: int nRows = getRowDimension();
0909: int nCols = getColumnDimension();
0910: if (m.getColumnDimension() != nCols
0911: || m.getRowDimension() != nRows) {
0912: return false;
0913: }
0914: for (int row = 0; row < nRows; row++) {
0915: for (int col = 0; col < nCols; col++) {
0916: if (Double.doubleToLongBits(data[row][col]) != Double
0917: .doubleToLongBits(m.getEntry(row, col))) {
0918: return false;
0919: }
0920: }
0921: }
0922: return true;
0923: }
0924:
0925: /**
0926: * Computes a hashcode for the matrix.
0927: *
0928: * @return hashcode for matrix
0929: */
0930: public int hashCode() {
0931: int ret = 7;
0932: int nRows = getRowDimension();
0933: int nCols = getColumnDimension();
0934: ret = ret * 31 + nRows;
0935: ret = ret * 31 + nCols;
0936: for (int row = 0; row < nRows; row++) {
0937: for (int col = 0; col < nCols; col++) {
0938: ret = ret * 31 + (11 * (row + 1) + 17 * (col + 1))
0939: * MathUtils.hash(data[row][col]);
0940: }
0941: }
0942: return ret;
0943: }
0944:
0945: //------------------------ Protected methods
0946:
0947: /**
0948: * Returns <code>dimension x dimension</code> identity matrix.
0949: *
0950: * @param dimension dimension of identity matrix to generate
0951: * @return identity matrix
0952: * @throws IllegalArgumentException if dimension is not positive
0953: * @deprecated use {@link MatrixUtils#createRealIdentityMatrix}
0954: */
0955: protected RealMatrix getIdentity(int dimension) {
0956: return MatrixUtils.createRealIdentityMatrix(dimension);
0957: }
0958:
0959: /**
0960: * Returns the LU decomposition as a RealMatrix.
0961: * Returns a fresh copy of the cached LU matrix if this has been computed;
0962: * otherwise the composition is computed and cached for use by other methods.
0963: * Since a copy is returned in either case, changes to the returned matrix do not
0964: * affect the LU decomposition property.
0965: * <p>
0966: * The matrix returned is a compact representation of the LU decomposition.
0967: * Elements below the main diagonal correspond to entries of the "L" matrix;
0968: * elements on and above the main diagonal correspond to entries of the "U"
0969: * matrix.
0970: * <p>
0971: * Example: <pre>
0972: *
0973: * Returned matrix L U
0974: * 2 3 1 1 0 0 2 3 1
0975: * 5 4 6 5 1 0 0 4 6
0976: * 1 7 8 1 7 1 0 0 8
0977: * </pre>
0978: *
0979: * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
0980: * where permuteRows reorders the rows of the matrix to follow the order determined
0981: * by the <a href=#getPermutation()>permutation</a> property.
0982: *
0983: * @return LU decomposition matrix
0984: * @throws InvalidMatrixException if the matrix is non-square or singular.
0985: */
0986: protected RealMatrix getLUMatrix() throws InvalidMatrixException {
0987: if (lu == null) {
0988: luDecompose();
0989: }
0990: return new RealMatrixImpl(lu);
0991: }
0992:
0993: /**
0994: * Returns the permutation associated with the lu decomposition.
0995: * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
0996: * <p>
0997: * Example:
0998: * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
0999: * and current first row is last.
1000: * <p>
1001: * Returns a fresh copy of the array.
1002: *
1003: * @return the permutation
1004: */
1005: protected int[] getPermutation() {
1006: int[] out = new int[permutation.length];
1007: System.arraycopy(permutation, 0, out, 0, permutation.length);
1008: return out;
1009: }
1010:
1011: //------------------------ Private methods
1012:
1013: /**
1014: * Returns a fresh copy of the underlying data array.
1015: *
1016: * @return a copy of the underlying data array.
1017: */
1018: private double[][] copyOut() {
1019: int nRows = this .getRowDimension();
1020: double[][] out = new double[nRows][this .getColumnDimension()];
1021: // can't copy 2-d array in one shot, otherwise get row references
1022: for (int i = 0; i < nRows; i++) {
1023: System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1024: }
1025: return out;
1026: }
1027:
1028: /**
1029: * Replaces data with a fresh copy of the input array.
1030: * <p>
1031: * Verifies that the input array is rectangular and non-empty
1032: *
1033: * @param in data to copy in
1034: * @throws IllegalArgumentException if input array is emtpy or not
1035: * rectangular
1036: * @throws NullPointerException if input array is null
1037: */
1038: private void copyIn(double[][] in) {
1039: setSubMatrix(in, 0, 0);
1040: }
1041:
1042: /**
1043: * Tests a given coordinate as being valid or invalid
1044: *
1045: * @param row the row index.
1046: * @param col the column index.
1047: * @return true if the coordinate is with the current dimensions
1048: */
1049: private boolean isValidCoordinate(int row, int col) {
1050: int nRows = this .getRowDimension();
1051: int nCols = this .getColumnDimension();
1052:
1053: return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols - 1);
1054: }
1055:
1056: }
|