/** * The O(n^3) multiplication algorithm that is here for if BLAS isn't * available, or if the dimensions wouldn't allow 1-d array indexing. * * @param m The right matrix to multiply by this * @return The result of multiplying this times m */ private Matrix slowMult( final DenseMatrix m) { DenseMatrix result = new DenseMatrix(getNumRows(), m.getNumColumns()); for (int i = 0; i < getNumRows(); i++) { for (int j = 0; j < m.getNumColumns(); j++) { double val = 0; for (int k = 0; k < getNumColumns(); k++) { val += get(i, k) * m.get(k, j); } result.setElement(i, j, val); } } return result; }
@Override public final Matrix times( final DenseMatrix other) { this.assertMultiplicationDimensions(other); // TODO: Make sure this BLAS is truly faster than slow version if (canUseBlasForMult(getNumRows(), getNumColumns(), other.getNumRows(), other.getNumColumns())) { double[] output = new double[getNumRows() * other.getNumColumns()]; BLAS.getInstance().dgemm("N", "N", getNumRows(), other.getNumColumns(), getNumColumns(), 1.0, this.toBlas(), getNumRows(), other.toBlas(), other.getNumRows(), 0.0, output, getNumRows()); return createFromBlas(output, getNumRows(), other.getNumColumns()); } else { return slowMult(other); } }
@Override public final void plusEquals( final DenseMatrix other) { this.assertSameDimensions(other); final int numRows = this.getNumRows(); final int numColumns = this.getNumColumns(); for (int i = 0; i < numRows; ++i) { for (int j = 0; j < numColumns; ++j) { this.rows[i].values[j] += other.rows[i].values[j]; } } }
/** * Initializes Q and R to the correct size (based on A), but leaves them * as zeroes. * * @param A The matrix to be factored into the Q and R stored herein. */ private QR( final DenseMatrix A) { Q = new DenseMatrix(A.getNumRows(), A.getNumRows()); R = new DenseMatrix(A.getNumRows(), A.getNumColumns()); }
/** * Reads the input double array (dense, column-major as output by BLAS) and * creates a new dense matrix with those values. * * @param d The BLAS representation of the matrix (must be length * numRows*numCols) * @param numRows The number of rows in d * @param numCols The number of columns in d * @return A new DenseMatrix with the specified size and elements stored in * d (column major order) */ static DenseMatrix createFromBlas( final double d[], final int numRows, final int numCols) { DenseMatrix result = new DenseMatrix(numRows, numCols); result.fromBlas(d, numRows, numCols); return result; }
@Override final public Matrix solve( final Matrix B) { checkSolveDimensions(B); if (!isSquare()) { throw new IllegalStateException("Solve only works on square " + "matrices (this is " + getNumRows() + " x " + getNumColumns()); } QR qr = qrDecompose(); // I'll only use it as the transpose qr.Q = qr.Q.transpose(); final int numColumns = B.getNumColumns(); DenseMatrix X = new DenseMatrix(getNumColumns(), numColumns); for (int i = 0; i < numColumns; ++i) { X.setColumn(i, upperTriangularSolve(qr.R, qr.Q.times(B.getColumn(i)))); } return X; }
@Override final public Vector solve( final Vector b) { checkSolveDimensions(b); if (!isSquare()) { throw new IllegalStateException("Solve only works on square " + "matrices (this is " + getNumRows() + " x " + getNumColumns()); } QR qr = qrDecompose(); return upperTriangularSolve(qr.R, qr.Q.transpose().times(b)); }
int m = getNumRows(); int n = getNumColumns(); if ((m == 0) || (n == 0)) double[] A = this.toBlas(); int ipivdim = Math.min(m, n); int[] ipiv = new int[ipivdim]; result.U.setElement(i, j, blasElement(i, j, A, m, n)); result.L.setElement(i, j, blasElement(i, j, A, m, n));
@Override public final Matrix times( final DiagonalMatrix other) { this.assertMultiplicationDimensions(other); final int numRows = this.getNumRows(); final DenseVector[] resultRows = new DenseVector[numRows]; for (int i = 0; i < numRows; ++i) { resultRows[i] = (DenseVector) other.preTimes(this.rows[i]); } return new DenseMatrix(resultRows); }
int m = getNumRows(); int n = getNumColumns(); if ((m == 0) || (n == 0)) double[] A = this.toBlas(); int lda = m; int mindim = Math.min(m, n); ((DenseMatrix) result.U).fromBlas(u, ldu, ucol); ((DenseMatrix) result.V).fromBlas(vt, ldvt, n); result.V = result.V.transpose(); for (int i = 0; i < mindim; ++i)
@Override public final void plusEquals( final DiagonalMatrix other) { this.assertSameDimensions(other); final int numRows = this.getNumRows(); for (int i = 0; i < numRows; ++i) { this.rows[i].values[i] += other.get(i, i); } }
/** * {@inheritDoc} * * NOTE: This is inclusive on both end points. * @param minRow {@inheritDoc} * @param maxRow {@inheritDoc} * @param minColumn {@inheritDoc} * @param maxColumn {@inheritDoc} * @return {@inheritDoc} */ @Override final public Matrix getSubMatrix( final int minRow, final int maxRow, final int minColumn, final int maxColumn) { checkSubmatrixRange(minRow, maxRow, minColumn, maxColumn); DenseMatrix result = new DenseMatrix( maxRow - minRow + 1, maxColumn - minColumn + 1); for (int i = minRow; i <= maxRow; ++i) { for (int j = minColumn; j <= maxColumn; ++j) { result.set(i - minRow, j - minColumn, this.rows[i].values[j]); } } return result; }
/** * Loads the input values from the BLAS-ordered (dense, column-major) vector * d into this. * * @param d The values to load into this * @param numRows The number of rows in d * @param numCols The number of columns in d * @throws IllegalArgumentException if the input dimensions don't match * this's dimensions (that's why this is a private method -- it should only * be called by people who know what they're doing). */ private void fromBlas( final double d[], final int numRows, final int numCols) { if ((getNumRows() != numRows) || (getNumColumns() != numCols)) { throw new IllegalArgumentException("Unable to convert from BLAS of " + "different size: (" + getNumRows() + " != " + numRows + ") || (" + getNumColumns() + " != " + numCols + ")"); } for (int i = 0; i < numRows; ++i) { for (int j = 0; j < numCols; ++j) { this.rows[i].values[j] = blasElement(i, j, d, numRows, numCols); } } }
@Override public final Matrix times( final DenseMatrix other) { this.assertMultiplicationDimensions(other); final DenseVector[] rows = new DenseVector[diagonal.length]; for (int i = 0; i < diagonal.length; ++i) { DenseVector v = other.row(i); rows[i] = (DenseVector) v.scale(diagonal[i]); } return new DenseMatrix(rows); }
/** * {@inheritDoc} * * NOTE: Upon completion this is in the compressed Yale format. * * This is implemented by creating a new dense matrix version of this and * calling its logDeterminant method -- It requires factoring the matrix * which is going to be memory intense anyway. * * @return {@inheritDoc} */ @Override final public ComplexNumber logDeterminant() { if (!isCompressed()) { compress(); } // NOTE: I'm using LU decomposition to do this, so ... to Dense! return (new DenseMatrix(this)).logDeterminant(); }
/** * Creates a deep copy of m into a DenseMatrix and returns it. * * @param m The matrix to copy * @return The DenseMatrix deep copy of m. */ @Override final public DenseMatrix copyMatrix( final Matrix m) { return new DenseMatrix(m); }
@Override public final void minusEquals( final SparseMatrix other) { this.assertSameDimensions(other); if (!other.isCompressed()) { other.compress(); } double[] ovals = other.getValues(); int[] ocolIdxs = other.getColumnIndices(); int[] ofirstRows = other.getFirstInRows(); int rownum = 0; for (int i = 0; i < ovals.length; ++i) { while (ofirstRows[rownum + 1] <= i) { ++rownum; } this.rows[rownum].values[ocolIdxs[i]] -= ovals[i]; } }
/** * {@inheritDoc} * * NOTE: Upon completion this is in the compressed Yale format. * * This is implemented by creating a new dense matrix version of this and * calling its inverse method -- inverting a sparse matrix is likely to * generate a dense matrix anyway. We would recommend using an iterative * solver (like Conjugate Gradient). * * @return {@inheritDoc} */ @Override final public Matrix inverse() { if (!isCompressed()) { compress(); } // NOTE: Inverting a sparse matrix generally creates a dense matrix // anyway. Sparse matrices make the most sense for iterative solvers return (new DenseMatrix(this)).inverse(); }
/** * {@inheritDoc} * * NOTE: Upon completion this is in the compressed Yale format. */ @Override public final void dotTimesEquals( final DenseMatrix other) { this.assertSameDimensions(other); if (!isCompressed()) { compress(); } // The shortcut I'll take here is that few of the dense matrix's values // are zero where I'm not int rownum = 0; for (int i = 0; i < values.length; ++i) { while (i >= firstIndicesForRows[rownum + 1]) { ++rownum; } values[i] *= other.get(rownum, columnIndices[i]); } }
@Override final public Matrix solve( final Matrix B) { checkSolveDimensions(B); if (!isSquare()) { throw new IllegalStateException("Solve only works on square " + "matrices (this is " + getNumRows() + " x " + getNumColumns()); } QR qr = qrDecompose(); // I'll only use it as the transpose qr.Q = qr.Q.transpose(); final int numColumns = B.getNumColumns(); DenseMatrix X = new DenseMatrix(getNumColumns(), numColumns); for (int i = 0; i < numColumns; ++i) { X.setColumn(i, upperTriangularSolve(qr.R, qr.Q.times(B.getColumn(i)))); } return X; }