Linear algebra utilities
Construct and convert
SparseMatrix<T> identityMatrix(size_t N)
Construct and N x N identity matrix of the requested type.
void shiftDiagonal(SparseMatrix<T>& m, T shiftAmount = 1e-4)
Shift the diagonal of matrix, by adding A + shiftDiagonal * identityMatrix().
SparseMatrix<T> verticalStack(const std::vector<SparseMatrix<T>>& mats)
Vertically stack sparse matrices like
all matrices must have the same number of columns.
Example:
SparseMatrix<double> matA = /* 35 x N */
SparseMatrix<double> matB = /* 10 x N */
SparseMatrix<double> matC = /* 2 x N */
SparseMatrix<double> stacked = verticalStack<double>({matA, matB, matC});
SparseMatrix<T> horizontalStack(const std::vector<SparseMatrix<T>>& mats)
Vertically stack sparse matrices like
all matrices must have the same number of rows..
Example:
SparseMatrix<double> matA = /* N x 35 */
SparseMatrix<double> matB = /* N x 10 */
SparseMatrix<double> matC = /* N x 2 */
SparseMatrix<double> stacked = horizontalStack<double>({matA, matB, matC});
SparseMatrix<double> complexToReal(const SparseMatrix<std::complex<double>>& m)
Convert an N x M complex matrix to a 2N x 2M real matrix, expanding each complex component in to a 2 x 2 block to evaluate the complex product.
Vector<double> complexToReal(const Vector<std::complex<double>>& v)
Convert an length N complex vector to a length 2N real vector, expanding each complex component in to consecutive real and imaginary components.
Validate matrix properties
void checkFinite(const Eigen::Matrix<>& m)
Verify that all entries in an matrix are finite, throwing if not. Defined for all Eigen matrix, vector, and sparse matrix types.
void checkSymmetric(const Eigen::SparseMatrix<>& m, double absoluteEPS=-1.)
Verify that a matrix is symmetric, throwing if not. Defined for all Eigen sparse matrix types.
absoluteEPS is an epsilon to use for the element-wise comparison test. If the default value of -1 is given, a reasonable epsilon is automatically computed from the matrix entries.
void checkHermitian(const Eigen::SparseMatrix<>& m, double absoluteEPS=-1.)
Verify that a matrix is Hermitian, throwing if not. Defined for all Eigen sparse matrix types.
For real matrices, identical to check symmetric.
absoluteEPS is an epsilon to use for the element-wise comparison test. If the default value of -1 is given, a reasonable epsilon is automatically computed from the matrix entries.
Block decomposition
These routines assist with decomposing a square matrix in to interleaved submatrix blocks, where the blocks might not necessarily be contiguous. One common usage is extracting boundary components of a finite element matrix to apply boundary conditions, as in the example below.
Example usage:
// Hypothetical input data
SparseMatrix<double> mat = /* your square matrix */;
size_t N = mat.rows();
size_t NBoundary = /* ... */;
Vector<double> rhsVals = Vector<double>::Zero(N); // rhs for the system
Vector<double> bcVals = Vector<double>::Ones(NBoundary); // boundary values at
// some nodes
// Build the membership vector, which indicates which entries should be separated
// in to set "A" (others are in "B")
Vector<bool> setAMembership(N);
for(size_t i = 0; i < N; i++) {
if(/* element i is boundary */) {
setAMembership(i) = true;
} else {
setAMembership(i) = false;
}
}
// Construct the decomposition
BlockDecompositionResult<double> decomp =
blockDecomposeSquare(mat, setAMembership, true);
// The four sub-blocks of the matrix are now in
// decomp.AA, decomp.AB, decomp.BA, decomp.BB
// Split up the rhs vector
Vector<double> rhsValsA, rhsValsB;
decomposeVector(decomp, rhsVals, rhsValsA, rhsValsB);
// Solve problem
Vector<double> combinedRHS = rhsValsA - decomp.AB * bcVals;
Vector<double> Aresult = solve(decomp.AA, combinedRHS);
// Combine the two boundary conditions and interior solution to a full vector
Vector<double> result = reassembleVector(decomp, Aresult, bcVals);
BlockDecompositionResult<T> blockDecomposeSquare(const SparseMatrix<T>& m, const Vector<bool>& Aset, bool buildBuildBside = true)
Build a block decomposition of a matrix.
void decomposeVector(BlockDecompositionResult<T>& decomp, const Vector<T>& vec, Vector<T>& vecAOut, Vector<T>& vecBOut)
Use an existing block decomposition to partition a vector.
Vector<T> reassembleVector(BlockDecompositionResult<T>& decomp, const Vector<T>& vecA, const Vector<T>& vecB)
Use an existing block decomposition to build a vector from partitioned pieces.
Saving and loading matrices
These routines save and load matrices from ASCII files.
Example
SparseMatrix<double> spMat = loadSparseMatrix<double>(sparseFilename);
DenseMatrix<int> dMat = loadDenseMatrix<int>(denseFilename);
/* do something with matrices */
saveSparseMatrix(newSparseFilename, spMat);
saveDenseMatrix(newDenseFilename, dMat);
Sparse matrices
Sparse matrices are represented as ASCII files.
Each line gives the row, column, and value of an entry in the matrix.
The rows and columns are 1-indexed, following the matlab convention, and
can be loaded into matlab via spconvert.
Each file starts with a comment of the form # sparse rows cols giving the matrix dimensions.
Sparse matrices are 1-indexed
Example file
# sparse 1858 1858
1 1 3.969412682295582
2 1 -1.010447934542366
3 1 -0.5162236467210192
...
void saveSparseMatrix(std::string filename, SparseMatrix<T>& matrix);
Writes matrix to the file filename.
void saveSparseMatrix(std::ostream& out, SparseMatrix<T>& matrix);
Writes matrix to the stream out.
SparseMatrix<T> loadSparseMatrix(std::string filename);
Read a sparse matrix from file filename.
SparseMatrix<T> loadSparseMatrix(std::istream& in);
Read a sparse matrix from the stream in.
Dense matrices
Dense matrices are represented as ASCII files with one row per matrix row.
Each file starts with a comment of the form # dense rows cols giving the matrix dimensions.
Example file
# dense 1858 3
-0.095724 -0.9734159999999999 -0.244313
-0.105421 -0.984621 -0.301998
-0.003629 -0.979971 -0.245996
...
void saveDenseMatrix(std::string filename, DenseMatrix<T>& matrix);
Writes matrix to the file filename.
void saveDenseMatrix(std::ostream& out, DenseMatrix<T>& matrix);
Writes matrix to the stream out.
DenseMatrix<T> loadDenseMatrix(std::string filename);
Read a dense matrix from file filename.
DenseMatrix<T> loadDenseMatrix(std::istream& in);
Read a dense matrix from the stream in.