pub struct SparseLdlt { /* private fields */ }Expand description
An L D Lᵀ factorization of a symmetric matrix.
L is stored in CSC by column with an implicit unit diagonal (only the strictly
lower entries are kept); d is the signed diagonal of D.
Implementations§
Source§impl SparseLdlt
impl SparseLdlt
Sourcepub fn factor(
n: usize,
col_ptr: &[usize],
row_idx: &[usize],
values: &[f64],
) -> Result<Self, LdltError>
pub fn factor( n: usize, col_ptr: &[usize], row_idx: &[usize], values: &[f64], ) -> Result<Self, LdltError>
Factor a symmetric n x n matrix supplied in CSC form.
col_ptrhas lengthn + 1; columnkoccupiescol_ptr[k]..col_ptr[k+1].row_idxandvaluesare parallel arrays of the nonzeros (any row order).
Only the upper triangle (entries with row ≤ col) is read; a fully symmetric matrix works too. No fill-reducing reordering is applied - permute the matrix first if you want one (RCM, AMD, nested dissection, …).
Sourcepub fn factor_reporting_collapse(
n: usize,
col_ptr: &[usize],
row_idx: &[usize],
values: &[f64],
) -> Result<(Self, Vec<usize>), LdltError>
pub fn factor_reporting_collapse( n: usize, col_ptr: &[usize], row_idx: &[usize], values: &[f64], ) -> Result<(Self, Vec<usize>), LdltError>
Like SparseLdlt::factor, but a near-zero pivot is RECORDED and the factorization
continues, instead of aborting at the first one.
This exists for RANK CHECKS, not for solves. A caller asking “which directions of this
matrix are null” needs the elimination to run to the end and name every column that
collapsed - a structure with fifteen mechanisms has fifteen of them, and stopping at the
first would report one. The returned factor is NOT fit to solve or to sign-count with:
every collapsed column’s pivot is rounding noise, exactly the value SparseLdlt::factor
refuses to return. Use the column list; discard d() for anything but structure.
An exact zero pivot still aborts, as it must: the elimination cannot proceed through it.
Sourcepub fn factor_perm(
n: usize,
col_ptr: &[usize],
row_idx: &[usize],
values: &[f64],
order: &[usize],
) -> Result<Self, LdltError>
pub fn factor_perm( n: usize, col_ptr: &[usize], row_idx: &[usize], values: &[f64], order: &[usize], ) -> Result<Self, LdltError>
Factor P A Pᵀ for a symmetric permutation P given as order, where
order[k] is the original index eliminated k-th (e.g. the output of amd).
The returned factorization solves A x = b DIRECTLY - the permutation is stored and
solve maps the right-hand side in and the solution back out, so callers that just
want answers use it exactly like SparseLdlt::factor. Fill-in drops because the
elimination order follows the ordering: on a random 2%-dense 1024 matrix the plain
factor carries ~9x the nonzeros of the AMD-ordered one.
Inertia is untouched by a symmetric permutation (Sylvester’s law: P A Pᵀ is a
congruence of A), so Sturm counts are identical with or without ordering.
§Errors
LdltError::InvalidInput if order is not a permutation of 0..n, plus
everything SparseLdlt::factor can return.
Sourcepub fn factor_shifted(
n: usize,
col_ptr: &[usize],
row_idx: &[usize],
values: &[f64],
) -> Result<Self, LdltError>
pub fn factor_shifted( n: usize, col_ptr: &[usize], row_idx: &[usize], values: &[f64], ) -> Result<Self, LdltError>
Like SparseLdlt::factor, but on a breakdown it retries with a positive diagonal
shift instead of giving up.
The unshifted factorization is tried first, so a well-conditioned matrix costs nothing
extra and comes back with SparseLdlt::shift == 0.0. On LdltError::ZeroPivot
or LdltError::NearZeroPivot the matrix is refactored as A + shift * I, starting
from the suggested shift and multiplying by 8 each attempt, at most 8 attempts; if none
succeeds the last error is returned.
THE RESULT IS AN EXACT FACTORIZATION OF A NEARBY MATRIX, NOT OF A. Its pivots are the
pivots of A + shift * I, so its inertia is that matrix’s inertia and a Sturm count
taken from it is a count at a sigma moved by shift. A solve against it is a solve of
the shifted system. Ignoring SparseLdlt::shift is a bug in the caller.
Sourcepub fn factor_perm_shifted(
n: usize,
col_ptr: &[usize],
row_idx: &[usize],
values: &[f64],
order: &[usize],
) -> Result<Self, LdltError>
pub fn factor_perm_shifted( n: usize, col_ptr: &[usize], row_idx: &[usize], values: &[f64], order: &[usize], ) -> Result<Self, LdltError>
SparseLdlt::factor_perm with the shifted-retry behaviour of
SparseLdlt::factor_shifted. The same warning applies: a non-zero
SparseLdlt::shift means this factored A + shift * I, not A.
Sourcepub fn shift(&self) -> f64
pub fn shift(&self) -> f64
The diagonal shift actually applied. 0.0 for SparseLdlt::factor /
SparseLdlt::factor_perm, which never shift. Non-zero means this is an exact
factorization of A + shift * I, NOT of A: its inertia is the inertia of the shifted
matrix, so a Sturm count taken from it is a count for the caller’s sigma moved by this
much, and the caller must correct for it.
Sourcepub fn d(&self) -> &[f64]
pub fn d(&self) -> &[f64]
The signed diagonal D. The count of negative entries is the matrix inertia
(number of negative eigenvalues), e.g. for a Sturm eigenvalue count.
Sourcepub fn flops(&self) -> u64
pub fn flops(&self) -> u64
Floating-point operation count of the factorization: for each column of L with
c stored entries, c*c + 3*c (the column-update arithmetic). Deterministic, so
two factorizations of the same sparsity pattern report identical counts - callers
(e.g. the supernodal equivalence gates in FEM Studio) assert on exactly that.
Sourcepub fn solve(&self, b: &[f64]) -> Result<Vec<f64>, LdltError>
pub fn solve(&self, b: &[f64]) -> Result<Vec<f64>, LdltError>
Solve A x = b for a single right-hand side, returning x.
Works for both SparseLdlt::factor and SparseLdlt::factor_perm - the stored
elimination order is applied to the right-hand side and inverted on the solution,
so the caller never sees the permutation.
§Errors
Returns LdltError::SizeMismatch if b.len() != self.dim().