/*
* Copyright (C) 2003-2006 Bjørn-Ove Heimsund
*
* This file is part of MTJ.
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation; either version 2.1 of the License, or (at your
* option) any later version.
*
* This library is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
* for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
package no.uib.cipr.matrix;
import java.io.IOException;
import no.uib.cipr.matrix.io.MatrixInfo;
import no.uib.cipr.matrix.io.MatrixSize;
import no.uib.cipr.matrix.io.MatrixVectorReader;
import org.netlib.blas.BLAS;
import org.netlib.lapack.LAPACK;
import org.netlib.util.intW;
/**
* Dense matrix. It is a good all-round matrix structure, with fast access and
* efficient algebraic operations. The matrix
* <p>
* <table border="1">
* <tr>
* <td>a<sub>11</sub></td>
* <td>a<sub>12</sub></td>
* <td>a<sub>13</sub></td>
* <td>a<sub>14</sub></td>
* </tr>
* <tr>
* <td>a<sub>21</sub></td>
* <td>a<sub>22</sub></td>
* <td>a<sub>23</sub></td>
* <td>a<sub>24</sub></td>
* </tr>
* <tr>
* <td>a<sub>31</sub></td>
* <td>a<sub>32</sub></td>
* <td>a<sub>33</sub></td>
* <td>a<sub>34</sub></td>
* </tr>
* <tr>
* <td>a<sub>41</sub></td>
* <td>a<sub>42</sub></td>
* <td>a<sub>43</sub></td>
* <td>a<sub>44</sub></td>
* </tr>
* </table>
* </p>
* <p>
* is stored column major in a single array, as follows:
* </p>
* <p>
* <table border="1">
* <tr>
* <td>a<sub>11</sub></td>
* <td>a<sub>21</sub></td>
* <td>a<sub>31</sub></td>
* <td>a<sub>41</sub></td>
* <td>a<sub>12</sub></td>
* <td>a<sub>22</sub></td>
* <td>a<sub>32</sub></td>
* <td>a<sub>42</sub></td>
* <td>a<sub>13</sub></td>
* <td>a<sub>23</sub></td>
* <td>a<sub>33</sub></td>
* <td>a<sub>43</sub></td>
* <td>a<sub>14</sub></td>
* <td>a<sub>24</sub></td>
* <td>a<sub>34</sub></td>
* <td>a<sub>44</sub></td>
* </tr>
* </table>
* </p>
*/
public class DenseMatrix extends AbstractDenseMatrix {
/**
* Constructor for DenseMatrix
*
* @param r
* Reader to get the matrix from
*/
public DenseMatrix(MatrixVectorReader r) throws IOException {
// Start with a zero-sized matrix
super(0, 0);
// Get matrix information. Use the header if present, else use a safe
// default
MatrixInfo info = null;
if (r.hasInfo())
info = r.readMatrixInfo();
else
info = new MatrixInfo(true, MatrixInfo.MatrixField.Real,
MatrixInfo.MatrixSymmetry.General);
MatrixSize size = r.readMatrixSize(info);
// Resize the matrix to correct size
numRows = size.numRows();
numColumns = size.numColumns();
data = new double[numRows * numColumns];
// Check that the matrix is in an acceptable format
if (info.isPattern())
throw new UnsupportedOperationException(
"Pattern matrices are not supported");
if (info.isComplex())
throw new UnsupportedOperationException(
"Complex matrices are not supported");
// Read the entries, in either coordinate or array format
if (info.isCoordinate()) {
// Read coordinate data
int nz = size.numEntries();
int[] row = new int[nz];
int[] column = new int[nz];
double[] entry = new double[nz];
r.readCoordinate(row, column, entry);
// Shift indices from 1-offset to 0-offset
r.add(-1, row);
r.add(-1, column);
// Store them
for (int i = 0; i < nz; ++i)
set(row[i], column[i], entry[i]);
} else
// info.isArray()
r.readArray(data);
// Put in missing entries from symmetry or skew symmetry
if (info.isSymmetric())
for (int i = 0; i < numRows; ++i)
for (int j = 0; j < i; ++j)
set(j, i, get(i, j));
else if (info.isSkewSymmetric())
for (int i = 0; i < numRows; ++i)
for (int j = 0; j < i; ++j)
set(j, i, -get(i, j));
}
/**
* Constructor for DenseMatrix
*
* @param numRows
* Number of rows
* @param numColumns
* Number of columns
*/
public DenseMatrix(int numRows, int numColumns) {
super(numRows, numColumns);
}
/**
* Constructor for DenseMatrix
*
* @param A
* Matrix to copy. A deep copy is made
*/
public DenseMatrix(Matrix A) {
super(A);
}
/**
* Constructor for DenseMatrix
*
* @param A
* Matrix to copy contents from
* @param deep
* If true, <code>A</code> is copied, else a shallow copy is
* made and the matrices share underlying storage. For this,
* <code>A</code> must be a dense matrix
*/
public DenseMatrix(Matrix A, boolean deep) {
super(A, deep);
}
/**
* Constructor for DenseMatrix. Builds the matrix from a vector
*
* @param x
* Vector to copy from. This will form this matrix' single column
* @param deep
* If true, x is copied, if false, the internal storage of this
* matrix is the same as that of the vector. In that case,
* <code>x</code> must be a <code>DenseVector</code>
*/
public DenseMatrix(Vector x, boolean deep) {
super(x.size(), 1);
if (deep)
for (VectorEntry e : x)
set(e.index(), 0, e.get());
else {
if (!(x instanceof DenseVector))
throw new IllegalArgumentException("x must be a DenseVector");
data = ((DenseVector) x).getData();
}
}
/**
* Constructor for DenseMatrix. Builds the matrix from a vector
*
* @param x
* The vector which forms this matrix' single column. It is
* copied, not referenced
*/
public DenseMatrix(Vector x) {
this(x, true);
}
/**
* Constructor for DenseMatrix. Builds the matrix from vectors. Each vector
* will correspond to a column of the matrix
*
* @param x
* Vectors which forms the columns of this matrix. Every vector
* must have the same size
*/
public DenseMatrix(Vector[] x) {
super(x[0].size(), x.length);
// Ensure correct sizes
for (Vector v : x)
if (v.size() != numRows)
throw new IllegalArgumentException(
"All vectors must be of the same size");
// Copy the contents
for (int j = 0; j < x.length; ++j)
for (VectorEntry e : x[j])
set(e.index(), j, e.get());
}
/**
* Constructor for DenseMatrix. Copies from the passed array
*
* @param values
* Arrays to copy from. Every sub-array must have the same size
*/
public DenseMatrix(double[][] values) {
super(values.length, values[0].length);
// Copy the contents
for (int i = 0; i < values.length; ++i) {
if (values[i].length != numColumns)
throw new IllegalArgumentException("Array cannot be jagged");
for (int j = 0; j < values[i].length; ++j)
set(i, j, values[i][j]);
}
}
@Override
public DenseMatrix copy() {
return new DenseMatrix(this);
}
@Override
void copy(Matrix A) {
for (MatrixEntry e : A)
set(e.row(), e.column(), e.get());
}
@Override
public Matrix multAdd(double alpha, Matrix B, Matrix C) {
if (!(B instanceof DenseMatrix) || !(C instanceof DenseMatrix))
return super.multAdd(alpha, B, C);
checkMultAdd(B, C);
double[] Bd = ((DenseMatrix) B).getData(), Cd = ((DenseMatrix) C)
.getData();
BLAS.getInstance().dgemm(Transpose.NoTranspose.netlib(), Transpose.NoTranspose.netlib(),
C.numRows(), C.numColumns(), numColumns, alpha, data,
Math.max(1, numRows), Bd, Math.max(1, B.numRows()), 1, Cd,
Math.max(1, C.numRows()));
return C;
}
@Override
public Matrix transAmultAdd(double alpha, Matrix B, Matrix C) {
if (!(B instanceof DenseMatrix) || !(C instanceof DenseMatrix))
return super.transAmultAdd(alpha, B, C);
checkTransAmultAdd(B, C);
double[] Bd = ((DenseMatrix) B).getData(), Cd = ((DenseMatrix) C)
.getData();
BLAS.getInstance().dgemm(Transpose.Transpose.netlib(), Transpose.NoTranspose.netlib(),
C.numRows(), C.numColumns(), numRows, alpha, data,
Math.max(1, numRows), Bd, Math.max(1, B.numRows()), 1, Cd,
Math.max(1, C.numRows()));
return C;
}
@Override
public Matrix transBmultAdd(double alpha, Matrix B, Matrix C) {
if (!(B instanceof DenseMatrix) || !(C instanceof DenseMatrix))
return super.transBmultAdd(alpha, B, C);
checkTransBmultAdd(B, C);
double[] Bd = ((DenseMatrix) B).getData(), Cd = ((DenseMatrix) C)
.getData();
BLAS.getInstance().dgemm(Transpose.NoTranspose.netlib(), Transpose.Transpose.netlib(),
C.numRows(), C.numColumns(), numColumns, alpha, data,
Math.max(1, numRows), Bd, Math.max(1, B.numRows()), 1, Cd,
Math.max(1, C.numRows()));
return C;
}
@Override
public Matrix transABmultAdd(double alpha, Matrix B, Matrix C) {
if (!(B instanceof DenseMatrix) || !(C instanceof DenseMatrix))
return super.transABmultAdd(alpha, B, C);
checkTransABmultAdd(B, C);
double[] Bd = ((DenseMatrix) B).getData(), Cd = ((DenseMatrix) C)
.getData();
BLAS.getInstance().dgemm(Transpose.Transpose.netlib(), Transpose.Transpose.netlib(),
C.numRows(), C.numColumns(), numRows, alpha, data,
Math.max(1, numRows), Bd, Math.max(1, B.numRows()), 1, Cd,
Math.max(1, C.numRows()));
return C;
}
@Override
public Matrix rank1(double alpha, Vector x, Vector y) {
if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
return super.rank1(alpha, x, y);
checkRank1(x, y);
double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
.getData();
BLAS.getInstance().dger(numRows, numColumns, alpha, xd, 1, yd, 1, data,
Math.max(1, numRows));
return this;
}
@Override
public Vector multAdd(double alpha, Vector x, Vector y) {
if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
return super.multAdd(alpha, x, y);
checkMultAdd(x, y);
double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
.getData();
BLAS.getInstance().dgemv(Transpose.NoTranspose.netlib(), numRows, numColumns,
alpha, data, Math.max(numRows, 1), xd, 1, 1, yd, 1);
return y;
}
@Override
public Vector transMultAdd(double alpha, Vector x, Vector y) {
if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
return super.transMultAdd(alpha, x, y);
checkTransMultAdd(x, y);
double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
.getData();
BLAS.getInstance().dgemv(Transpose.Transpose.netlib(), numRows, numColumns, alpha,
data, Math.max(numRows, 1), xd, 1, 1, yd, 1);
return y;
}
@Override
public Matrix solve(Matrix B, Matrix X) {
// We allow non-square matrices, as we then use a least-squares solver
if (numRows != B.numRows())
throw new IllegalArgumentException("numRows != B.numRows() ("
+ numRows + " != " + B.numRows() + ")");
if (numColumns != X.numRows())
throw new IllegalArgumentException("numColumns != X.numRows() ("
+ numColumns + " != " + X.numRows() + ")");
if (X.numColumns() != B.numColumns())
throw new IllegalArgumentException(
"X.numColumns() != B.numColumns() (" + X.numColumns()
+ " != " + B.numColumns() + ")");
if (isSquare())
return LUsolve(B, X);
else
return QRsolve(B, X, Transpose.NoTranspose);
}
@Override
public Vector solve(Vector b, Vector x) {
DenseMatrix B = new DenseMatrix(b, false), X = new DenseMatrix(x, false);
solve(B, X);
return x;
}
@Override
public Matrix transSolve(Matrix B, Matrix X) {
// We allow non-square matrices, as we then use a least-squares solver
if (numColumns != B.numRows())
throw new IllegalArgumentException("numColumns != B.numRows() ("
+ numColumns + " != " + B.numRows() + ")");
if (numRows != X.numRows())
throw new IllegalArgumentException("numRows != X.numRows() ("
+ numRows + " != " + X.numRows() + ")");
if (X.numColumns() != B.numColumns())
throw new IllegalArgumentException(
"X.numColumns() != B.numColumns() (" + X.numColumns()
+ " != " + B.numColumns() + ")");
return QRsolve(B, X, Transpose.Transpose);
}
@Override
public Vector transSolve(Vector b, Vector x) {
DenseMatrix B = new DenseMatrix(b, false), X = new DenseMatrix(x, false);
transSolve(B, X);
return x;
}
Matrix LUsolve(Matrix B, Matrix X) {
if (!(X instanceof DenseMatrix))
throw new UnsupportedOperationException("X must be a DenseMatrix");
double[] Xd = ((DenseMatrix) X).getData();
X.set(B);
int[] piv = new int[numRows];
intW info = new intW(0);
LAPACK.getInstance().dgesv(numRows, B.numColumns(),
data.clone(), Matrices.ld(numRows), piv, Xd, Matrices.ld(numRows), info);
if (info.val > 0)
throw new MatrixSingularException();
else if (info.val < 0)
throw new IllegalArgumentException();
return X;
}
Matrix QRsolve(Matrix B, Matrix X, Transpose trans) {
int nrhs = B.numColumns();
// Allocate temporary solution matrix
DenseMatrix Xtmp = new DenseMatrix(Math.max(numRows, numColumns), nrhs);
int M = trans == Transpose.NoTranspose ? numRows : numColumns;
for (int j = 0; j < nrhs; ++j)
for (int i = 0; i < M; ++i)
Xtmp.set(i, j, B.get(i, j));
double[] newData = data.clone();
// Query optimal workspace
double[] work = new double[1];
intW info = new intW(0);
LAPACK.getInstance().dgels(trans.netlib(), numRows, numColumns, nrhs,
newData, Matrices.ld(numRows), Xtmp.getData(), Matrices.ld(numRows, numColumns),
work, -1, info);
// Allocate workspace
int lwork = -1;
if (info.val != 0)
lwork = Math.max(1, Math.min(numRows, numColumns)
+ Math.max(Math.min(numRows, numColumns), nrhs));
else
lwork = Math.max((int) work[0], 1);
work = new double[lwork];
// Compute the factorization
info.val = 0;
LAPACK.getInstance().dgels(trans.netlib(), numRows, numColumns, nrhs,
newData, Matrices.ld(numRows), Xtmp.getData(), Matrices.ld(numRows, numColumns),
work, lwork, info);
if (info.val < 0)
throw new IllegalArgumentException();
// Extract the solution
int N = trans == Transpose.NoTranspose ? numColumns : numRows;
for (int j = 0; j < nrhs; ++j)
for (int i = 0; i < N; ++i)
X.set(i, j, Xtmp.get(i, j));
return X;
}
}