You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by ps...@apache.org on 2004/06/06 06:20:45 UTC
cvs commit: jakarta-commons/math/src/test/org/apache/commons/math/linear BigMatrixImplTest.java
psteitz 2004/06/05 21:20:45
Added: math/src/java/org/apache/commons/math/linear BigMatrix.java
BigMatrixImpl.java
math/src/test/org/apache/commons/math/linear
BigMatrixImplTest.java
Log:
Initial Commit of BigMatrix classes.
PR #28819
Submitted by Matthew Inger
Revision Changes Path
1.1 jakarta-commons/math/src/java/org/apache/commons/math/linear/BigMatrix.java
Index: BigMatrix.java
===================================================================
/*
* Copyright 2003-2004 The Apache Software Foundation.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import java.math.BigDecimal;
/**
* Interface defining a real-valued matrix with basic algebraic operations, using
* BigDecimal representations for the entries.
*
* @version $Revision: 1.1 $ $Date: 2004/06/06 04:20:45 $
*/
public interface BigMatrix {
/**
* Returns a (deep) copy of this.
*
* @return matrix copy
*/
BigMatrix copy();
/**
* Compute the sum of this and m.
*
* @param m matrix to be added
* @return this + m
* @exception IllegalArgumentException if m is not the same size as this
*/
BigMatrix add(BigMatrix m) throws IllegalArgumentException;
/**
* Compute this minus m.
*
* @param m matrix to be subtracted
* @return this + m
* @exception IllegalArgumentException if m is not the same size as this
*/
BigMatrix subtract(BigMatrix m) throws IllegalArgumentException;
/**
* Returns the result of adding d to each entry of this.
*
* @param d value to be added to each entry
* @return d + this
*/
BigMatrix scalarAdd(BigDecimal d);
/**
* Returns the result multiplying each entry of this by d.
*
* @param d value to multiply all entries by
* @return d * this
*/
BigMatrix scalarMultiply(BigDecimal d);
/**
* Returns the result postmultiplying this by m.
*
* @param m matrix to postmultiply by
* @return this * m
* @throws IllegalArgumentException
* if columnDimension(this) != rowDimension(m)
*/
BigMatrix multiply(BigMatrix m) throws IllegalArgumentException;
/**
* Returns the result premultiplying this by <code>m</code>.
* @param m matrix to premultiply by
* @return m * this
* @throws IllegalArgumentException
* if rowDimension(this) != columnDimension(m)
*/
public BigMatrix preMultiply(BigMatrix m) throws IllegalArgumentException;
/**
* Returns matrix entries as a two-dimensional array.
*
* @return 2-dimensional array of entries
*/
BigDecimal[][] getData();
/**
* Returns matrix entries as a two-dimensional array.
*
* @return 2-dimensional array of entries
*/
double [][] getDataAsDoubleArray();
/**
* Overwrites the underlying data for the matrix with
* a fresh copy of <code>data</code>.
*
* @param data 2-dimensional array of entries
*/
void setData(BigDecimal[][] data);
/**
* Overwrites the underlying data for the matrix with
* a fresh copy of <code>data</code>.
*
* @param data 2-dimensional array of entries
*/
void setData(double[][] data);
/***
* Sets the rounding mode to use when dividing values
* @see java.math.BigDecimal
* @param roundingMode
*/
void setRoundingMode(int roundingMode);
/***
* Gets the rounding mode
* @return
*/
int getRoundingMode();
/**
* Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html">
* maximum absolute row sum norm</a> of the matrix.
*
* @return norm
*/
BigDecimal getNorm();
/**
* Returns the entries in row number <code>row</code> as an array.
*
* @param row the row to be fetched
* @return array of entries in the row
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified row is greater
* than the number of rows in this matrix
*/
BigDecimal[] getRow(int row) throws MatrixIndexException;
/**
* Returns the entries in row number <code>row</code> as an array
* of double values.
*
* @param row the row to be fetched
* @return array of entries in the row
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified row is greater
* than the number of rows in this matrix
*/
double [] getRowAsDoubleArray(int row) throws MatrixIndexException;
/**
* Returns the entries in column number <code>col</code> as an array.
*
* @param col column to fetch
* @return array of entries in the column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified column is greater
* than the number of columns in this matrix
*/
BigDecimal[] getColumn(int col) throws MatrixIndexException;
/**
* Returns the entries in column number <code>col</code> as an array
* of double values.
*
* @param col column to fetch
* @return array of entries in the column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified column is greater
* than the number of columns in this matrix
*/
double [] getColumnAsDoubleArray(int col) throws MatrixIndexException;
/**
* Returns the entry in the specified row and column.
*
* @param row row location of entry to be fetched
* @param column column location of entry to be fetched
* @return matrix entry in row,column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* the dimensions of this matrix
*/
BigDecimal getEntry(int row, int column) throws MatrixIndexException;
/**
* Returns the entry in the specified row and column as a double
*
* @param row row location of entry to be fetched
* @param column column location of entry to be fetched
* @return matrix entry in row,column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* the dimensions of this matrix
*/
double getEntryAsDouble(int row, int column) throws MatrixIndexException;
/**
* Sets the entry in the specified row and column to the specified value.
*
* @param row row location of entry to be set
* @param column column location of entry to be set
* @param value value to set
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* he dimensions of this matrix
*/
void setEntry(int row, int column, BigDecimal value)
throws MatrixIndexException;
/**
* Sets the entry in the specified row and column to the specified value.
*
* @param row row location of entry to be set
* @param column column location of entry to be set
* @param value value to set
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* he dimensions of this matrix
*/
void setEntry(int row, int column, double value)
throws MatrixIndexException;
/**
* Returns the transpose of this matrix.
*
* @return transpose matrix
*/
BigMatrix transpose();
/**
* Returns the inverse of this matrix.
*
* @return inverse matrix
* @throws org.apache.commons.math.linear.InvalidMatrixException if this is not invertible
*/
BigMatrix inverse() throws InvalidMatrixException;
/**
* Returns the determinant of this matrix.
*
* @return determinant
*@throws InvalidMatrixException if matrix is not square
*/
BigDecimal getDeterminant();
/**
* Is this a square matrix?
* @return true if the matrix is square (rowDimension = columnDimension)
*/
boolean isSquare();
/**
* Is this a singular matrix?
* @return true if the matrix is singular
*/
boolean isSingular();
/**
* Returns the number of rows in the matrix.
*
* @return rowDimension
*/
int getRowDimension();
/**
* Returns the number of columns in the matrix.
*
* @return columnDimension
*/
int getColumnDimension();
/**
* Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
* trace</a> of the matrix (the sum of the elements on the main diagonal).
*
* @return trace
*/
BigDecimal getTrace();
/**
* Returns the result of multiplying this by the vector <code>v</code>.
*
* @param v the vector to operate on
* @return this*v
* @throws IllegalArgumentException if columnDimension != v.size()
*/
BigDecimal[] operate(BigDecimal[] v) throws IllegalArgumentException;
/**
* Returns the (row) vector result of premultiplying this by the vector <code>v</code>.
*
* @param v the row vector to premultiply by
* @return v*this
* @throws IllegalArgumentException if rowDimension != v.size()
*/
BigDecimal[] preMultiply(BigDecimal[] v) throws IllegalArgumentException;
/**
* Returns the solution vector for a linear system with coefficient
* matrix = this and constant vector = <code>b</code>.
*
* @param b constant vector
* @return vector of solution values to AX = b, where A is *this
* @throws IllegalArgumentException if this.rowDimension != b.length
* @throws org.apache.commons.math.linear.InvalidMatrixException if this matrix is not square or is singular
*/
BigDecimal[] solve(BigDecimal[] b) throws IllegalArgumentException, InvalidMatrixException;
/**
* Returns a matrix of (column) solution vectors for linear systems with
* coefficient matrix = this and constant vectors = columns of
* <code>b</code>.
*
* @param b matrix of constant vectors forming RHS of linear systems to
* to solve
* @return matrix of solution vectors
* @throws IllegalArgumentException if this.rowDimension != row dimension
* @throws org.apache.commons.math.linear.InvalidMatrixException if this matrix is not square or is singular
*/
BigMatrix solve(BigMatrix b) throws IllegalArgumentException, InvalidMatrixException;
}
1.1 jakarta-commons/math/src/java/org/apache/commons/math/linear/BigMatrixImpl.java
Index: BigMatrixImpl.java
===================================================================
/*
* Copyright 2004 The Apache Software Foundation.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import java.io.Serializable;
import java.math.BigDecimal;
/**
* Implementation for {@link BigMatrix} using a BigDecimal[][] array to store entries
* and <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
* LU decompostion</a> to support linear system
* solution and inverse.
* <p>
* The LU decompostion is performed as needed, to support the following operations: <ul>
* <li>solve</li>
* <li>isSingular</li>
* <li>getDeterminant</li>
* <li>inverse</li> </ul>
* <p>
* <strong>Usage note</strong>:<br>
* The LU decomposition is stored and reused on subsequent calls. If matrix
* data are modified using any of the public setXxx methods, the saved
* decomposition is discarded. If data are modified via references to the
* underlying array obtained using <code>getDataRef()</code>, then the stored
* LU decomposition will not be discarded. In this case, you need to
* explicitly invoke <code>LUDecompose()</code> to recompute the decomposition
* before using any of the methods above.
*
* @version $Revision: 1.1 $ $Date: 2004/06/06 04:20:45 $
*/
public class BigMatrixImpl implements BigMatrix, Serializable {
/** Serialization id */
static final long serialVersionUID = -1011428905656140431L;
private static final BigDecimal ZERO = new BigDecimal(0);
private static final BigDecimal ONE = new BigDecimal(1);
/** Entries of the matrix */
private BigDecimal data[][] = null;
/** Entries of cached LU decomposition.
* All updates to data (other than luDecompose()) *must* set this to null
*/
private BigDecimal lu[][] = null;
/** Permutation associated with LU decompostion */
private int[] permutation = null;
/** Parity of the permutation associated with the LU decomposition */
private int parity = 1;
/** Rounding mode for divisions **/
private int roundingMode = BigDecimal.ROUND_HALF_UP;
/*** BigDecimal scale ***/
private int scale = 64;
/** Bound to determine effective singularity in LU decomposition */
protected static BigDecimal TOO_SMALL = new BigDecimal(10E-12);
/**
* Creates a matrix with no data
*/
public BigMatrixImpl() {
}
/**
* Create a new BigMatrix with the supplied row and column dimensions.
*
* @param rowDimension the number of rows in the new matrix
* @param columnDimension the number of columns in the new matrix
*/
public BigMatrixImpl(int rowDimension, int columnDimension) {
data = new BigDecimal[rowDimension][columnDimension];
lu = null;
}
/**
* Create a new BigMatrix using the <code>data</code> as the underlying
* data array.
* <p>
* The input array is copied, not referenced.
*
* @param d data for new matrix
*/
public BigMatrixImpl(BigDecimal[][] d) {
this.copyIn(d);
lu = null;
}
/**
* Create a new BigMatrix using the <code>data</code> as the underlying
* data array.
* <p>
* The input array is copied, not referenced.
*
* @param d data for new matrix
*/
public BigMatrixImpl(double[][] d) {
this.copyIn(d);
lu = null;
}
/**
* Create a new (column) BigMatrix using <code>v</code> as the
* data for the unique column of the <code>v.length x 1</code> matrix
* created.
* <p>
* The input array is copied, not referenced.
*
* @param v column vector holding data for new matrix
*/
public BigMatrixImpl(BigDecimal[] v) {
int nRows = v.length;
data = new BigDecimal[nRows][1];
for (int row = 0; row < nRows; row++) {
data[row][0] = v[row];
}
}
/**
* Create a new BigMatrix which is a copy of this.
*
* @return the cloned matrix
*/
public BigMatrix copy() {
return new BigMatrixImpl(this.copyOut());
}
/**
* Compute the sum of this and <code>m</code>.
*
* @param m matrix to be added
* @return this + m
* @exception IllegalArgumentException if m is not the same size as this
*/
public BigMatrix add(BigMatrix m) throws IllegalArgumentException {
if (this.getColumnDimension() != m.getColumnDimension() ||
this.getRowDimension() != m.getRowDimension()) {
throw new IllegalArgumentException("matrix dimension mismatch");
}
int rowCount = this.getRowDimension();
int columnCount = this.getColumnDimension();
BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
BigDecimal[][] mData = m.getData();
for (int row = 0; row < rowCount; row++) {
for (int col = 0; col < columnCount; col++) {
outData[row][col] = data[row][col].add(mData[row][col]);
}
}
return new BigMatrixImpl(outData);
}
/**
* Compute this minus <code>m</code>.
*
* @param m matrix to be subtracted
* @return this + m
* @exception IllegalArgumentException if m is not the same size as *this
*/
public BigMatrix subtract(BigMatrix m) throws IllegalArgumentException {
if (this.getColumnDimension() != m.getColumnDimension() ||
this.getRowDimension() != m.getRowDimension()) {
throw new IllegalArgumentException("matrix dimension mismatch");
}
int rowCount = this.getRowDimension();
int columnCount = this.getColumnDimension();
BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
BigDecimal[][] mData = m.getData();
for (int row = 0; row < rowCount; row++) {
for (int col = 0; col < columnCount; col++) {
outData[row][col] = data[row][col].subtract(mData[row][col]);
}
}
return new BigMatrixImpl(outData);
}
/**
* Returns the result of adding d to each entry of this.
*
* @param d value to be added to each entry
* @return d + this
*/
public BigMatrix scalarAdd(BigDecimal d) {
int rowCount = this.getRowDimension();
int columnCount = this.getColumnDimension();
BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
for (int row = 0; row < rowCount; row++) {
for (int col = 0; col < columnCount; col++) {
outData[row][col] = data[row][col].add(d);
}
}
return new BigMatrixImpl(outData);
}
/**
* Returns the result multiplying each entry of this by <code>d</code>
* @param d value to multiply all entries by
* @return d * this
*/
public BigMatrix scalarMultiply(BigDecimal d) {
int rowCount = this.getRowDimension();
int columnCount = this.getColumnDimension();
BigDecimal[][] outData = new BigDecimal[rowCount][columnCount];
for (int row = 0; row < rowCount; row++) {
for (int col = 0; col < columnCount; col++) {
outData[row][col] = data[row][col].multiply(d);
}
}
return new BigMatrixImpl(outData);
}
/**
* Returns the result postmultiplying this by <code>m</code>.
* @param m matrix to postmultiply by
* @return this*m
* @throws IllegalArgumentException
* if columnDimension(this) != rowDimension(m)
*/
public BigMatrix multiply(BigMatrix m) throws IllegalArgumentException {
if (this.getColumnDimension() != m.getRowDimension()) {
throw new IllegalArgumentException("Matrices are not multiplication compatible.");
}
int nRows = this.getRowDimension();
int nCols = m.getColumnDimension();
int nSum = this.getColumnDimension();
BigDecimal[][] mData = m.getData();
BigDecimal[][] outData = new BigDecimal[nRows][nCols];
BigDecimal sum = ZERO;
for (int row = 0; row < nRows; row++) {
for (int col = 0; col < nCols; col++) {
sum = ZERO;
for (int i = 0; i < nSum; i++) {
sum = sum.add(data[row][i].multiply(mData[i][col]));
}
outData[row][col] = sum;
}
}
return new BigMatrixImpl(outData);
}
/**
* Returns the result premultiplying this by <code>m</code>.
* @param m matrix to premultiply by
* @return m * this
* @throws IllegalArgumentException
* if rowDimension(this) != columnDimension(m)
*/
public BigMatrix preMultiply(BigMatrix m) throws IllegalArgumentException {
return m.multiply(this);
}
/**
* Returns matrix entries as a two-dimensional array.
* <p>
* Makes a fresh copy of the underlying data.
*
* @return 2-dimensional array of entries
*/
public BigDecimal[][] getData() {
return copyOut();
}
public double[][] getDataAsDoubleArray() {
int nRows = getRowDimension();
int nCols = getColumnDimension();
double d[][] = new double[nRows][nCols];
for (int i = 0; i < nRows; i++) {
for (int j=0; j<nCols;j++) {
d[i][j] = data[i][j].doubleValue();
}
}
return d;
}
/**
* Overwrites the underlying data for the matrix
* with a fresh copy of <code>inData</code>.
*
* @param inData 2-dimensional array of entries
*/
public void setData(BigDecimal[][] inData) {
copyIn(inData);
lu = null;
}
public void setData(double[][] inData) {
copyIn(inData);
lu = null;
}
/**
* Returns a reference to the underlying data array.
* <p>
* Does not make a fresh copy of the underlying data.
*
* @return 2-dimensional array of entries
*/
public BigDecimal[][] getDataRef() {
return data;
}
/**
* Overwrites the underlying data for the matrix
* with a reference to <code>inData</code>.
* <p>
* Does not make a fresh copy of <code>data</code>.
*
* @param inData 2-dimensional array of entries
*/
public void setDataRef(BigDecimal[][] inData) {
this.data = inData;
lu = null;
}
/***
* Gets the rounding mode for division operations
* The default is {@link BigDecimal.ROUND_HALF_UP}
* @see BigDecimal
* @return
*/
public int getRoundingMode() {
return roundingMode;
}
/***
* Sets the rounding mode for decimal divisions.
* @see BigDecimal
* @param roundingMode
*/
public void setRoundingMode(int roundingMode) {
this.roundingMode = roundingMode;
}
/***
* Sets the scale for division operations.
* The default is 64
* @see BigDecimal
* @return
*/
public int getScale() {
return scale;
}
/***
* Sets the scale for division operations.
* @see BigDecimal
* @param scale
*/
public void setScale(int scale) {
this.scale = scale;
}
/**
* Returns the <a href="http://mathworld.wolfram.com/MaximumAbsoluteRowSumNorm.html">
* maximum absolute row sum norm</a> of the matrix.
*
* @return norm
*/
public BigDecimal getNorm() {
BigDecimal maxColSum = ZERO;
for (int col = 0; col < this.getColumnDimension(); col++) {
BigDecimal sum = ZERO;
for (int row = 0; row < this.getRowDimension(); row++) {
sum = sum.add(data[row][col].abs());
}
maxColSum = maxColSum.max(sum);
}
return maxColSum;
}
/**
* Returns the entries in row number <code>row</code> as an array.
*
* @param row the row to be fetched
* @return array of entries in the row
* @throws MatrixIndexException if the specified row is greater
* than the number of rows in this matrix
*/
public BigDecimal[] getRow(int row) throws MatrixIndexException {
if ( !isValidCoordinate( row, 1 ) ) {
throw new MatrixIndexException("illegal row argument");
}
int ncols = this.getColumnDimension();
BigDecimal[] out = new BigDecimal[ncols];
System.arraycopy(data[row - 1], 0, out, 0, ncols);
return out;
}
/**
* Returns the entries in row number <code>row</code> as an array
* of double values.
*
* @param row the row to be fetched
* @return array of entries in the row
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified row is greater
* than the number of rows in this matrix
*/
public double[] getRowAsDoubleArray(int row) throws MatrixIndexException {
if ( !isValidCoordinate( row, 1 ) ) {
throw new MatrixIndexException("illegal row argument");
}
int ncols = this.getColumnDimension();
double[] out = new double[ncols];
for (int i=0;i<ncols;i++) {
out[i] = data[row-1][i].doubleValue();
}
return out;
}
/**
* Returns the entries in column number <code>col</code> as an array.
*
* @param col column to fetch
* @return array of entries in the column
* @throws MatrixIndexException if the specified column is greater
* than the number of columns in this matrix
*/
public BigDecimal[] getColumn(int col) throws MatrixIndexException {
if ( !isValidCoordinate(1, col) ) {
throw new MatrixIndexException("illegal column argument");
}
int nRows = this.getRowDimension();
BigDecimal[] out = new BigDecimal[nRows];
for (int i = 0; i < nRows; i++) {
out[i] = data[i][col - 1];
}
return out;
}
/**
* Returns the entries in column number <code>col</code> as an array
* of double values.
*
* @param col column to fetch
* @return array of entries in the column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified column is greater
* than the number of columns in this matrix
*/
public double[] getColumnAsDoubleArray(int col) throws MatrixIndexException {
if ( !isValidCoordinate( 1, col ) ) {
throw new MatrixIndexException("illegal column argument");
}
int nrows = this.getRowDimension();
double[] out = new double[nrows];
for (int i=0;i<nrows;i++) {
out[i] = data[i][col-1].doubleValue();
}
return out;
}
/**
* Returns the entry in the specified row and column.
*
* @param row row location of entry to be fetched
* @param column column location of entry to be fetched
* @return matrix entry in row,column
* @throws MatrixIndexException if the specified coordinate is outside
* the dimensions of this matrix
*/
public BigDecimal getEntry(int row, int column)
throws MatrixIndexException {
if (!isValidCoordinate(row,column)) {
throw new MatrixIndexException("matrix entry does not exist");
}
return data[row - 1][column - 1];
}
/**
* Returns the entry in the specified row and column as a double
*
* @param row row location of entry to be fetched
* @param column column location of entry to be fetched
* @return matrix entry in row,column
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* the dimensions of this matrix
*/
public double getEntryAsDouble(int row, int column) throws MatrixIndexException {
return getEntry(row,column).doubleValue();
}
/**
* Sets the entry in the specified row and column to the specified value.
*
* @param row row location of entry to be set
* @param column column location of entry to be set
* @param value value to set
* @throws MatrixIndexException if the specified coordinate is outside
* he dimensions of this matrix
*/
public void setEntry(int row, int column, BigDecimal value)
throws MatrixIndexException {
if (!isValidCoordinate(row,column)) {
throw new MatrixIndexException("matrix entry does not exist");
}
data[row - 1][column - 1] = value;
lu = null;
}
/**
* Sets the entry in the specified row and column to the specified value.
*
* @param row row location of entry to be set
* @param column column location of entry to be set
* @param value value to set
* @throws org.apache.commons.math.linear.MatrixIndexException if the specified coordinate is outside
* he dimensions of this matrix
*/
public void setEntry(int row, int column, double value) throws MatrixIndexException {
setEntry(row, column, new BigDecimal(value));
}
/**
* Returns the transpose matrix.
*
* @return transpose matrix
*/
public BigMatrix transpose() {
int nRows = this.getRowDimension();
int nCols = this.getColumnDimension();
BigMatrixImpl out = new BigMatrixImpl(nCols, nRows);
BigDecimal[][] outData = out.getDataRef();
for (int row = 0; row < nRows; row++) {
for (int col = 0; col < nCols; col++) {
outData[col][row] = data[row][col];
}
}
return out;
}
/**
* Returns the inverse matrix if this matrix is invertible.
*
* @return inverse matrix
* @throws InvalidMatrixException if this is not invertible
*/
public BigMatrix inverse() throws InvalidMatrixException {
return solve(getIdentity(this.getRowDimension()));
}
/**
* Returns the determinant of this matrix.
*
* @return determinant
* @throws InvalidMatrixException if matrix is not square
*/
public BigDecimal getDeterminant() throws InvalidMatrixException {
if (!isSquare()) {
throw new InvalidMatrixException("matrix is not square");
}
if (isSingular()) { // note: this has side effect of attempting LU decomp if lu == null
return ZERO;
} else {
BigDecimal det = (parity == 1) ? ONE : ONE.negate();
for (int i = 0; i < this.getRowDimension(); i++) {
det = det.multiply(lu[i][i]);
}
return det;
}
}
/**
* Is this a square matrix?
* @return true if the matrix is square (rowDimension = columnDimension)
*/
public boolean isSquare() {
return (this.getColumnDimension() == this.getRowDimension());
}
/**
* Is this a singular matrix?
* @return true if the matrix is singular
*/
public boolean isSingular() {
if (lu == null) {
try {
luDecompose();
return false;
} catch (InvalidMatrixException ex) {
return true;
}
} else { // LU decomp must have been successfully performed
return false; // so the matrix is not singular
}
}
/**
* Returns the number of rows in the matrix.
*
* @return rowDimension
*/
public int getRowDimension() {
return data.length;
}
/**
* Returns the number of columns in the matrix.
*
* @return columnDimension
*/
public int getColumnDimension() {
return data[0].length;
}
/**
* Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
* trace</a> of the matrix (the sum of the elements on the main diagonal).
*
* @return trace
*/
public BigDecimal getTrace() throws IllegalArgumentException {
if (!isSquare()) {
throw new IllegalArgumentException("matrix is not square");
}
BigDecimal trace = data[0][0];
for (int i = 1; i < this.getRowDimension(); i++) {
trace = trace.add(data[i][i]);
}
return trace;
}
/**
* Returns the result of multiplying this by the vector <code>v</code>.
*
* @param v the vector to operate on
* @return this*v
* @throws IllegalArgumentException if columnDimension != v.size()
*/
public BigDecimal[] operate(BigDecimal[] v) throws IllegalArgumentException {
if (v.length != this.getColumnDimension()) {
throw new IllegalArgumentException("vector has wrong length");
}
int nRows = this.getRowDimension();
int nCols = this.getColumnDimension();
BigDecimal[] out = new BigDecimal[v.length];
for (int row = 0; row < nRows; row++) {
BigDecimal sum = ZERO;
for (int i = 0; i < nCols; i++) {
sum = sum.add(data[row][i].multiply(v[i]));
}
out[row] = sum;
}
return out;
}
/**
* Returns the result of multiplying this by the vector <code>v</code>.
*
* @param v the vector to operate on
* @return this*v
* @throws IllegalArgumentException if columnDimension != v.size()
*/
public BigDecimal[] operate(double[] v) throws IllegalArgumentException {
BigDecimal bd[] = new BigDecimal[v.length];
for (int i=0;i<bd.length;i++) {
bd[i] = new BigDecimal(v[i]);
}
return operate(bd);
}
/**
* Returns the (row) vector result of premultiplying this by the vector <code>v</code>.
*
* @param v the row vector to premultiply by
* @return v*this
* @throws IllegalArgumentException if rowDimension != v.size()
*/
public BigDecimal[] preMultiply(BigDecimal[] v) throws IllegalArgumentException {
int nRows = this.getRowDimension();
if (v.length != nRows) {
throw new IllegalArgumentException("vector has wrong length");
}
int nCols = this.getColumnDimension();
BigDecimal[] out = new BigDecimal[nCols];
for (int col = 0; col < nCols; col++) {
BigDecimal sum = ZERO;
for (int i = 0; i < nRows; i++) {
sum = sum.add(data[i][col].multiply(v[i]));
}
out[col] = sum;
}
return out;
}
/**
* Returns a matrix of (column) solution vectors for linear systems with
* coefficient matrix = this and constant vectors = columns of
* <code>b</code>.
*
* @param b array of constants forming RHS of linear systems to
* to solve
* @return solution array
* @throws IllegalArgumentException if this.rowDimension != row dimension
* @throws InvalidMatrixException if this matrix is not square or is singular
*/
public BigDecimal[] solve(BigDecimal[] b) throws IllegalArgumentException, InvalidMatrixException {
int nRows = this.getRowDimension();
if (b.length != nRows) {
throw new IllegalArgumentException("constant vector has wrong length");
}
BigMatrix bMatrix = new BigMatrixImpl(b);
BigDecimal[][] solution = ((BigMatrixImpl) (solve(bMatrix))).getDataRef();
BigDecimal[] out = new BigDecimal[nRows];
for (int row = 0; row < nRows; row++) {
out[row] = solution[row][0];
}
return out;
}
/**
* Returns a matrix of (column) solution vectors for linear systems with
* coefficient matrix = this and constant vectors = columns of
* <code>b</code>.
*
* @param b array of constants forming RHS of linear systems to
* to solve
* @return solution array
* @throws IllegalArgumentException if this.rowDimension != row dimension
* @throws InvalidMatrixException if this matrix is not square or is singular
*/
public BigDecimal[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
BigDecimal bd[] = new BigDecimal[b.length];
for (int i=0;i<bd.length;i++) {
bd[i] = new BigDecimal(b[i]);
}
return solve(bd);
}
/**
* Returns a matrix of (column) solution vectors for linear systems with
* coefficient matrix = this and constant vectors = columns of
* <code>b</code>.
*
* @param b matrix of constant vectors forming RHS of linear systems to
* to solve
* @return matrix of solution vectors
* @throws IllegalArgumentException if this.rowDimension != row dimension
* @throws InvalidMatrixException if this matrix is not square or is singular
*/
public BigMatrix solve(BigMatrix b) throws IllegalArgumentException, InvalidMatrixException {
if (b.getRowDimension() != this.getRowDimension()) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (!this.isSquare()) {
throw new InvalidMatrixException("coefficient matrix is not square");
}
if (this.isSingular()) { // side effect: compute LU decomp
throw new InvalidMatrixException("Matrix is singular.");
}
int nCol = this.getColumnDimension();
int nColB = b.getColumnDimension();
int nRowB = b.getRowDimension();
// Apply permutations to b
BigDecimal[][] bv = b.getData();
BigDecimal[][] bp = new BigDecimal[nRowB][nColB];
for (int row = 0; row < nRowB; row++) {
for (int col = 0; col < nColB; col++) {
bp[row][col] = bv[permutation[row]][col];
}
}
bv = null;
// Solve LY = b
for (int col = 0; col < nCol; col++) {
for (int i = col + 1; i < nCol; i++) {
for (int j = 0; j < nColB; j++) {
bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
}
}
}
// Solve UX = Y
for (int col = nCol - 1; col >= 0; col--) {
for (int j = 0; j < nColB; j++) {
bp[col][j] = bp[col][j].divide(lu[col][col], scale, roundingMode);
}
for (int i = 0; i < col; i++) {
for (int j = 0; j < nColB; j++) {
bp[i][j] = bp[i][j].subtract(bp[col][j].multiply(lu[i][col]));
}
}
}
BigMatrixImpl outMat = new BigMatrixImpl(bp);
return outMat;
}
/**
* Computes a new
* <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
* LU decompostion</a> for this matrix, storing the result for use by other methods.
* <p>
* <strong>Implementation Note</strong>:<br>
* Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
* Crout's algortithm</a>, with partial pivoting.
* <p>
* <strong>Usage Note</strong>:<br>
* This method should rarely be invoked directly. Its only use is
* to force recomputation of the LU decomposition when changes have been
* made to the underlying data using direct array references. Changes
* made using setXxx methods will trigger recomputation when needed
* automatically.
*
* @throws InvalidMatrixException if the matrix is non-square or singular.
*/
public void luDecompose() throws InvalidMatrixException {
int nRows = this.getRowDimension();
int nCols = this.getColumnDimension();
if (nRows != nCols) {
throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
}
lu = this.getData();
// Initialize permutation array and parity
permutation = new int[nRows];
for (int row = 0; row < nRows; row++) {
permutation[row] = row;
}
parity = 1;
// Loop over columns
for (int col = 0; col < nCols; col++) {
BigDecimal sum = ZERO;
// upper
for (int row = 0; row < col; row++) {
sum = lu[row][col];
for (int i = 0; i < row; i++) {
sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
}
lu[row][col] = sum;
}
// lower
int max = col; // permutation row
BigDecimal largest = ZERO;
for (int row = col; row < nRows; row++) {
sum = lu[row][col];
for (int i = 0; i < col; i++) {
sum = sum.subtract(lu[row][i].multiply(lu[i][col]));
}
lu[row][col] = sum;
// maintain best permutation choice
if (sum.abs().compareTo(largest) == 1) {
largest = sum.abs();
max = row;
}
}
// Singularity check
if (lu[max][col].abs().compareTo(TOO_SMALL) <= 0) {
lu = null;
throw new InvalidMatrixException("matrix is singular");
}
// Pivot if necessary
if (max != col) {
BigDecimal tmp = ZERO;
for (int i = 0; i < nCols; i++) {
tmp = lu[max][i];
lu[max][i] = lu[col][i];
lu[col][i] = tmp;
}
int temp = permutation[max];
permutation[max] = permutation[col];
permutation[col] = temp;
parity = -parity;
}
//Divide the lower elements by the "winning" diagonal elt.
for (int row = col + 1; row < nRows; row++) {
lu[row][col] = lu[row][col].divide(lu[col][col], scale, roundingMode);
}
}
}
/**
*
* @see Object#toString()
*/
public String toString() {
StringBuffer res = new StringBuffer();
res.append("BigMatrixImpl{");
for (int i = 0; i < data.length; i++) {
if (i > 0)
res.append(",");
res.append("{");
for (int j = 0; j < data[0].length; j++) {
if (j > 0)
res.append(",");
res.append(data[i][j]);
} //for
res.append("}");
} //for
res.append("}");
return res.toString();
} //toString
//------------------------ Protected methods
/**
* Returns <code>dimension x dimension</code> identity matrix.
*
* @param dimension dimension of identity matrix to generate
* @return identity matrix
*/
protected BigMatrix getIdentity(int dimension) {
BigMatrixImpl out = new BigMatrixImpl(dimension, dimension);
BigDecimal[][] d = out.getDataRef();
for (int row = 0; row < dimension; row++) {
for (int col = 0; col < dimension; col++) {
d[row][col] = row == col ? ONE : ZERO;
}
}
return out;
}
/**
* Returns the LU decomposition as a BigMatrix.
* Returns a fresh copy of the cached LU matrix if this has been computed;
* otherwise the composition is computed and cached for use by other methods.
* Since a copy is returned in either case, changes to the returned matrix do not
* affect the LU decomposition property.
* <p>
* The matrix returned is a compact representation of the LU decomposition.
* Elements below the main diagonal correspond to entries of the "L" matrix;
* elements on and above the main diagonal correspond to entries of the "U"
* matrix.
* <p>
* Example: <pre>
*
* Returned matrix L U
* 2 3 1 1 0 0 2 3 1
* 5 4 6 5 1 0 0 4 6
* 1 7 8 1 7 1 0 0 8
* </pre>
*
* The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
* where permuteRows reorders the rows of the matrix to follow the order determined
* by the <a href=#getPermutation()>permutation</a> property.
*
* @return LU decomposition matrix
* @throws InvalidMatrixException if the matrix is non-square or singular.
*/
protected BigMatrix getLUMatrix() throws InvalidMatrixException {
if (lu == null) {
luDecompose();
}
return new BigMatrixImpl(lu);
}
/**
* Returns the permutation associated with the lu decomposition.
* The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
* <p>
* Example:
* permutation = [1, 2, 0] means current 2nd row is first, current third row is second
* and current first row is last.
* <p>
* Returns a fresh copy of the array.
*
* @return the permutation
*/
protected int[] getPermutation() {
int[] out = new int[permutation.length];
System.arraycopy(permutation, 0, out, 0, permutation.length);
return out;
}
//------------------------ Private methods
/**
* Returns a fresh copy of the underlying data array.
*
* @return a copy of the underlying data array.
*/
private BigDecimal[][] copyOut() {
int nRows = this.getRowDimension();
BigDecimal[][] out = new BigDecimal[nRows][this.getColumnDimension()];
// can't copy 2-d array in one shot, otherwise get row references
for (int i = 0; i < nRows; i++) {
System.arraycopy(data[i], 0, out[i], 0, data[i].length);
}
return out;
}
/**
* Replaces data with a fresh copy of the input array.
*
* @param in data to copy in
*/
private void copyIn(BigDecimal[][] in) {
int nRows = in.length;
int nCols = in[0].length;
data = new BigDecimal[nRows][nCols];
System.arraycopy(in, 0, data, 0, in.length);
for (int i = 0; i < nRows; i++) {
System.arraycopy(in[i], 0, data[i], 0, nCols);
}
lu = null;
}
/**
* Replaces data with a fresh copy of the input array.
*
* @param in data to copy in
*/
private void copyIn(double[][] in) {
int nRows = in.length;
int nCols = in[0].length;
data = new BigDecimal[nRows][nCols];
for (int i = 0; i < nRows; i++) {
for (int j=0; j < nCols; j++) {
data[i][j] = new BigDecimal(in[i][j]);
}
}
lu = null;
}
/**
* Tests a given coordinate as being valid or invalid
*
* @param row the row index.
* @param col the column index.
* @return true if the coordinate is with the current dimensions
*/
private boolean isValidCoordinate(int row, int col) {
int nRows = this.getRowDimension();
int nCols = this.getColumnDimension();
return !(row < 1 || row > nRows || col < 1 || col > nCols);
}
}
1.1 jakarta-commons/math/src/test/org/apache/commons/math/linear/BigMatrixImplTest.java
Index: BigMatrixImplTest.java
===================================================================
/*
* Copyright 2004 The Apache Software Foundation.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
import java.math.BigDecimal;
/**
* Test cases for the {@link BigMatrixImpl} class.
*
* @version $Revision: 1.1 $ $Date: 2004/06/06 04:20:45 $
*/
public final class BigMatrixImplTest extends TestCase {
private double[][] testData = { {1d,2d,3d}, {2d,5d,3d}, {1d,0d,8d} };
private double[][] testDataLU = {{2d, 5d, 3d}, {.5d, -2.5d, 6.5d}, {0.5d, 0.2d, .2d}};
private double[][] testDataPlus2 = { {3d,4d,5d}, {4d,7d,5d}, {3d,2d,10d} };
private double[][] testDataMinus = { {-1d,-2d,-3d}, {-2d,-5d,-3d},
{-1d,0d,-8d} };
private double[] testDataRow1 = {1d,2d,3d};
private double[] testDataCol3 = {3d,3d,8d};
private double[][] testDataInv =
{ {-40d,16d,9d}, {13d,-5d,-3d}, {5d,-2d,-1d} };
private double[] preMultTest = {8,12,33};
private double[][] testData2 ={ {1d,2d,3d}, {2d,5d,3d}};
private double[][] testData2T = { {1d,2d}, {2d,5d}, {3d,3d}};
private double[][] testDataPlusInv =
{ {-39d,18d,12d}, {15d,0d,0d}, {6d,-2d,7d} };
private double[][] id = { {1d,0d,0d}, {0d,1d,0d}, {0d,0d,1d} };
private double[][] luData = { {2d,3d,3d}, {0d,5d,7d}, {6d,9d,8d} };
private double[][] luDataLUDecomposition = { {6d,9d,8d}, {0d,5d,7d}, {0.33333333333333,0d,0.33333333333333} };
private double[][] singular = { {2d,3d}, {2d,3d} };
private double[][] bigSingular = {{1d,2d,3d,4d}, {2d,5d,3d,4d},
{7d,3d,256d,1930d}, {3d,7d,6d,8d}}; // 4th row = 1st + 2nd
private double[][] detData = { {1d,2d,3d}, {4d,5d,6d}, {7d,8d,10d} };
private double[][] detData2 = { {1d, 3d}, {2d, 4d}};
private double[] testVector = {1,2,3};
private double[] testVector2 = {1,2,3,4};
private double entryTolerance = 10E-16;
private double normTolerance = 10E-14;
public BigMatrixImplTest(String name) {
super(name);
}
public void setUp() {
}
public static Test suite() {
TestSuite suite = new TestSuite(BigMatrixImplTest.class);
suite.setName("BigMatrixImpl Tests");
return suite;
}
public static final double[] asDouble(BigDecimal[] data) {
double d[] = new double[data.length];
for (int i=0;i<d.length;i++) {
d[i] = data[i].doubleValue();
}
return d;
}
public static final double[][] asDouble(BigDecimal[][] data) {
double d[][] = new double[data.length][data[0].length];
for (int i=0;i<d.length;i++) {
for (int j=0;j<d[i].length;j++)
d[i][j] = data[i][j].doubleValue();
}
return d;
}
public static final BigDecimal[] asBigDecimal(double [] data) {
BigDecimal d[] = new BigDecimal[data.length];
for (int i=0;i<d.length;i++) {
d[i] = new BigDecimal(data[i]);
}
return d;
}
public static final BigDecimal[][] asBigDecimal(double [][] data) {
BigDecimal d[][] = new BigDecimal[data.length][data[0].length];
for (int i=0;i<d.length;i++) {
for (int j=0;j<data[i].length;j++) {
d[i][j] = new BigDecimal(data[i][j]);
}
}
return d;
}
/** test dimensions */
public void testDimensions() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
assertEquals("testData row dimension",3,m.getRowDimension());
assertEquals("testData column dimension",3,m.getColumnDimension());
assertTrue("testData is square",m.isSquare());
assertEquals("testData2 row dimension",m2.getRowDimension(),2);
assertEquals("testData2 column dimension",m2.getColumnDimension(),3);
assertTrue("testData2 is not square",!m2.isSquare());
BigMatrixImpl m3 = new BigMatrixImpl();
m3.setData(testData);
}
/** test copy functions */
public void testCopyFunctions() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
m2.setData(m.getData());
assertClose("getData",m2,m,entryTolerance);
// no dangling reference...
m2.setEntry(1,1,2000d);
BigMatrixImpl m3 = new BigMatrixImpl(testData);
assertClose("no getData side effect",m,m3,entryTolerance);
m3 = (BigMatrixImpl) m.copy();
double[][] stompMe = {{1d,2d,3d}};
m3.setDataRef(asBigDecimal(stompMe));
assertClose("no copy side effect",m,new BigMatrixImpl(testData),
entryTolerance);
}
/** test add */
public void testAdd() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl mInv = new BigMatrixImpl(testDataInv);
BigMatrixImpl mPlusMInv = (BigMatrixImpl)m.add(mInv);
double[][] sumEntries = asDouble(mPlusMInv.getData());
for (int row = 0; row < m.getRowDimension(); row++) {
for (int col = 0; col < m.getColumnDimension(); col++) {
assertEquals("sum entry entry",
testDataPlusInv[row][col],sumEntries[row][col],
entryTolerance);
}
}
}
/** test add failure */
public void testAddFail() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
try {
BigMatrixImpl mPlusMInv = (BigMatrixImpl)m.add(m2);
fail("IllegalArgumentException expected");
} catch (IllegalArgumentException ex) {
;
}
}
/** test norm */
public void testNorm() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
assertEquals("testData norm",14d,m.getNorm().doubleValue(),entryTolerance);
assertEquals("testData2 norm",7d,m2.getNorm().doubleValue(),entryTolerance);
}
/** test m-n = m + -n */
public void testPlusMinus() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl m2 = new BigMatrixImpl(testDataInv);
assertClose("m-n = m + -n",m.subtract(m2),
m2.scalarMultiply(new BigDecimal(-1d)).add(m),entryTolerance);
try {
BigMatrix a = m.subtract(new BigMatrixImpl(testData2));
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
/** test multiply */
public void testMultiply() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl mInv = new BigMatrixImpl(testDataInv);
BigMatrixImpl identity = new BigMatrixImpl(id);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
assertClose("inverse multiply",m.multiply(mInv),
identity,entryTolerance);
assertClose("inverse multiply",mInv.multiply(m),
identity,entryTolerance);
assertClose("identity multiply",m.multiply(identity),
m,entryTolerance);
assertClose("identity multiply",identity.multiply(mInv),
mInv,entryTolerance);
assertClose("identity multiply",m2.multiply(identity),
m2,entryTolerance);
try {
BigMatrix a = m.multiply(new BigMatrixImpl(bigSingular));
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
//Additional Test for BigMatrixImplTest.testMultiply
private double[][] d3 = new double[][] {{1,2,3,4},{5,6,7,8}};
private double[][] d4 = new double[][] {{1},{2},{3},{4}};
private double[][] d5 = new double[][] {{30},{70}};
public void testMultiply2() {
BigMatrix m3 = new BigMatrixImpl(d3);
BigMatrix m4 = new BigMatrixImpl(d4);
BigMatrix m5 = new BigMatrixImpl(d5);
assertClose("m3*m4=m5", m3.multiply(m4), m5, entryTolerance);
}
/** test isSingular */
public void testIsSingular() {
BigMatrixImpl m = new BigMatrixImpl(singular);
assertTrue("singular",m.isSingular());
m = new BigMatrixImpl(bigSingular);
assertTrue("big singular",m.isSingular());
m = new BigMatrixImpl(id);
assertTrue("identity nonsingular",!m.isSingular());
m = new BigMatrixImpl(testData);
assertTrue("testData nonsingular",!m.isSingular());
}
/** test inverse */
public void testInverse() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrix mInv = new BigMatrixImpl(testDataInv);
assertClose("inverse",mInv,m.inverse(),normTolerance);
assertClose("inverse^2",m,m.inverse().inverse(),10E-12);
// Not square
m = new BigMatrixImpl(testData2);
try {
m.inverse();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
// expected
}
// Singular
m = new BigMatrixImpl(singular);
try {
m.inverse();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
// expected
}
}
/** test solve */
public void testSolve() {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrix mInv = new BigMatrixImpl(testDataInv);
// being a bit slothful here -- actually testing that X = A^-1 * B
assertClose("inverse-operate",
asDouble(mInv.operate(asBigDecimal(testVector))),
asDouble(m.solve(asBigDecimal(testVector))),
normTolerance);
try {
double[] x = asDouble(m.solve(asBigDecimal(testVector2)));
fail("expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
BigMatrix bs = new BigMatrixImpl(bigSingular);
try {
BigMatrix a = bs.solve(bs);
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
;
}
try {
BigMatrix a = m.solve(bs);
fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
try {
BigMatrix a = (new BigMatrixImpl(testData2)).solve(bs);
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
try {
(new BigMatrixImpl(testData2)).luDecompose();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
;
}
}
/** test determinant */
public void testDeterminant() {
BigMatrix m = new BigMatrixImpl(bigSingular);
assertEquals("singular determinant",0,m.getDeterminant().doubleValue(),0);
m = new BigMatrixImpl(detData);
assertEquals("nonsingular test",-3d,m.getDeterminant().doubleValue(),normTolerance);
// Examples verified against R (version 1.8.1, Red Hat Linux 9)
m = new BigMatrixImpl(detData2);
assertEquals("nonsingular R test 1",-2d,m.getDeterminant().doubleValue(),normTolerance);
m = new BigMatrixImpl(testData);
assertEquals("nonsingular R test 2",-1d,m.getDeterminant().doubleValue(),normTolerance);
try {
double a = new BigMatrixImpl(testData2).getDeterminant().doubleValue();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
;
}
}
/** test trace */
public void testTrace() {
BigMatrix m = new BigMatrixImpl(id);
assertEquals("identity trace",3d,m.getTrace().doubleValue(),entryTolerance);
m = new BigMatrixImpl(testData2);
try {
double x = m.getTrace().doubleValue();
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
/** test sclarAdd */
public void testScalarAdd() {
BigMatrix m = new BigMatrixImpl(testData);
assertClose("scalar add",new BigMatrixImpl(testDataPlus2),
m.scalarAdd(new BigDecimal(2d)),entryTolerance);
}
/** test operate */
public void testOperate() {
BigMatrix m = new BigMatrixImpl(id);
double[] x = asDouble(m.operate(asBigDecimal(testVector)));
assertClose("identity operate",testVector,x,entryTolerance);
m = new BigMatrixImpl(bigSingular);
try {
x = asDouble(m.operate(asBigDecimal(testVector)));
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
/** test transpose */
public void testTranspose() {
BigMatrix m = new BigMatrixImpl(testData);
assertClose("inverse-transpose",m.inverse().transpose(),
m.transpose().inverse(),normTolerance);
m = new BigMatrixImpl(testData2);
BigMatrix mt = new BigMatrixImpl(testData2T);
assertClose("transpose",mt,m.transpose(),normTolerance);
}
/** test preMultiply by vector */
public void testPremultiplyVector() {
BigMatrix m = new BigMatrixImpl(testData);
assertClose("premultiply",asDouble(m.preMultiply(asBigDecimal(testVector))),preMultTest,normTolerance);
m = new BigMatrixImpl(bigSingular);
try {
m.preMultiply(asBigDecimal(testVector));
fail("expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
public void testPremultiply() {
BigMatrix m3 = new BigMatrixImpl(d3);
BigMatrix m4 = new BigMatrixImpl(d4);
BigMatrix m5 = new BigMatrixImpl(d5);
assertClose("m3*m4=m5", m4.preMultiply(m3), m5, entryTolerance);
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrixImpl mInv = new BigMatrixImpl(testDataInv);
BigMatrixImpl identity = new BigMatrixImpl(id);
BigMatrixImpl m2 = new BigMatrixImpl(testData2);
assertClose("inverse multiply",m.preMultiply(mInv),
identity,entryTolerance);
assertClose("inverse multiply",mInv.preMultiply(m),
identity,entryTolerance);
assertClose("identity multiply",m.preMultiply(identity),
m,entryTolerance);
assertClose("identity multiply",identity.preMultiply(mInv),
mInv,entryTolerance);
try {
BigMatrix a = m.preMultiply(new BigMatrixImpl(bigSingular));
fail("Expecting illegalArgumentException");
} catch (IllegalArgumentException ex) {
;
}
}
public void testGetVectors() {
BigMatrix m = new BigMatrixImpl(testData);
assertClose("get row",m.getRowAsDoubleArray(1),testDataRow1,entryTolerance);
assertClose("get col",m.getColumnAsDoubleArray(3),testDataCol3,entryTolerance);
try {
double[] x = m.getRowAsDoubleArray(10);
fail("expecting MatrixIndexException");
} catch (MatrixIndexException ex) {
;
}
try {
double[] x = m.getColumnAsDoubleArray(-1);
fail("expecting MatrixIndexException");
} catch (MatrixIndexException ex) {
;
}
}
public void testEntryMutators() {
BigMatrix m = new BigMatrixImpl(testData);
assertEquals("get entry",m.getEntry(1,2).doubleValue(),2d,entryTolerance);
m.setEntry(1,2,100d);
assertEquals("get entry",m.getEntry(1,2).doubleValue(),100d,entryTolerance);
try {
double x = m.getEntry(0,2).doubleValue();
fail("expecting MatrixIndexException");
} catch (MatrixIndexException ex) {
;
}
try {
m.setEntry(1,4,200d);
fail("expecting MatrixIndexException");
} catch (MatrixIndexException ex) {
;
}
}
public void testLUDecomposition() throws Exception {
BigMatrixImpl m = new BigMatrixImpl(testData);
BigMatrix lu = m.getLUMatrix();
assertClose("LU decomposition", lu, (BigMatrix) new BigMatrixImpl(testDataLU), normTolerance);
verifyDecomposition(m, lu);
m = new BigMatrixImpl(luData);
lu = m.getLUMatrix();
assertClose("LU decomposition", lu, (BigMatrix) new BigMatrixImpl(luDataLUDecomposition), normTolerance);
verifyDecomposition(m, lu);
m = new BigMatrixImpl(testDataMinus);
lu = m.getLUMatrix();
verifyDecomposition(m, lu);
m = new BigMatrixImpl(id);
lu = m.getLUMatrix();
verifyDecomposition(m, lu);
try {
m = new BigMatrixImpl(bigSingular); // singular
lu = m.getLUMatrix();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
// expected
}
try {
m = new BigMatrixImpl(testData2); // not square
lu = m.getLUMatrix();
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ex) {
// expected
}
}
//--------------- -----------------Protected methods
/** verifies that two matrices are close (1-norm) */
protected void assertClose(String msg, BigMatrix m, BigMatrix n,
double tolerance) {
assertTrue(msg,m.subtract(n).getNorm().doubleValue() < tolerance);
}
/** verifies that two vectors are close (sup norm) */
protected void assertClose(String msg, double[] m, double[] n,
double tolerance) {
if (m.length != n.length) {
fail("vectors not same length");
}
for (int i = 0; i < m.length; i++) {
assertEquals(msg + " " + i + " elements differ",
m[i],n[i],tolerance);
}
}
/** extracts the l and u matrices from compact lu representation */
protected void splitLU(BigMatrix lu, BigMatrix lower, BigMatrix upper) throws InvalidMatrixException {
if (!lu.isSquare() || !lower.isSquare() || !upper.isSquare() ||
lower.getRowDimension() != upper.getRowDimension()
|| lower.getRowDimension() != lu.getRowDimension()) {
throw new InvalidMatrixException("incorrect dimensions");
}
int n = lu.getRowDimension();
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
if (j < i) {
lower.setEntry(i, j, lu.getEntry(i, j));
upper.setEntry(i, j, 0d);
} else if (i == j) {
lower.setEntry(i, j, 1d);
upper.setEntry(i, j, lu.getEntry(i, j));
} else {
lower.setEntry(i, j, 0d);
upper.setEntry(i, j, lu.getEntry(i, j));
}
}
}
}
/** Returns the result of applying the given row permutation to the matrix */
protected BigMatrix permuteRows(BigMatrix matrix, int[] permutation) {
if (!matrix.isSquare() || matrix.getRowDimension() != permutation.length) {
throw new IllegalArgumentException("dimension mismatch");
}
int n = matrix.getRowDimension();
BigMatrix out = new BigMatrixImpl(n, n);
for (int i =1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
out.setEntry(i, j, matrix.getEntry(permutation[i -1] + 1, j));
}
}
return out;
}
/** Extracts l and u matrices from lu and verifies that matrix = l times u modulo permutation */
protected void verifyDecomposition(BigMatrix matrix, BigMatrix lu) throws Exception{
int n = matrix.getRowDimension();
BigMatrix lower = new BigMatrixImpl(n, n);
BigMatrix upper = new BigMatrixImpl(n, n);
splitLU(lu, lower, upper);
int[] permutation = ((BigMatrixImpl) matrix).getPermutation();
BigMatrix permuted = permuteRows(matrix, permutation);
assertClose("lu decomposition does not work", permuted, lower.multiply(upper), normTolerance);
}
/** Useful for debugging */
private void dumpMatrix(BigMatrix m) {
for (int i = 0; i < m.getRowDimension(); i++) {
String os = "";
for (int j = 0; j < m.getColumnDimension(); j++) {
os += m.getEntry(i+1, j+1) + " ";
}
System.out.println(os);
}
}
}
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org