Modules | C++
Linalg
Sparse symmetric positive definite solves — one pattern, many value sets, many right-hand sides.
The Linalg module solves sparse symmetric positive definite systems by Cholesky factorization. Its three steps are split by what they read: the pattern is analyzed once, every value set on that pattern is factored, and every factor solves as many right-hand sides as are given. The Dirichlet field and the fill's fairing solve through it.
Include the module with:
#include <trueform/linalg.hpp>
Types
| Type | Holds |
|---|---|
tf::symmetric_lower<Index> | A symmetric matrix by its lower triangle: row_offsets, columns and values, row i holding its columns j <= i ascending, so its diagonal is its last entry |
tf::cholesky_pattern<Index> | The factor's structure: the fill-reducing elimination order, the supernodes, their rows and the schedule the numeric factor takes them in |
tf::cholesky_factor | The factor's values, one dense block per supernode where the pattern places it |
tf::cholesky_scratch<Index> | The workspace every step walks on, owned by the caller and reused across calls |
Functions
cholesky.cpp
// A 1D Laplacian with a positive diagonal: row i holds (i - 1, i)
const std::size_t n = 1000;
tf::symmetric_lower<int> matrix;
matrix.row_offsets.push_back(0);
for (std::size_t i = 0; i < n; ++i) {
if (i > 0) {
matrix.columns.push_back(int(i - 1));
matrix.values.push_back(-1.0);
}
matrix.columns.push_back(int(i));
matrix.values.push_back(2.5);
matrix.row_offsets.push_back(int(matrix.columns.size()));
}
tf::cholesky_scratch<int> scratch;
tf::cholesky_pattern<int> pattern;
tf::analyze_cholesky(matrix, scratch, pattern); // once per pattern
tf::cholesky_factor factor;
if (tf::factor_cholesky(pattern, matrix.values, scratch, factor)) {
// Two right-hand sides, interleaved: entry (i, c) at i * 2 + c
tf::buffer<double> rhs;
rhs.allocate(2 * n);
// ... fill rhs ...
tf::solve_cholesky(pattern, factor, rhs, 2, scratch); // rhs receives x
}
| Function | Returns | Does |
|---|---|---|
tf::analyze_cholesky(matrix, scratch, pattern) | bool | Reads the factor's structure from the matrix's pattern alone; false where the index type cannot name the factor's values |
tf::factor_cholesky(pattern, values, scratch, factor) | bool | Factors one value set, values one per entry of the matrix the pattern was read from, in its order; false where a pivot is not positive or not finite |
tf::solve_cholesky(pattern, factor, rhs, r, scratch) | — | Solves r interleaved right-hand sides in place, in the matrix's own row order |
A refactor of new values on the same pattern skips the analysis, and a factor serves any number of solves.
The factor is the same bits at every thread count, partition and run, and so is every solve; a right-hand side's answer is the same bits whatever the number of right-hand sides solved beside it. A factor is not the bits of a column-by-column factorization of the same matrix, which sums in another order — the two agree to the condition number times the unit roundoff.