edu.mines.jtk.lapack.DMatrix Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2005, Colorado School of Mines and others.
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 edu.mines.jtk.lapack;
import org.netlib.blas.BLAS;
import static edu.mines.jtk.util.ArrayMath.*;
import edu.mines.jtk.util.Check;
/**
* A double-precision matrix.
* Elements of an m-by-n matrix are stored contiguously in an array a[m*n],
* such that array element a[i+j*m] corresponds to the i'th row and the
* j'th column of the m-by-n matrix. For example, a 3-by-3 matrix is stored
* as:
*
* a[0] a[3] a[6]
* a[1] a[4] a[7]
* a[2] a[5] a[8]
*
* This is the column-major order required by LAPACK.
*
* Decompositions of a matrix A facilate the solutions of various problems.
* For example, the QR decomposition A = Q*R of a rectangular matrix A
* (where Q is orthogonal and R is upper right triangular) is useful in
* least-square approximations. Decompositions currently provided by this
* class include
*
* - LU decomposition of rectangular matrices
* - QR decomposition of rectangular matrices
* - singular value decomposition of rectangular matrices
* - eigenvalue and eigenvector decomposition of square matrices
* - Cholesky decomposition of symmetric positive-definite matrices
*
* @author Dave Hale, Colorado School of Mines
* @version 2005.12.12
*/
public class DMatrix {
/**
* Constructs an m-by-n matrix of zeros.
* @param m the number of rows.
* @param n the number of columns.
*/
public DMatrix(int m, int n) {
_m = m;
_n = n;
_a = new double[m*n];
}
/**
* Constructs an m-by-n matrix filled with the specified value.
* @param m the number of rows.
* @param n the number of columns.
* @param v the value.
*/
public DMatrix(int m, int n, double v) {
this(m,n);
fill(v,_a);
}
/**
* Constructs a matrix from the specified array. Does not copy array
* elements into a new array. Rather, this matrix simply references
* the specified array.
*
* That array contains packed columns of this matrix. In other words,
* array element a[i+j*m] corresponds to the i'th row and j'th column
* of this matrix.
* @param m the number of rows.
* @param n the number of columns.
* @param a the array.
*/
public DMatrix(int m, int n, double[] a) {
Check.argument(m*n<=a.length,"m*n <= a.length");
_m = m;
_n = n;
_a = a;
}
/**
* Constructs a matrix from the specified array. Copies array elements
* a[i][j] to the i'th row and j'th column of this matrix. In other
* words, the specified array contains matrix elements ordered by rows.
*
* The specified array must be regular. That is, each row much contain
* the same number of columns, and each column must contain the same
* number of rows.
* @param a the array.
*/
public DMatrix(double[][] a) {
Check.argument(isRegular(a),"array a is regular");
_m = a.length;
_n = a[0].length;
_a = new double[_m*_n];
set(a);
}
/**
* Constructs a copy of the specified matrix.
* @param a the matrix.
*/
public DMatrix(DMatrix a) {
this(a._m,a._n, copy(a._a));
}
/**
* Gets the number of rows in this matrix.
* @return the number of rows.
*/
public int getM() {
return _m;
}
/**
* Gets the number of rows in this matrix.
* @return the number of rows.
*/
public int getRowCount() {
return _m;
}
/**
* Gets the number of columns in this matrix.
* @return the number of columns.
*/
public int getN() {
return _n;
}
/**
* Gets the number of columns in this matrix.
* @return the number of columns.
*/
public int getColumnCount() {
return _n;
}
/**
* Gets the array in which matrix elements are stored.
* @return the array; by reference, not by copy.
*/
public double[] getArray() {
return _a;
}
/**
* Determines whether this matrix is square.
* @return true, if square; false, otherwise.
*/
public boolean isSquare() {
return _m==_n;
}
/**
* Determines whether this matrix is symmetric (and square).
* @return true, if symmetric (and square); false, otherwise.
*/
public boolean isSymmetric() {
if (!isSquare())
return false;
for (int i=0; i<_n; ++i)
for (int j=i+1; j<_n; ++j)
if (_a[i+j*_m]!=_a[j+i*_m])
return false;
return true;
}
/**
* Gets all elements of this matrix into a new array.
* @return the array.
*/
public double[][] get() {
double[][] a = new double[_m][_n];
get(a);
return a;
}
/**
* Gets all elements of this matrix into the specified array.
* @param a the array.
*/
public void get(double[][] a) {
for (int i=0; i<_m; ++i) {
double[] ai = a[i];
for (int j=0; j<_n; ++j) {
ai[j] = _a[i+j*_m];
}
}
}
/**
* Gets a matrix element.
* @param i the row index.
* @param j the column index.
* @return the element.
*/
public double get(int i, int j) {
return _a[i+j*_m];
}
/**
* Gets the specified submatrix a[i0:i1][j0:j1] of this matrix.
* @param i0 the index of first row.
* @param i1 the index of last row.
* @param j0 the index of first column.
* @param j1 the index of last column.
*/
public DMatrix get(int i0, int i1, int j0, int j1) {
checkI(i0,i1);
checkJ(j0,j1);
int m = i1-i0+1;
int n = j1-j0+1;
double[] b = new double[m*n];
for (int j=0; j=_n,"number of rows is not less than number of columns");
if (_m==_n) {
return new DMatrixLud(this).solve(b);
} else {
return new DMatrixQrd(this).solve(b);
}
}
/**
* Returns the inverse of this matrix.
* If m>n for this m-by-n matrix A, then returns the pseudo-inverse.
*/
public DMatrix inverse() {
return solve(identity(_m,_m));
}
/**
* Returns C = -A, where A is this matrix.
* @return C = -A.
*/
public DMatrix negate() {
DMatrix c = new DMatrix(_m,_n);
neg(_a,c._a);
return c;
}
/**
* Returns C = A + B, where A is this matrix.
* @param b the matrix B.
* @return C = A + B.
*/
public DMatrix plus(DMatrix b) {
DMatrix c = new DMatrix(_m,_n);
add(_a,b._a,c._a);
return c;
}
/**
* Returns A = A + B, where A is this matrix.
* @param b the matrix B.
* @return A = A + B.
*/
public DMatrix plusEquals(DMatrix b) {
add(_a,b._a,_a);
return this;
}
/**
* Returns C = A - B, where A is this matrix.
* @param b the matrix B.
* @return C = A - B.
*/
public DMatrix minus(DMatrix b) {
DMatrix c = new DMatrix(_m,_n);
sub(_a,b._a,c._a);
return c;
}
/**
* Returns A = A - B, where A is this matrix.
* @param b the matrix B.
* @return A = A - B.
*/
public DMatrix minusEquals(DMatrix b) {
sub(_a,b._a,_a);
return this;
}
/**
* Returns C = A .* B, where A is this matrix.
* The symbol .* denotes element-by-element multiplication.
* @param b the matrix B.
* @return C = A .* B.
*/
public DMatrix arrayTimes(DMatrix b) {
DMatrix c = new DMatrix(_m,_n);
mul(_a,b._a,c._a);
return c;
}
/**
* Returns A = A .* B, where A is this matrix.
* The symbol .* denotes element-by-element multiplication.
* @param b the matrix B.
* @return A = A .* B.
*/
public DMatrix arrayTimesEquals(DMatrix b) {
mul(_a,b._a,_a);
return this;
}
/**
* Returns C = A ./ B, where A is this matrix.
* The symbol ./ denotes element-by-element right division.
* @param b the matrix B.
* @return C = A ./ B.
*/
public DMatrix arrayRightDivide(DMatrix b) {
DMatrix c = new DMatrix(_m,_n);
div(_a,b._a,c._a);
return c;
}
/**
* Returns A = A ./ B, where A is this matrix.
* The symbol ./ denotes element-by-element right division.
* @param b the matrix B.
* @return A = A ./ B.
*/
public DMatrix arrayRightDivideEquals(DMatrix b) {
div(_a,b._a,_a);
return this;
}
/**
* Returns C = A .\ B, where A is this matrix.
* The symbol .\ denotes element-by-element left division.
* @param b the matrix B.
* @return C = A .\ B.
*/
public DMatrix arrayLeftDivide(DMatrix b) {
DMatrix c = new DMatrix(_m,_n);
div(b._a,_a,c._a);
return c;
}
/**
* Returns A = A .\ B, where A is this matrix.
* The symbol .\ denotes element-by-element left division.
* @param b the matrix B.
* @return A = A .\ B.
*/
public DMatrix arrayLeftDivideEquals(DMatrix b) {
div(b._a,_a,_a);
return this;
}
/**
* Returns C = A * s, where A is this matrix, and s is a scalar.
* @param s the scalar s.
* @return C = A * s.
*/
public DMatrix times(double s) {
DMatrix c = new DMatrix(_m,_n);
mul(_a,s,c._a);
return c;
}
/**
* Returns A = A * s, where A is this matrix, and s is a scalar.
* @param s the scalar s.
* @return A = A * s.
*/
public DMatrix timesEquals(double s) {
mul(_a,s,_a);
return this;
}
/**
* Returns C = A * B, where A is this matrix. The number of columns in
* this matrix A must equal the number of rows in the specified matrix B.
* @param b the matrix B.
* @return C = A * B.
*/
public DMatrix times(DMatrix b) {
Check.argument(_n==b._m,
"number of columns in A equals number of rows in B");
DMatrix c = new DMatrix(_m,b._n);
_blas.dgemm("N","N",_m,b._n,_n,1.0,_a,_m,b._a,b._m,1.0,c._a,c._m);
return c;
}
/**
* Returns C = A * B', where A is this matrix and B' is B transposed.
* The number of columns in this matrix A must equal the number of
* columns in the specified matrix B.
* @param b the matrix B.
* @return C = A * B'.
*/
public DMatrix timesTranspose(DMatrix b) {
Check.argument(_n==b._n,
"number of columns in A equals number of columns in B");
DMatrix c = new DMatrix(_m,b._m);
_blas.dgemm("N","T",_m,b._n,_n,1.0,_a,_m,b._a,b._m,1.0,c._a,c._m);
return c;
}
/**
* Returns C = A' * B, where A' is this matrix transposed.
* The number of rows in this matrix A must equal the number of
* rows in the specified matrix B.
* @param b the matrix B.
* @return C = A' * B.
*/
public DMatrix transposeTimes(DMatrix b) {
Check.argument(_m==b._m,
"number of rows in A equals number of rows in B");
DMatrix c = new DMatrix(_n,b._n);
_blas.dgemm("T","N",_n,b._n,_m,1.0,_a,_m,b._a,b._m,1.0,c._a,c._m);
return c;
}
/**
* Returns a new matrix with random elements. The distribution of the
* random numbers is uniform in the interval [0,1).
* @param m the number of rows.
* @param n the number of columns.
* @return the random matrix.
*/
public static DMatrix random(int m, int n) {
DMatrix x = new DMatrix(m,n);
rand(x._a);
return x;
}
/**
* Returns a new square matrix with random elements. The distribution
* of the random numbers is uniform in the interval [0,1).
* @param n the number of rows and columns.
* @return the random matrix.
*/
public static DMatrix random(int n) {
return DMatrix.random(n,n);
}
/**
* Returns a new identity matrix.
* @param m the number of rows.
* @param n the number of columns.
* @return the identity matrix.
*/
public static DMatrix identity(int m, int n) {
DMatrix x = new DMatrix(m,n);
double[] xa = x._a;
int mn = min(m,n);
for (int i=0; i>>32));
}
return h;
}
public String toString() {
String ls = System.getProperty("line.separator");
StringBuilder sb = new StringBuilder();
String[][] s = format(_m,_n,_a);
int max = maxlen(s);
String format = "%"+max+"s";
sb.append("[[");
int ncol = 77/(max+2);
if (ncol>=5)
ncol = (ncol/5)*5;
for (int i=0; i<_m; ++i) {
int nrow = 1+(_n-1)/ncol;
if (i>0)
sb.append(" [");
for (int irow=0,j=0; irow=_m)
Check.argument(0<=i && i<_m,"row index i="+i+" is in bounds");
}
private void checkJ(int j) {
if (j<0 || j>=_n)
Check.argument(0<=j && j<_n,"column index j="+j+" is in bounds");
}
private void checkI(int i0, int i1) {
checkI(i0); checkI(i1); Check.argument(i0<=i1,"i0<=i1");
}
private void checkJ(int j0, int j1) {
checkJ(j0); checkJ(j1); Check.argument(j0<=j1,"j0<=j1");
}
// Used in implementation of toString() above.
private static String[][] format(int m, int n, double[] d) {
int pg = 6;
String fg = "% ."+pg+"g";
int pemax = -1;
int pfmax = -1;
for (int i=0; i7)?ls-7:0;
if (pemax=0)?ls-1-ip:0;
if (pfmax=0) {
if (pfmax>pg-1)
pfmax = pg-1;
int pe = (pemax>pfmax)?pemax:pfmax;
f = "% ."+pe+"e";
} else {
int pf = pfmax;
f = "% ."+pf+"f";
}
for (int i=0; i0) {
while (ibeg>0 && s.charAt(ibeg-1)=='0')
--ibeg;
if (ibeg>0 && s.charAt(ibeg-1)=='.')
--ibeg;
}
if (ibeg