|
import java.text.NumberFormat;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.Locale;
import java.text.FieldPosition;
import java.io.PrintWriter;
import java.io.BufferedReader;
import java.io.StreamTokenizer;
/**
Jama = Java Matrix class.
The Java Matrix Class provides the fundamental operations of numerical
linear algebra. Various constructors create Matrices from two dimensional
arrays of double precision floating point numbers. Various "gets" and
"sets" provide access to submatrices and matrix elements. Several methods
implement basic matrix arithmetic, including matrix addition and
multiplication, matrix norms, and element-by-element array operations.
Methods for reading and printing matrices are also included. All the
operations in this version of the Matrix Class involve real matrices.
Complex matrices may be handled in a future version.
Five fundamental matrix decompositions, which consist of pairs or triples
of matrices, permutation vectors, and the like, produce results in five
decomposition classes. These decompositions are accessed by the Matrix
class to compute solutions of simultaneous linear equations, determinants,
inverses and other matrix functions. The five decompositions are:
Cholesky Decomposition of symmetric, positive definite matrices.
LU Decomposition of rectangular matrices.
QR Decomposition of rectangular matrices.
Singular Value Decomposition of rectangular matrices.
Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
Example of use:
Solve a linear system A x = b and compute the residual norm, ||b - A x||.
double[][] vals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
Matrix A = new Matrix(vals);
Matrix b = Matrix.random(3,1);
Matrix x = A.solve(b);
Matrix r = A.times(x).minus(b);
double rnorm = r.normInf();
@author The MathWorks, Inc. and the National Institute of Standards and
Technology.
@version 5 August 1998
*/
public class MatrixJama
//implements Cloneable, java.io.Serializable {
{
/* ------------------------
Class variables
* ------------------------ */
/** Array for internal storage of elements.
@serial internal array storage.
*/
private double[][] A;
/** Row and column dimensions.
@serial row dimension.
@serial column dimension.
*/
private int m, n;
/* ------------------------
Constructors
* ------------------------ */
/** Construct an m-by-n matrix of zeros.
@param m Number of rows.
@param n Number of colums.
*/
public MatrixJama (int m, int n)
{
this.m = m;
this.n = n;
A = new double[m][n];
}
/** Construct a matrix from a 2-D array.
@param A Two-dimensional array of doubles.
@exception IllegalArgumentException All rows must have the same
length
@see #constructWithCopy
*/
public MatrixJama (double[][] A)
{
m = A.length;
n = A[0].length;
for (int i = 0; i < m; i++)
{
if (A[i].length != n)
{
throw new
IllegalArgumentException("All rows must have the same length.");
}
}
this.A = A;
}
/** Construct a matrix quickly without
checking arguments.
@param A Two-dimensional array of doubles.
@param m Number of rows.
@param n Number of colums.
*/
public MatrixJama(double[][] A, int m, int n)
{
this.A = A;
this.m = m;
this.n = n;
}
/* ------------------------
Public Methods
* ------------------------ */
public double[][] getArray()
{
return A;
}
/** Copy the internal two-dimensional array.
@return Two-dimensional array copy of matrix elements.
*/
public double[][] getArrayCopy()
{
double[][] C = new double[m][n];
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[i][j] = A[i][j];
}
}
return C;
}
/** Get row dimension.
@return m, the number of rows.
*/
public int getRowDimension()
{
return m;
}
/** Get column dimension.
@return n, the number of columns.
*/
public int getColumnDimension()
{
return n;
}
/** Get a single element.
@param i Row index.
@param j Column index.
@return A(i,j)
@exception ArrayIndexOutOfBoundsException
*/
public double get(int i, int j)
{
return A[i][j];
}
/** Get a submatrix.
@param i0 Initial row index
@param i1 Final row index
@param j0 Initial column index
@param j1 Final column index
@return A(i0:i1,j0:j1)
@exception ArrayIndexOutOfBoundsException Submatrix indices
*/
public MatrixJama getMatrix (int i0, int i1, int j0, int j1)
{
MatrixJama X = new MatrixJama(i1-i0+1,j1-j0+1);
double[][] B = X.getArray();
try
{
for (int i = i0; i <= i1; i++)
{
for (int j = j0; j <= j1; j++)
{
B[i-i0][j-j0] =
A[i][j];
}
}
}
catch(ArrayIndexOutOfBoundsException e)
{
throw new
ArrayIndexOutOfBoundsException("Submatrix indices");
}
return X;
}
/** Set a single element.
@param i Row index.
@param j Column index.
@param s A(i,j).
@exception ArrayIndexOutOfBoundsException
*/
public void set (int i, int j, double s)
{
A[i][j] = s;
}
/** Set a submatrix.
@param i0 Initial row index
@param i1 Final row index
@param j0 Initial column index
@param j1 Final column index
@param X A(i0:i1,j0:j1)
@exception ArrayIndexOutOfBoundsException Submatrix indices
*/
public void setMatrix(int i0, int i1, int j0, int j1, MatrixJama X)
{
try
{
for (int i = i0; i <= i1; i++)
{
for (int j = j0; j <= j1; j++)
{
A[i][j] =
X.get(i-i0,j-j0);
}
}
}
catch(ArrayIndexOutOfBoundsException e)
{
throw new
ArrayIndexOutOfBoundsException("Submatrix indices");
}
}
/** Matrix transpose.
@return A'
*/
public MatrixJama transpose ()
{
MatrixJama X = new MatrixJama(n,m);
double[][] C = X.getArray();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[j][i] = A[i][j];
}
}
return X;
}
/** C = A + B
@param B another matrix
@return A + B
*/
public MatrixJama plus(MatrixJama B)
{
checkMatrixDimensions(B);
MatrixJama X = new MatrixJama(m,n);
double[][] C = X.getArray();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[i][j] = A[i][j] + B.A[i][j];
}
}
return X;
}
/** C = A - B
@param B another matrix
@return A - B
*/
public MatrixJama minus(MatrixJama B)
{
checkMatrixDimensions(B);
MatrixJama X = new MatrixJama(m,n);
double[][] C = X.getArray();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[i][j] = A[i][j] - B.A[i][j];
}
}
return X;
}
/** Multiply a matrix by a scalar, C = s*A
@param s scalar
@return s*A
*/
public MatrixJama times(double s)
{
MatrixJama X = new MatrixJama(m,n);
double[][] C = X.getArray();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[i][j] = s*A[i][j];
}
}
return X;
}
/** Linear algebraic matrix multiplication, A
* B
@param B another matrix
@return Matrix product, A * B
@exception IllegalArgumentException Matrix inner dimensions must
agree.
*/
public MatrixJama times(MatrixJama B)
{
if (B.m != n)
{
throw new IllegalArgumentException("Matrix
inner dimensions must agree.");
}
MatrixJama X = new MatrixJama(m,B.n);
double[][] C = X.getArray();
double[] Bcolj = new double[n];
for (int j = 0; j < B.n; j++)
{
for (int k = 0; k < n; k++)
{
Bcolj[k] = B.A[k][j];
}
for (int i = 0; i < m; i++)
{
double[] Arowi = A[i];
double s = 0;
for (int k = 0; k < n; k++)
{
s += Arowi[k]*Bcolj[k];
}
C[i][j] = s;
}
}
return X;
}
public EigenvalueDecomposition eig()
{
return new EigenvalueDecomposition(this);
}
/** Generate identity matrix
@param m Number of rows.
@param n Number of colums.
@return An m-by-n matrix with ones on the diagonal and zeros
elsewhere.
*/
public static MatrixJama identity(int m, int n)
{
MatrixJama A = new MatrixJama(m,n);
double[][] X = A.getArray();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
X[i][j] = (i == j ? 1.0 : 0.0);
}
}
return A;
}
/** Print the matrix to stdout. Line the
elements up in columns
* with a Fortran-like 'Fw.d' style format.
@param w Column width.
@param d Number of digits after the decimal.
*/
public void print(int w, int d)
{
print(new PrintWriter(System.out,true),w,d);
}
/** Print the matrix to the output stream.
Line the elements up in
* columns with a Fortran-like 'Fw.d' style format.
@param output Output stream.
@param w Column width.
@param d Number of digits after the decimal.
*/
public void print(PrintWriter output, int w, int d)
{
DecimalFormat format = new DecimalFormat();
format.setDecimalFormatSymbols(new
DecimalFormatSymbols(Locale.US));
format.setMinimumIntegerDigits(1);
format.setMaximumFractionDigits(d);
format.setMinimumFractionDigits(d);
format.setGroupingUsed(false);
print(output,format,w+2);
}
/** Print the matrix to stdout. Line the
elements up in columns.
* Use the format object, and right justify within columns of width
* characters.
* Note that is the matrix is to be read back in, you probably will
want
* to use a NumberFormat that is set to US Locale.
@param format A Formatting object for individual elements.
@param width Field width for each column.
@see java.text.DecimalFormat#setDecimalFormatSymbols
*/
public void print(NumberFormat format, int width)
{
print(new PrintWriter(System.out,true),format,width);
}
// DecimalFormat is a little disappointing
coming from Fortran or C's printf.
// Since it doesn't pad on the left, the elements will come out
different
// widths. Consequently, we'll pass the desired column width in as
an
// argument and do the extra padding ourselves.
/** Print the matrix to the output stream.
Line the elements up in columns.
* Use the format object, and right justify within columns of width
* characters.
* Note that is the matrix is to be read back in, you probably will
want
* to use a NumberFormat that is set to US Locale.
@param output the output stream.
@param format A formatting object to format the matrix elements
@param width Column width.
@see java.text.DecimalFormat#setDecimalFormatSymbols
*/
public void print(PrintWriter output, NumberFormat format, int
width)
{
output.println(); // start on new line.
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
String s = format.format(A[i][j]);
// format the number
int padding =
Math.max(1,width-s.length()); // At _least_ 1 space
for (int k = 0; k < padding;
k++)
output.print(' ');
output.print(s);
}
output.println();
}
output.println(); // end with blank line.
}
/** Read a matrix from a stream. The format
is the same the print method,
* so printed matrices can be read back in (provided they were
printed using
* US Locale). Elements are separated by
* whitespace, all the elements for each row appear on a single
line,
* the last row is followed by a blank line.
@param input the input stream.
*/
public static MatrixJama read(BufferedReader input) throws
java.io.IOException
{
StreamTokenizer tokenizer= new StreamTokenizer(input);
// Although StreamTokenizer will parse
numbers, it doesn't recognize
// scientific notation (E or D); however, Double.valueOf does.
// The strategy here is to disable StreamTokenizer's number parsing.
// We'll only get whitespace delimited words, EOL's and EOF's.
// These words should all be numbers, for Double.valueOf to parse.
tokenizer.resetSyntax();
tokenizer.wordChars(0,255);
tokenizer.whitespaceChars(0, ' ');
tokenizer.eolIsSignificant(true);
java.util.Vector v = new java.util.Vector();
// Ignore initial empty lines
while (tokenizer.nextToken() == StreamTokenizer.TT_EOL);
if (tokenizer.ttype == StreamTokenizer.TT_EOF)
throw new java.io.IOException("Unexpected
EOF on matrix read.");
do
{
v.addElement(Double.valueOf(tokenizer.sval));
// Read & store 1st row.
}
while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
int n = v.size(); // Now we've got the
number of columns!
double row[] = new double[n];
for (int j=0; j<n; j++) // extract the elements of the
1st row.
row[j]=((Double)v.elementAt(j)).doubleValue();
v.removeAllElements();
v.addElement(row); // Start storing rows instead of
columns.
while (tokenizer.nextToken() == StreamTokenizer.TT_WORD)
{
// While non-empty lines
v.addElement(row = new double[n]);
int j = 0;
do
{
if (j >= n) throw new
java.io.IOException("Row " + v.size() + " is too long.");
row[j++] =
Double.valueOf(tokenizer.sval).doubleValue();
}
while (tokenizer.nextToken() ==
StreamTokenizer.TT_WORD);
if (j < n)
throw new
java.io.IOException("Row " + v.size() + " is too short.");
}
int m = v.size(); // Now we've got the number of rows.
double[][] A = new double[m][];
v.copyInto(A); // copy the rows out of the vector
return new MatrixJama(A);
}
/* ------------------------
Private Methods
* ------------------------ */
/** Check if size(A) == size(B) **/
private void checkMatrixDimensions(MatrixJama B)
{
if (B.m != m || B.n != n)
{
throw new IllegalArgumentException("Matrix
dimensions must agree.");
}
}
}