diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index dbbcb23..ae82633 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -265,6 +265,12 @@ pub(crate) fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) - "the stiffness is singular at {at}: the model can move without straining (missing or insufficient boundary conditions, a part not connected to the rest, or a reference point direction nothing holds)" )) } + FactorError::Singular { column } => { + let at = reduction.eq_dof.get(column).map_or_else(|| "an unknown degree of freedom".to_string(), |d| mesh.describe_dof(*d)); + Error::Singular(format!( + "the matrix is singular at {at}: no pivot of its supernode is more than round-off (a mechanism, or a shift or load at an eigenvalue)" + )) + } } } diff --git a/crates/solver/src/cholesky.rs b/crates/solver/src/cholesky.rs index 54334fa..6b2bc73 100644 --- a/crates/solver/src/cholesky.rs +++ b/crates/solver/src/cholesky.rs @@ -211,6 +211,10 @@ fn postorder(parent: &[usize]) -> Vec { pub enum FactorError { /// The pivot of this column is not positive (relative to the diagonal). NotPositiveDefinite { column: usize }, + /// The symmetric indefinite factorization ([`crate::ldlt`]) found no pivot in this column's + /// supernode that is not negligible next to the column's entries: the matrix is singular to + /// working precision there. + Singular { column: usize }, /// The factor would hold more entries than allowed. TooLarge { entries: usize, limit: usize }, /// The memory for the factor could not be had. @@ -219,6 +223,25 @@ pub enum FactorError { Internal(String), } +impl FactorError { + /// The column a numeric failure names. + fn column(&self) -> Option { + match self { + FactorError::NotPositiveDefinite { column } | FactorError::Singular { column } => Some(*column), + _ => None, + } + } + + /// The same failure at another column. + fn at(self, column: usize) -> FactorError { + match self { + FactorError::NotPositiveDefinite { .. } => FactorError::NotPositiveDefinite { column }, + FactorError::Singular { .. } => FactorError::Singular { column }, + e => e, + } + } +} + /// One supernode's update to a later one: supernode `k`'s rows from position `p1` on reach the /// later one, whose columns they meet at positions `p1..p2`. #[derive(Debug, Clone, Copy)] @@ -230,7 +253,7 @@ struct Update { /// The supernodal structure (postorder numbering). #[derive(Debug, Clone)] -struct Structure { +pub(crate) struct Structure { /// First column of each supernode, and `n` at the end. first: Vec, /// Row structure of each supernode (its own columns first, then the rows below), sorted. @@ -239,64 +262,108 @@ struct Structure { } impl Structure { - fn cols(&self, s: usize) -> std::ops::Range { + pub(crate) fn cols(&self, s: usize) -> std::ops::Range { self.first.get(s).copied().unwrap_or(0)..self.first.get(s + 1).copied().unwrap_or(0) } - fn rows(&self, s: usize) -> &[usize] { + pub(crate) fn rows(&self, s: usize) -> &[usize] { let (a, b) = (self.rptr.get(s).copied().unwrap_or(0), self.rptr.get(s + 1).copied().unwrap_or(0)); self.rows.get(a..b).unwrap_or_default() } /// Entries of the supernode's lower trapezoid. - fn size(&self, s: usize) -> usize { + pub(crate) fn size(&self, s: usize) -> usize { trapezoid(self.rows(s).len(), self.cols(s).len()) } } /// Entries of an m-row, w-column lower trapezoid (column c holds rows c..m). -fn trapezoid(m: usize, w: usize) -> usize { +pub(crate) fn trapezoid(m: usize, w: usize) -> usize { (0..w).map(|c| m.saturating_sub(c)).sum() } /// Where column c of an m-row trapezoid starts. -fn offset(m: usize, c: usize) -> usize { +pub(crate) fn offset(m: usize, c: usize) -> usize { c * m - c * c.saturating_sub(1) / 2 } /// A dense lower trapezoid of m rows: one supernode's values, column by column, column c /// holding rows c..m of the supernode's row structure. -struct Panel<'a> { - m: usize, - data: &'a [f64], +pub(crate) struct Panel<'a> { + pub(crate) m: usize, + pub(crate) data: &'a [f64], } impl Panel<'_> { /// L(r, c) for r ≥ c. - fn at(&self, r: usize, c: usize) -> f64 { + pub(crate) fn at(&self, r: usize, c: usize) -> f64 { self.data.get(offset(self.m, c) + r - c).copied().unwrap_or(0.0) } /// Rows r..r+len of column c (r ≥ c). - fn column(&self, c: usize, r: usize, len: usize) -> &[f64] { + pub(crate) fn column(&self, c: usize, r: usize, len: usize) -> &[f64] { let start = offset(self.m, c) + r - c; self.data.get(start..start + len).unwrap_or_default() } } -/// S(i, c) = Σ_k L(p + i, k) L(p + c, k) over the columns k < `kw` of a trapezoid, for rows i in -/// 0..mu and columns c in 0..wu (meaningful for i ≥ c), each sum in column order from zero. -/// Returned by row tiles: tile t holds rows t·TILE.. as `tile[c * TILE + i]`. -fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec> { - // Rows p..p+wu of the kw columns, transposed so that column c's coefficients are contiguous. +/// The pivots of a symmetric indefinite factor's block diagonal D, per column of a supernode +/// ([`crate::ldlt`]): a 1 × 1 pivot, or the first or second column of a 2 × 2 one. D's entries +/// sit on the trapezoid's diagonal, and a 2 × 2 pivot's off-diagonal entry below its first +/// column's diagonal (where L's unit block has a zero). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Pivot { + One, + First, + Second, +} + +/// `row` (entries for the columns k0.. of a trapezoid) ← `row` D, with D the block diagonal the +/// trapezoid holds for those columns and `pivots` its pivots from column k0. +pub(crate) fn times_d(row: &mut [f64], l: &Panel, k0: usize, pivots: &[Pivot]) { + let mut k = 0; + while k < row.len() { + let d = |r: usize, c: usize| l.at(k0 + r, k0 + c); + if pivots.get(k) == Some(&Pivot::First) { + let (a, b) = (row.get(k).copied().unwrap_or(0.0), row.get(k + 1).copied().unwrap_or(0.0)); + let (d11, d21, d22) = (d(k, k), d(k + 1, k), d(k + 1, k + 1)); + if let Some(x) = row.get_mut(k) { + *x = a * d11 + b * d21; + } + if let Some(x) = row.get_mut(k + 1) { + *x = a * d21 + b * d22; + } + k += 2; + } else { + if let Some(x) = row.get_mut(k) { + *x *= d(k, k); + } + k += 1; + } + } +} + +/// S(i, c) = Σ_k L(p + i, k) Y(p + c, k) over the columns k in `cols` of a trapezoid, for rows i +/// in 0..mu and columns c in 0..wu (meaningful for i ≥ c), each sum in column order. Y is L +/// (Cholesky), or L D with D the block diagonal the trapezoid holds when `pivots` (of the +/// columns `cols`) are given. Returned by row tiles: tile t holds rows t·TILE.. as +/// `tile[c * TILE + i]`. +pub(crate) fn products(l: &Panel, cols: std::ops::Range, p: usize, mu: usize, wu: usize, pivots: Option<&[Pivot]>) -> Vec> { + let (k0, kw) = (cols.start, cols.len()); + // Rows p..p+wu of the columns, transposed so that column c's coefficients are contiguous. let mut yt = vec![0.0; wu * kw]; for k in 0..kw { - for (c, y) in l.column(k, p, wu).iter().enumerate() { + for (c, y) in l.column(k0 + k, p, wu).iter().enumerate() { if let Some(slot) = yt.get_mut(c * kw + k) { *slot = *y; } } } + if let Some(pivots) = pivots { + for row in yt.chunks_mut(kw.max(1)) { + times_d(row, l, k0, pivots); + } + } let tiles = mu.div_ceil(TILE); let tile = |t: usize| -> Vec { let i0 = t * TILE; @@ -305,7 +372,7 @@ fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec Vec], mu: usize, wu: usize, rel: impl Fn(usize) -> usize) { +pub(crate) fn subtract(own: &mut [f64], m: usize, tiles: &[Vec], mu: usize, wu: usize, rel: impl Fn(usize) -> usize) { for (t, tile) in tiles.iter().enumerate() { let i0 = t * TILE; let h = TILE.min(mu - i0); @@ -379,11 +446,13 @@ fn solve_below(l: &Panel, c0: usize, c1: usize) -> Vec> { } /// The matrix in postorder by lower columns: column j holds A(i, j) for i ≥ j, rows sorted. -struct Lower { +pub(crate) struct Lower { ptr: Vec, idx: Vec, val: Vec, diag: Vec, + /// The largest magnitude in each row (= column) of the whole symmetric matrix. + pub(crate) scale: Vec, } impl Lower { @@ -407,6 +476,7 @@ impl Lower { let mut pairs = vec![(0usize, 0.0f64); ptr.last().copied().unwrap_or(0)]; let mut next = ptr.clone(); let mut diag = vec![0.0; n]; + let mut scale = vec![0.0f64; n]; for k in 0..n { let (rows, vals) = a.column(k); for (&i, &v) in rows.iter().zip(vals) { @@ -416,6 +486,11 @@ impl Lower { { *d = v; } + for x in [r, c] { + if let Some(s) = scale.get_mut(x) { + *s = s.max(v.abs()); + } + } if let Some(slot) = next.get_mut(c) { if let Some(p) = pairs.get_mut(*slot) { *p = (r, v); @@ -431,7 +506,7 @@ impl Lower { } } let (idx, val) = pairs.into_iter().unzip(); - Lower { ptr, idx, val, diag } + Lower { ptr, idx, val, diag, scale } } fn column(&self, j: usize) -> (&[usize], &[f64]) { @@ -442,17 +517,38 @@ impl Lower { /// What a supernode's factorization reads: the matrix, the structure and the updates each /// supernode receives. -struct Context<'a> { - a: &'a Lower, - s: &'a Structure, +pub(crate) struct Context<'a> { + pub(crate) a: &'a Lower, + pub(crate) s: &'a Structure, updates: &'a [Vec], } +/// What the numeric factorization of one kind keeps per supernode. +pub(crate) trait Supernode: Send + Sync { + /// The supernode's values: its lower trapezoid by columns. + fn values(&self) -> &[f64]; + /// The pivots of D (symmetric indefinite factors), `None` for Cholesky. + fn pivots(&self) -> Option<&[Pivot]>; +} + +impl Supernode for Vec { + fn values(&self) -> &[f64] { + self + } + fn pivots(&self) -> Option<&[Pivot]> { + None + } +} + +fn internal(what: &str) -> FactorError { + FactorError::Internal(format!("supernodal factorization: {what}")) +} + impl Context<'_> { - /// Factor supernode `j` into `own`; its descendants' values are `done[k - base]`. Fails with - /// the first column (postorder) whose pivot is not positive. - fn supernode(&self, j: usize, done: &[Vec], base: usize, own: &mut [f64]) -> std::result::Result<(), FactorError> { - let internal = |what: &str| FactorError::Internal(format!("supernodal factorization: {what}")); + /// Supernode `j` assembled into `own`: its columns of A minus the update L_K Y_Kᵀ of every + /// descendant K that reaches it (Y = L, or L D), in supernode order; the descendants are + /// `done[k - base]`. + pub(crate) fn assemble(&self, j: usize, done: &[T], base: usize, own: &mut [f64]) -> std::result::Result<(), FactorError> { let rows = self.s.rows(j); let cols = self.s.cols(j); let (m, w, f) = (rows.len(), cols.len(), cols.start); @@ -476,7 +572,8 @@ impl Context<'_> { // The descendants' updates, in supernode order. let mut relative = Vec::new(); for u in self.updates.get(j).map(Vec::as_slice).unwrap_or_default() { - let lk = done.get(u.k.wrapping_sub(base)).ok_or_else(|| internal("an update from outside the subtree"))?; + let done_k = done.get(u.k.wrapping_sub(base)).ok_or_else(|| internal("an update from outside the subtree"))?; + let lk = done_k.values(); let (kr, kw) = (self.s.rows(u.k), self.s.cols(u.k).len()); let below = kr.get(u.p1..).unwrap_or_default(); let (mu, wu) = (below.len(), u.p2 - u.p1); @@ -492,16 +589,26 @@ impl Context<'_> { } relative.push(at); } - let tiles = products(&Panel { m: kr.len(), data: lk }, kw, u.p1, mu, wu); + let tiles = products(&Panel { m: kr.len(), data: lk }, 0..kw, u.p1, mu, wu, done_k.pivots()); subtract(own, m, &tiles, mu, wu, |i| relative.get(i).copied().unwrap_or(0)); } + Ok(()) + } + + /// Factor supernode `j` into `own`; its descendants' values are `done[k - base]`. Fails with + /// the first column (postorder) whose pivot is not positive. + fn supernode(&self, j: usize, done: &[Vec], base: usize, own: &mut [f64]) -> std::result::Result<(), FactorError> { + self.assemble(j, done, base, own)?; + let rows = self.s.rows(j); + let cols = self.s.cols(j); + let (m, w, f) = (rows.len(), cols.len(), cols.start); // Dense factorization by panels. for c0 in (0..w).step_by(PANEL) { let c1 = (c0 + PANEL).min(w); let pw = c1 - c0; if c0 > 0 { // The panel minus its products with the columns before it. - let tiles = products(&Panel { m, data: own }, c0, c0, m - c0, pw); + let tiles = products(&Panel { m, data: own }, 0..c0, c0, m - c0, pw, None); subtract(own, m, &tiles, m - c0, pw, |i| c0 + i); } // The diagonal block, column by column. @@ -549,20 +656,28 @@ impl Context<'_> { } } -/// A Cholesky factor: L in supernodes, for the matrix with its columns in postorder. -#[derive(Debug, Clone)] -pub struct Factor { - pub n: usize, +/// The symbolic factorization shared by the Cholesky and the symmetric indefinite (LDLᵀ, +/// [`crate::ldlt`]) factors: the postorder, the supernodes and their row structures, the updates +/// each supernode receives, the matrix in that order, and the subtree tasks of the numeric +/// phase. +pub(crate) struct Symbolic { + pub(crate) n: usize, /// `perm[k]`: the original column that is column k of the factor. - perm: Vec, - structure: Structure, - /// The values of each supernode (a lower trapezoid by columns). - values: Vec>, - entries: usize, + pub(crate) perm: Vec, + pub(crate) structure: Structure, + a: Lower, + updates: Vec>, + /// Subtrees factored as parallel tasks, and the supernodes above them, in order. + tasks: Vec>, + top: Vec, + /// Entries of L. + pub(crate) entries: usize, } -impl Factor { - pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result { +impl Symbolic { + /// The structure of the factor of `a`, refused before any numeric memory is taken when the + /// factor would hold more than `max_entries` entries. + pub(crate) fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result { let n = a.n; let internal = |what: &str| FactorError::Internal(format!("supernodal structure: {what}")); let (parent0, counts0) = symbolic(a); @@ -682,16 +797,6 @@ impl Factor { } } - // Memory for the values, taken without aborting. - let mut values: Vec> = Vec::new(); - values.try_reserve_exact(ns).map_err(|_| FactorError::OutOfMemory { entries: total })?; - for s in 0..ns { - let mut v = Vec::new(); - v.try_reserve_exact(structure.size(s)).map_err(|_| FactorError::OutOfMemory { entries: total })?; - v.resize(structure.size(s), 0.0); - values.push(v); - } - // Tasks: subtrees (contiguous ranges of supernodes in postorder) of at most a share of the // work each; the supernodes above them are factored afterwards, in order. let mut subtree: Vec = (0..ns) @@ -737,12 +842,49 @@ impl Factor { } } - let context = Context { a: &a, s: &structure, updates: &updates }; + Ok(Symbolic { n, perm, structure, a, updates, tasks, top, entries: total }) + } + + /// Supernodes. + pub(crate) fn supernodes(&self) -> usize { + self.structure.first.len().saturating_sub(1) + } + + /// The original column of factor column `column`. + pub(crate) fn original(&self, column: usize) -> usize { + self.perm.get(column).copied().unwrap_or(column) + } + + /// Zeroed memory for every supernode's trapezoid, taken without aborting. + pub(crate) fn allocate(&self) -> std::result::Result>, FactorError> { + let ns = self.supernodes(); + let mut values: Vec> = Vec::new(); + values.try_reserve_exact(ns).map_err(|_| FactorError::OutOfMemory { entries: self.entries })?; + for s in 0..ns { + let mut v = Vec::new(); + v.try_reserve_exact(self.structure.size(s)).map_err(|_| FactorError::OutOfMemory { entries: self.entries })?; + v.resize(self.structure.size(s), 0.0); + values.push(v); + } + Ok(values) + } + + /// Run `supernode(context, j, done, base, own)` over every supernode, children before + /// parents: the subtree tasks in parallel, each through its supernodes in order, then the + /// supernodes above them in order. A numeric failure is the first failing column in + /// postorder (the one a pass in order would meet first, whichever thread meets it), named by + /// its original column. + pub(crate) fn factor( + &self, + values: &mut [T], + supernode: impl Fn(&Context, usize, &[T], usize, &mut T) -> std::result::Result<(), FactorError> + Sync, + ) -> std::result::Result<(), FactorError> { + let context = Context { a: &self.a, s: &self.structure, updates: &self.updates }; // The tasks run in parallel, each through its supernodes in order. - let mut pieces: Vec<(std::ops::Range, &mut [Vec], Option)> = Vec::with_capacity(tasks.len()); - let mut rest: &mut [Vec] = &mut values; + let mut pieces: Vec<(std::ops::Range, &mut [T], Option)> = Vec::with_capacity(self.tasks.len()); + let mut rest: &mut [T] = values; let mut at = 0; - for range in &tasks { + for range in &self.tasks { let (_, tail) = std::mem::take(&mut rest).split_at_mut(range.start - at); let (piece, tail) = tail.split_at_mut(range.len()); pieces.push((range.clone(), piece, None)); @@ -750,14 +892,14 @@ impl Factor { at = range.end; } let broken = std::sync::Mutex::new(None::); - parallel::for_each_mut(&mut pieces, |_, (range, piece, singular)| { + parallel::for_each_mut(&mut pieces, |_, (range, piece, failed)| { for j in range.clone() { let (done, own) = piece.split_at_mut(j - range.start); let Some(own) = own.first_mut() else { break }; - match context.supernode(j, done, range.start, own) { + match supernode(&context, j, done, range.start, own) { Ok(()) => {} - Err(FactorError::NotPositiveDefinite { column }) => { - *singular = Some(column); + Err(e) if e.column().is_some() => { + *failed = Some(e); break; } Err(e) => { @@ -772,27 +914,50 @@ impl Factor { if let Some(e) = broken.into_inner().unwrap_or_else(std::sync::PoisonError::into_inner) { return Err(e); } - // The first singular column in postorder is the one a pass in order would meet first: + // The first failing column in postorder is the one a pass in order would meet first: // the supernodes above the tasks run up to it. - let earliest = pieces.iter().filter_map(|p| p.2).min(); + let earliest = pieces.iter_mut().filter_map(|p| p.2.take()).min_by_key(|e| e.column()); drop(pieces); - let original = |column: usize| perm.get(column).copied().unwrap_or(column); - for &j in &top { - if earliest.is_some_and(|c| structure.cols(j).start > c) { + let named = |e: FactorError| match e.column() { + Some(c) => e.at(self.original(c)), + None => e, + }; + for &j in &self.top { + if earliest.as_ref().and_then(FactorError::column).is_some_and(|c| self.structure.cols(j).start > c) { break; } let (done, own) = values.split_at_mut(j); let Some(own) = own.first_mut() else { break }; - match context.supernode(j, done, 0, own) { - Ok(()) => {} - Err(FactorError::NotPositiveDefinite { column }) => return Err(FactorError::NotPositiveDefinite { column: original(column) }), - Err(e) => return Err(e), - } + supernode(&context, j, done, 0, own).map_err(named)?; } - if let Some(column) = earliest { - return Err(FactorError::NotPositiveDefinite { column: original(column) }); + match earliest { + Some(e) => Err(named(e)), + None => Ok(()), } - Ok(Factor { n, perm, structure, values, entries: total }) + } +} + +/// A Cholesky factor: L in supernodes, for the matrix with its columns in postorder. +#[derive(Debug, Clone)] +pub struct Factor { + pub n: usize, + /// `perm[k]`: the original column that is column k of the factor. + perm: Vec, + structure: Structure, + /// The values of each supernode (a lower trapezoid by columns). + values: Vec>, + entries: usize, +} + +impl Factor { + pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result { + Factor::from_symbolic(&Symbolic::new(a, max_entries)?) + } + + pub(crate) fn from_symbolic(symbolic: &Symbolic) -> std::result::Result { + let mut values = symbolic.allocate()?; + symbolic.factor(&mut values, |context, j, done, base, own| context.supernode(j, done, base, own))?; + Ok(Factor { n: symbolic.n, perm: symbolic.perm.clone(), structure: symbolic.structure.clone(), values, entries: symbolic.entries }) } /// Entries in L. diff --git a/crates/solver/src/heat.rs b/crates/solver/src/heat.rs index 5df9298..9b676fa 100644 --- a/crates/solver/src/heat.rs +++ b/crates/solver/src/heat.rs @@ -310,7 +310,7 @@ fn factor(a: &Csr, limits: &Limits, mesh: &HeatMesh, eqs: &Equations) -> Result< } FactorError::OutOfMemory { entries } => Error::Limit(format!("not enough memory for the factorization ({entries} entries)")), FactorError::Internal(what) => Error::Internal(what), - FactorError::NotPositiveDefinite { column } => { + FactorError::NotPositiveDefinite { column } | FactorError::Singular { column } => { let label = eqs.node.get(column).and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0); Error::Singular(format!( "the temperature is undetermined at node {label}: no temperature boundary condition or convection reaches its part of the model" diff --git a/crates/solver/src/ldlt.rs b/crates/solver/src/ldlt.rs new file mode 100644 index 0000000..1c99f4b --- /dev/null +++ b/crates/solver/src/ldlt.rs @@ -0,0 +1,886 @@ +//! Sparse symmetric indefinite factorization P A Pᵀ = L D Lᵀ: L unit lower triangular, D block +//! diagonal with 1 × 1 and 2 × 2 pivots; supernodal, left-looking between supernodes, +//! multithreaded, with a fixed operation order, on the symbolic factorization of the Cholesky +//! factor ([`crate::cholesky`]). +//! +//! **Pivoting** (J. R. Bunch, L. Kaufman, *Some stable methods for calculating inertia and solving +//! symmetric linear systems*, Math. Comp. 31(137):163–179, 1977), in its bounded form (C. +//! Ashcraft, R. G. Grimes, J. G. Lewis, *Accurate symmetric indefinite linear equation solvers*, +//! SIAM J. Matrix Anal. Appl. 20(2):513–561, 1998, the "bounded Bunch–Kaufman" or rook search): +//! with α = (1 + √17)/8 and ω_i the largest off-diagonal magnitude of column i, column k is a +//! 1 × 1 pivot when |a_kk| ≥ α ω_k; otherwise the search walks from column i to the row r of its +//! largest entry: r is a 1 × 1 pivot when |a_rr| ≥ α ω_r, (i, r) a 2 × 2 pivot when ω_r = ω_i +//! (the entry is the largest of both columns), else the search continues from r. Every entry of +//! L is then bounded by 1/(1 − α) ≈ 2.78 for 1 × 1 pivots (and a like constant for 2 × 2 ones), +//! which Bunch–Kaufman's single look does not guarantee; ω grows strictly, so the search ends. +//! A 2 × 2 pivot has |a_ii a_rr| < α² a_ir², a negative determinant: one positive and one +//! negative eigenvalue. +//! +//! **Where pivots come from: within supernodes** (the supernode Bunch–Kaufman of O. Schenk, K. +//! Gärtner, *On fast factorization pivoting methods for sparse symmetric indefinite systems*, +//! Electron. Trans. Numer. Anal. 23:158–179, 2006). The columns of a supernode are fully summed +//! once its descendants' updates are in (left-looking, as the Cholesky factor), so any of them +//! can be the next pivot: a symmetric interchange inside the supernode leaves the row structure +//! of every supernode unchanged, so the static symbolic factorization, the ordering, the +//! supernodal tasks and the fill are the Cholesky factor's exactly. The pivot search looks at +//! the supernode's own rows; entries in the rows below (later supernodes) are not candidates. +//! Multifrontal codes delay a column whose pivot fails into the parent (I. S. Duff, J. K. Reid, +//! *The multifrontal solution of indefinite sparse symmetric linear equations*, ACM TOMS +//! 9(3):302–325, 1983; J. W. H. Liu, *The multifrontal method for sparse matrix solution: theory +//! and practice*, SIAM Review 34(1):82–109, 1992), which changes the structure as it goes; this +//! factor keeps the static structure instead (deterministic, threads as Cholesky, the same +//! memory bound checked up front), and: +//! - a column with no acceptable pivot in its supernode — every remaining candidate negligible, +//! at most `PIVOT_TOLERANCE` times the largest entry of the matrix's row — is reported as +//! singular, naming the column; a nonsingular saddle point whose zero diagonal couples only to +//! rows of a later supernode would also stop there (none of FieldCraft's systems has one: the +//! hybrid elements' pressures are condensed per element); +//! - L's entries in the rows below the supernode are not bounded by the search; the stiffness, +//! shifted and mixed tangents solved here have sizeable diagonals, and +//! [`Symmetric::solve_refined`] adds iterative refinement against A where accuracy matters. +//! +//! **Numeric.** A supernode starts as its columns of A minus L_K D_K L_Kᵀ of every descendant K +//! that reaches it, in supernode order (the L D rows formed on the fly from L and D, so nothing +//! more is stored than Cholesky's L plus the pivot order). It is then factored right-looking by +//! panels of [`PANEL`] columns: inside a panel each column is brought up to date with the panel's +//! columns before it when the search looks at it (the partial factorization of LAPACK's +//! `dsytrf`, Anderson et al., *LAPACK Users' Guide*, SIAM 1999), and after the panel the rest of +//! the supernode receives its update by the tiled, parallel kernel of the Cholesky factor. A +//! 2 × 2 pivot may end a panel one column late. +//! +//! **Inertia** (Sylvester's law): the counts of positive and negative eigenvalues of A are those +//! of D — a 1 × 1 pivot by its sign, a 2 × 2 one by its eigenvalues — so a factorization of +//! K − σ M counts the eigenvalues below σ, and a tangent's negative pivots count its unstable +//! directions. A pivot that would be zero stops the factorization (above), so the count of zero +//! eigenvalues is never reported as anything but a located error. +//! +//! **Threads and determinism.** Every pivot decision and every entry is computed by one task in +//! an order fixed by the matrix, ties in the search broken by the lowest index, so the factor, +//! its pivots and its inertia are the same bit for bit on any number of threads. + +use serde::Serialize; + +use crate::cholesky::{ + Context, Factor, FactorError, PANEL, PIVOT_TOLERANCE, Panel, Pivot, Structure, Supernode, Symbolic, UpperColumns, offset, products, subtract, + times_d, +}; +use crate::sparse::Csr; +use crate::{Error, Result}; + +/// (1 + √17)/8: the growth bound of Bunch and Kaufman. +const ALPHA: f64 = 0.640_388_203_202_207_6; + +/// Most steps of one rook search (each strictly raises ω, so the cap only bounds a pathological +/// sequence of near ties); at the cap the last two columns are the 2 × 2 pivot. +const ROOK_STEPS: usize = 64; + +/// The counts of positive and negative eigenvalues of a factored symmetric matrix. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize)] +pub struct Inertia { + pub positive: usize, + pub negative: usize, +} + +/// One supernode of an LDLᵀ factor. +#[derive(Debug, Clone, Default)] +struct Block { + /// L below the pivots and D on them, the lower trapezoid by columns in pivot order. + values: Vec, + pivots: Vec, + /// `order[c]`: the supernode's own column (0..w, in the symbolic order) at pivot position c. + order: Vec, + inertia: Inertia, + two_by_two: usize, +} + +impl Supernode for Block { + fn values(&self) -> &[f64] { + &self.values + } + fn pivots(&self) -> Option<&[Pivot]> { + Some(&self.pivots) + } +} + +/// The pivot a search chose, at supernode-local positions not yet moved to k. +enum Choice { + One(usize), + Two(usize, usize), +} + +/// Symmetric interchange of local columns p < q (both among the supernode's own columns) in an +/// m-row lower trapezoid: the rows of the columns before p, the diagonals, the strip between them +/// and the rows below q. +fn interchange(v: &mut [f64], m: usize, p: usize, q: usize) { + let at = |r: usize, c: usize| offset(m, c) + r - c; + for c in 0..p { + v.swap(at(p, c), at(q, c)); + } + v.swap(at(p, p), at(q, q)); + for c in p + 1..q { + v.swap(at(c, p), at(q, c)); + } + for r in q + 1..m { + v.swap(at(r, p), at(r, q)); + } +} + +/// The largest magnitude of `v[i]` over `range` except `skip`, with the lowest such index; NaN +/// when an entry is not a number. +fn largest(v: &[f64], range: std::ops::Range, skip: Option) -> (usize, f64) { + let mut best = (range.start, 0.0f64); + for i in range { + if Some(i) == skip { + continue; + } + let x = v.get(i).copied().unwrap_or(0.0).abs(); + if x.is_nan() { + return (i, f64::NAN); + } + if x > best.1 { + best = (i, x); + } + } + best +} + +/// The work buffers of one supernode's panels. +struct Scratch { + /// The current values of the column being searched and of a candidate, by local row. + a: Vec, + b: Vec, + /// A row of L D over the panel's factored columns. + w: Vec, +} + +impl Block { + /// Column `c` (local, ≥ k) of the part still to factor, rows k..m, into `out`: the trapezoid's + /// entries less the products with the panel's factored columns c0..k. + fn current(&self, m: usize, c0: usize, k: usize, c: usize, out: &mut [f64], w: &mut Vec) { + let t = Panel { m, data: &self.values }; + for r in k..m { + if let Some(o) = out.get_mut(r) { + *o = if r >= c { t.at(r, c) } else { t.at(c, r) }; + } + } + w.clear(); + w.extend((c0..k).map(|q| t.at(c, q))); + times_d(w, &t, c0, self.pivots.get(c0..k).unwrap_or_default()); + for (q, wq) in (c0..k).zip(w.iter()) { + for (o, l) in out.get_mut(k..m).unwrap_or_default().iter_mut().zip(t.column(q, k, m - k)) { + *o -= l * wq; + } + } + } + + /// Swap local columns p and q (pivot positions ≥ k) in the trapezoid, the order and the + /// current columns. + fn swap(&mut self, m: usize, p: usize, q: usize, scratch: &mut Scratch) { + if p == q { + return; + } + let (p, q) = (p.min(q), p.max(q)); + interchange(&mut self.values, m, p, q); + self.order.swap(p, q); + scratch.a.swap(p, q); + scratch.b.swap(p, q); + } + + /// Store a 1 × 1 pivot at k from the current column `col`. + fn one(&mut self, m: usize, k: usize, col: &[f64]) { + let d = col.get(k).copied().unwrap_or(0.0); + let start = offset(m, k); + if let Some(x) = self.values.get_mut(start) { + *x = d; + } + for r in k + 1..m { + if let Some(x) = self.values.get_mut(start + r - k) { + *x = col.get(r).copied().unwrap_or(0.0) / d; + } + } + if let Some(p) = self.pivots.get_mut(k) { + *p = Pivot::One; + } + if d > 0.0 { + self.inertia.positive += 1; + } else { + self.inertia.negative += 1; + } + } + + /// Store a 2 × 2 pivot at k, k + 1 from the current columns `a` (k) and `b` (k + 1). + fn two(&mut self, m: usize, k: usize, a: &[f64], b: &[f64]) { + let get = |v: &[f64], r: usize| v.get(r).copied().unwrap_or(0.0); + let (d11, d21, d22) = (get(a, k), get(a, k + 1), get(b, k + 1)); + let det = d11 * d22 - d21 * d21; + let (s0, s1) = (offset(m, k), offset(m, k + 1)); + for (at, v) in [(s0, d11), (s0 + 1, d21), (s1, d22)] { + if let Some(x) = self.values.get_mut(at) { + *x = v; + } + } + for r in k + 2..m { + let (x, y) = (get(a, r), get(b, r)); + if let Some(l) = self.values.get_mut(s0 + r - k) { + *l = (x * d22 - y * d21) / det; + } + if let Some(l) = self.values.get_mut(s1 + r - k - 1) { + *l = (y * d11 - x * d21) / det; + } + } + if let Some(p) = self.pivots.get_mut(k) { + *p = Pivot::First; + } + if let Some(p) = self.pivots.get_mut(k + 1) { + *p = Pivot::Second; + } + // The eigenvalues' signs: one of each when det < 0, else both the sign of the trace. + if det < 0.0 { + self.inertia.positive += 1; + self.inertia.negative += 1; + } else if d11 + d22 > 0.0 { + self.inertia.positive += 2; + } else { + self.inertia.negative += 2; + } + self.two_by_two += 1; + } + + /// Factor the panel from column c0 of a supernode with m rows and w columns: pivots chosen + /// among the columns c0..w, each column brought up to date with the panel's columns before + /// it. Returns where the panel ended. `scale(local)` is the largest entry of the matrix in + /// that column; `column(local)` its postorder column. + fn panel( + &mut self, + m: usize, + w: usize, + c0: usize, + scale: &dyn Fn(usize) -> f64, + column: &dyn Fn(usize) -> usize, + scratch: &mut Scratch, + ) -> std::result::Result { + let end = (c0 + PANEL).min(w); + let mut k = c0; + while k < end { + let singular = |b: &Block, at: usize| FactorError::Singular { column: column(b.order.get(at).copied().unwrap_or(at)) }; + let negligible = |b: &Block, at: usize, x: f64| !(x > PIVOT_TOLERANCE * scale(b.order.get(at).copied().unwrap_or(at))); + let mut a = std::mem::take(&mut scratch.a); + self.current(m, c0, k, k, &mut a, &mut scratch.w); + scratch.a = a; + let akk = scratch.a.get(k).copied().unwrap_or(0.0).abs(); + let (r0, omega_k) = largest(&scratch.a, k + 1..w, None); + if !akk.is_finite() || !omega_k.is_finite() { + return Err(singular(self, k)); + } + let choice = if akk >= ALPHA * omega_k { + Choice::One(k) + } else { + // The rook search: `a` holds column i, `b` the candidate r. + let (mut i, mut omega_i, mut r) = (k, omega_k, r0); + let mut step = 0; + loop { + let mut b = std::mem::take(&mut scratch.b); + self.current(m, c0, k, r, &mut b, &mut scratch.w); + scratch.b = b; + let (s, omega_r) = largest(&scratch.b, k..w, Some(r)); + let arr = scratch.b.get(r).copied().unwrap_or(0.0).abs(); + if !arr.is_finite() || !omega_r.is_finite() { + return Err(singular(self, k)); + } + step += 1; + if arr >= ALPHA * omega_r { + std::mem::swap(&mut scratch.a, &mut scratch.b); + break Choice::One(r); + } + if omega_r <= omega_i || step >= ROOK_STEPS { + break Choice::Two(i, r); + } + std::mem::swap(&mut scratch.a, &mut scratch.b); + (i, omega_i, r) = (r, omega_r, s); + } + }; + match choice { + Choice::One(p) => { + self.swap(m, k, p, scratch); + if negligible(self, k, scratch.a.get(k).copied().unwrap_or(0.0).abs()) { + return Err(singular(self, k)); + } + let a = std::mem::take(&mut scratch.a); + self.one(m, k, &a); + scratch.a = a; + k += 1; + } + Choice::Two(i, r) => { + // i to k, then r to k + 1 (r was at k if i took its place). + self.swap(m, k, i, scratch); + let r = if r == k { i } else { r }; + self.swap(m, k + 1, r, scratch); + let offdiagonal = scratch.a.get(k + 1).copied().unwrap_or(0.0).abs(); + let size = scale(self.order.get(k).copied().unwrap_or(k)).max(scale(self.order.get(k + 1).copied().unwrap_or(k + 1))); + if !(offdiagonal > PIVOT_TOLERANCE * size) { + return Err(singular(self, k)); + } + let (a, b) = (std::mem::take(&mut scratch.a), std::mem::take(&mut scratch.b)); + self.two(m, k, &a, &b); + (scratch.a, scratch.b) = (a, b); + k += 2; + } + } + } + Ok(k) + } +} + +/// Factor supernode `j` into `own`; its descendants are `done[k - base]`. +fn supernode(context: &Context, j: usize, done: &[Block], base: usize, own: &mut Block) -> std::result::Result<(), FactorError> { + context.assemble(j, done, base, &mut own.values)?; + let (m, cols) = (context.s.rows(j).len(), context.s.cols(j)); + let (w, f) = (cols.len(), cols.start); + own.pivots = vec![Pivot::One; w]; + own.order = (0..w).collect(); + own.inertia = Inertia::default(); + own.two_by_two = 0; + let scale = |local: usize| context.a.scale.get(f + local).copied().unwrap_or(0.0); + let column = |local: usize| f + local; + let mut scratch = Scratch { a: vec![0.0; m], b: vec![0.0; m], w: Vec::with_capacity(PANEL + 1) }; + let mut c0 = 0; + while c0 < w { + let c1 = own.panel(m, w, c0, &scale, &column, &mut scratch)?; + if c1 < w { + // The rest of the supernode minus the panel's L D Lᵀ. + let tiles = products(&Panel { m, data: &own.values }, c0..c1, c1, m - c1, w - c1, own.pivots.get(c0..c1)); + subtract(&mut own.values, m, &tiles, m - c1, w - c1, |i| c1 + i); + } + c0 = c1; + } + Ok(()) +} + +/// A symmetric indefinite factor P A Pᵀ = L D Lᵀ in supernodes. +#[derive(Debug, Clone)] +pub struct LdltFactor { + pub n: usize, + perm: Vec, + structure: Structure, + blocks: Vec, + entries: usize, + inertia: Inertia, + two_by_two: usize, +} + +impl LdltFactor { + pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result { + LdltFactor::from_symbolic(&Symbolic::new(a, max_entries)?) + } + + pub(crate) fn from_symbolic(symbolic: &Symbolic) -> std::result::Result { + let mut blocks = symbolic.allocate()?.into_iter().map(|values| Block { values, ..Block::default() }).collect::>(); + symbolic.factor(&mut blocks, supernode)?; + let mut inertia = Inertia::default(); + let mut two_by_two = 0; + for b in &blocks { + inertia.positive += b.inertia.positive; + inertia.negative += b.inertia.negative; + two_by_two += b.two_by_two; + } + Ok(LdltFactor { + n: symbolic.n, + perm: symbolic.perm.clone(), + structure: symbolic.structure.clone(), + blocks, + entries: symbolic.entries, + inertia, + two_by_two, + }) + } + + /// Entries in L (as the Cholesky factor of the same pattern). + pub fn entries(&self) -> usize { + self.entries + } + + pub fn inertia(&self) -> Inertia { + self.inertia + } + + /// The number of 2 × 2 pivots. + pub fn two_by_two(&self) -> usize { + self.two_by_two + } + + /// Solve A x = b in place. + pub fn solve(&self, b: &mut [f64]) -> Result<()> { + if b.len() != self.n { + return Err(Error::Internal("right-hand side length".into())); + } + let mut x: Vec = self.perm.iter().map(|&j| b.get(j).copied().unwrap_or(0.0)).collect(); + let s = &self.structure; + let get = |v: &[f64], i: usize| v.get(i).copied().unwrap_or(0.0); + let mut y = Vec::new(); + // L z = P b, then D⁻¹, by supernodes; the supernode's own unknowns are left in pivot + // order for the backward pass. + for (k, blk) in self.blocks.iter().enumerate() { + let rows = s.rows(k); + let (m, f, w) = (rows.len(), s.cols(k).start, s.cols(k).len()); + let t = Panel { m, data: &blk.values }; + y.clear(); + y.extend(blk.order.iter().map(|c| get(&x, f + c))); + let mut c = 0; + while c < w { + let width = if blk.pivots.get(c) == Some(&Pivot::First) { 2 } else { 1 }; + for q in c..c + width { + let yq = get(&y, q); + for (r, l) in (c + width..m).zip(t.column(q, c + width, m - c - width)) { + let target = if r < w { y.get_mut(r) } else { rows.get(r).and_then(|g| x.get_mut(*g)) }; + if let Some(v) = target { + *v -= l * yq; + } + } + } + c += width; + } + let mut c = 0; + while c < w { + if blk.pivots.get(c) == Some(&Pivot::First) { + let (d11, d21, d22) = (t.at(c, c), t.at(c + 1, c), t.at(c + 1, c + 1)); + let det = d11 * d22 - d21 * d21; + let (p, q) = (get(&y, c), get(&y, c + 1)); + if let Some(v) = y.get_mut(c) { + *v = (d22 * p - d21 * q) / det; + } + if let Some(v) = y.get_mut(c + 1) { + *v = (d11 * q - d21 * p) / det; + } + c += 2; + } else { + if let Some(v) = y.get_mut(c) { + *v /= t.at(c, c); + } + c += 1; + } + } + if let Some(dst) = x.get_mut(f..f + w) { + dst.copy_from_slice(&y); + } + } + // Lᵀ x = z, supernodes in reverse; each scatters its unknowns back from pivot order. + for (k, blk) in self.blocks.iter().enumerate().rev() { + let rows = s.rows(k); + let (m, f, w) = (rows.len(), s.cols(k).start, s.cols(k).len()); + let t = Panel { m, data: &blk.values }; + y.clear(); + y.extend(x.get(f..f + w).unwrap_or_default()); + let mut c = w; + while c > 0 { + let width = if c >= 2 && blk.pivots.get(c - 1) == Some(&Pivot::Second) { 2 } else { 1 }; + let first = c - width; + for q in first..c { + let mut sum = get(&y, q); + for (r, l) in (c..m).zip(t.column(q, c, m - c)) { + let v = if r < w { get(&y, r) } else { rows.get(r).map_or(0.0, |g| get(&x, *g)) }; + sum -= l * v; + } + if let Some(v) = y.get_mut(q) { + *v = sum; + } + } + c = first; + } + for (pos, c) in blk.order.iter().enumerate() { + if let Some(v) = x.get_mut(f + c) { + *v = get(&y, pos); + } + } + } + for (&j, v) in self.perm.iter().zip(&x) { + if let Some(t) = b.get_mut(j) { + *t = *v; + } + } + Ok(()) + } +} + +/// Which factorization a symmetric solve uses. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum Kind { + /// Cholesky; a matrix that is not positive definite is refused. + #[default] + Definite, + /// LDLᵀ with pivoting. + Indefinite, + /// Cholesky, and LDLᵀ when the matrix is not positive definite (one symbolic factorization + /// serves both). + Auto, +} + +/// A factored symmetric matrix: Cholesky or LDLᵀ. +#[derive(Debug, Clone)] +pub enum Symmetric { + Cholesky(Factor), + Indefinite(LdltFactor), +} + +impl Symmetric { + pub fn new(a: &UpperColumns, max_entries: usize, kind: Kind) -> std::result::Result { + match kind { + Kind::Definite => Factor::new(a, max_entries).map(Symmetric::Cholesky), + Kind::Indefinite => LdltFactor::new(a, max_entries).map(Symmetric::Indefinite), + Kind::Auto => { + let symbolic = Symbolic::new(a, max_entries)?; + match Factor::from_symbolic(&symbolic) { + Ok(f) => Ok(Symmetric::Cholesky(f)), + Err(FactorError::NotPositiveDefinite { .. }) => LdltFactor::from_symbolic(&symbolic).map(Symmetric::Indefinite), + Err(e) => Err(e), + } + } + } + } + + pub fn n(&self) -> usize { + match self { + Symmetric::Cholesky(f) => f.n, + Symmetric::Indefinite(f) => f.n, + } + } + + pub fn entries(&self) -> usize { + match self { + Symmetric::Cholesky(f) => f.entries(), + Symmetric::Indefinite(f) => f.entries(), + } + } + + /// The inertia: a Cholesky factor's matrix is positive definite. + pub fn inertia(&self) -> Inertia { + match self { + Symmetric::Cholesky(f) => Inertia { positive: f.n, negative: 0 }, + Symmetric::Indefinite(f) => f.inertia(), + } + } + + pub fn solve(&self, b: &mut [f64]) -> Result<()> { + match self { + Symmetric::Cholesky(f) => f.solve(b), + Symmetric::Indefinite(f) => f.solve(b), + } + } + + /// Solve A x = b in place, then refine with the residual against `a` (the matrix that was + /// factored, both triangles) while it shrinks, at most `steps` times (iterative refinement + /// in working precision: N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd + /// ed., SIAM 2002, §12.2). Cholesky solutions are returned as solved. + pub fn solve_refined(&self, a: &Csr, b: &mut [f64], steps: usize) -> Result<()> { + let rhs = b.to_vec(); + self.solve(b)?; + if matches!(self, Symmetric::Cholesky(_)) { + return Ok(()); + } + let norm = |v: &[f64]| v.iter().fold(0.0f64, |m, x| m.max(x.abs())); + let mut ax = vec![0.0; a.n]; + let residual = |x: &[f64], ax: &mut Vec| -> Vec { + a.mul(x, ax); + rhs.iter().zip(ax.iter()).map(|(p, q)| p - q).collect() + }; + let mut r = residual(b, &mut ax); + let mut size = norm(&r); + for _ in 0..steps { + if !(size > 0.0) { + break; + } + let mut dx = r.clone(); + self.solve(&mut dx)?; + let trial: Vec = b.iter().zip(&dx).map(|(x, d)| x + d).collect(); + let next = residual(&trial, &mut ax); + let next_size = norm(&next); + if !(next_size < size) { + break; + } + b.copy_from_slice(&trial); + (r, size) = (next, next_size); + } + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use fieldcraft_parallel as parallel; + + fn csr(rows: &[Vec<(usize, f64)>]) -> Csr { + let mut ptr = vec![0]; + let mut idx = Vec::new(); + let mut val = Vec::new(); + for row in rows { + let mut row = row.clone(); + row.sort_by_key(|e| e.0); + for (c, v) in row { + idx.push(c); + val.push(v); + } + ptr.push(idx.len()); + } + Csr { n: rows.len(), ptr, idx, val } + } + + fn residual(a: &Csr, f: &impl Fn(&mut [f64]), seed: f64) -> f64 { + let n = a.n; + let x: Vec = (0..n).map(|i| ((i as f64 + seed) * 0.37).sin()).collect(); + let mut b = vec![0.0; n]; + a.mul(&x, &mut b); + f(&mut b); + b.iter().zip(&x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max) + } + + /// A 3D grid Laplacian (3 coupled unknowns per point) shifted by σ: indefinite for σ inside + /// its spectrum, with a known count of eigenvalues below σ. + fn shifted_grid(n: [usize; 3], sigma: f64) -> (Csr, usize) { + let [nx, ny, nz] = n; + let point = |i: usize, j: usize, k: usize| i + nx * (j + ny * k); + let size = 3 * nx * ny * nz; + let mut rows: Vec> = vec![Vec::new(); size]; + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let p = 3 * point(i, j, k); + for a in 0..3 { + rows[p + a].push((p + a, 6.0 - sigma)); + } + let next = [ + (i + 1 < nx).then(|| point(i + 1, j, k)), + (j + 1 < ny).then(|| point(i, j + 1, k)), + (k + 1 < nz).then(|| point(i, j, k + 1)), + ]; + for q in next.into_iter().flatten() { + for a in 0..3 { + rows[p + a].push((3 * q + a, -1.0)); + rows[3 * q + a].push((p + a, -1.0)); + } + } + } + } + } + // Eigenvalues of the Dirichlet grid Laplacian: Σ 2 − 2 cos(π m/(n + 1)) over the axes. + let axis = |n: usize| (1..=n).map(move |m| 2.0 - 2.0 * (std::f64::consts::PI * m as f64 / (n as f64 + 1.0)).cos()); + let mut below = 0; + for a in axis(nx) { + for b in axis(ny) { + for c in axis(nz) { + if a + b + c < sigma { + below += 3; + } + } + } + } + (csr(&rows), below) + } + + #[test] + fn indefinite_systems_solve_with_their_inertia() { + for sigma in [0.0, 0.9, 3.3, 7.1] { + let (a, below) = shifted_grid([7, 6, 5], sigma); + let upper = UpperColumns::from_symmetric(&a); + let f = LdltFactor::new(&upper, usize::MAX).unwrap(); + assert_eq!(f.inertia(), Inertia { positive: a.n - below, negative: below }, "σ = {sigma}"); + let err = residual(&a, &|b: &mut [f64]| f.solve(b).unwrap(), sigma); + assert!(err < 1e-10, "σ = {sigma}: {err}"); + } + // A positive definite matrix: the same solution as Cholesky, no negative pivot. + let (a, _) = shifted_grid([5, 5, 5], -0.5); + let upper = UpperColumns::from_symmetric(&a); + let f = LdltFactor::new(&upper, usize::MAX).unwrap(); + assert_eq!(f.inertia().negative, 0); + assert_eq!(f.entries(), Factor::new(&upper, usize::MAX).unwrap().entries(), "the same fill"); + assert!(residual(&a, &|b: &mut [f64]| f.solve(b).unwrap(), 1.0) < 1e-12); + } + + #[test] + fn two_by_two_pivots_for_zero_diagonals() { + // A saddle point [[K, B], [Bᵀ, 0]] inside one dense supernode: zero diagonals need 2 × 2 + // pivots, and the inertia is (positive: K's size, negative: the constraints'). + let (nk, nb) = (8, 3); + let n = nk + nb; + let mut rows: Vec> = vec![Vec::new(); n]; + for i in 0..nk { + for j in 0..nk { + // Small next to the constraints' entries, so the search pairs them. + let v = if i == j { 0.5 } else { 0.05 / (1.0 + (i + j) as f64) }; + rows[i].push((j, v)); + } + } + for c in 0..nb { + for i in 0..nk { + let v = ((i * 7 + c * 3) % 5) as f64 - 2.0 + if i == c { 3.0 } else { 0.0 }; + rows[i].push((nk + c, v)); + rows[nk + c].push((i, v)); + } + // Stored zeros: the constraints share one supernode with K (a zero diagonal whose + // couplings all lie in later supernodes is the restriction the module documents). + for d in 0..nb { + rows[nk + c].push((nk + d, 0.0)); + } + } + // The zero-diagonal columns first, so the search must pivot past them. + let order: Vec = (nk..n).chain(0..nk).collect(); + let mut permuted: Vec> = vec![Vec::new(); n]; + let inverse: Vec = (0..n).map(|i| order.iter().position(|&o| o == i).unwrap()).collect(); + for (i, row) in rows.iter().enumerate() { + permuted[inverse[i]] = row.iter().map(|(c, v)| (inverse[*c], *v)).collect(); + } + let a = csr(&permuted); + let f = LdltFactor::new(&UpperColumns::from_symmetric(&a), usize::MAX).unwrap(); + assert_eq!(f.inertia(), Inertia { positive: nk, negative: nb }); + assert!(f.two_by_two() > 0); + assert!(residual(&a, &|b: &mut [f64]| f.solve(b).unwrap(), 0.5) < 1e-12); + } + + #[test] + fn singular_and_hostile_matrices_are_located_errors() { + // A Laplacian with the constant in its null space. + let n = 12; + let rows: Vec> = (0..n) + .map(|i: usize| { + let mut r = vec![]; + let mut d = 0.0; + for j in [i.wrapping_sub(1), i + 1] { + if j < n { + r.push((j, -1.0)); + d += 1.0; + } + } + r.push((i, d)); + r + }) + .collect(); + let e = LdltFactor::new(&UpperColumns::from_symmetric(&csr(&rows)), usize::MAX).unwrap_err(); + assert!(matches!(e, FactorError::Singular { column } if column < n), "{e:?}"); + // NaN. + let mut bad = rows.clone(); + bad[5][0].1 = f64::NAN; + bad[4].iter_mut().filter(|e| e.0 == 5).for_each(|e| e.1 = f64::NAN); + assert!(matches!(LdltFactor::new(&UpperColumns::from_symmetric(&csr(&bad)), usize::MAX), Err(FactorError::Singular { .. }))); + // The cap is checked before the numeric factorization; empty and 1 × 1 matrices. + let (a, _) = shifted_grid([3, 3, 3], 2.0); + assert!(matches!(LdltFactor::new(&UpperColumns::from_symmetric(&a), 10), Err(FactorError::TooLarge { .. }))); + let empty = Csr { n: 0, ptr: vec![0], idx: vec![], val: vec![] }; + LdltFactor::new(&UpperColumns::from_symmetric(&empty), 10).unwrap().solve(&mut []).unwrap(); + let f = LdltFactor::new(&UpperColumns::from_symmetric(&csr(&[vec![(0, -4.0)]])), 10).unwrap(); + let mut b = [2.0]; + f.solve(&mut b).unwrap(); + assert_eq!((b, f.inertia()), ([-0.5], Inertia { positive: 0, negative: 1 })); + } + + #[test] + fn same_factor_on_any_number_of_threads() { + // Large enough for parallel subtrees and tiles, with wide supernodes spanning panels. + let (a, below) = shifted_grid([14, 12, 10], 2.7); + let upper = UpperColumns::from_symmetric(&a); + let one = parallel::with_threads(1, || LdltFactor::new(&upper, usize::MAX)).unwrap(); + assert_eq!(one.inertia().negative, below); + let widest = one.blocks.iter().map(|b| b.order.len()).max().unwrap(); + assert!(widest > 2 * PANEL, "a supernode spans panels: {widest}"); + for threads in [2, 7] { + let other = parallel::with_threads(threads, || LdltFactor::new(&upper, usize::MAX)).unwrap(); + for (x, y) in one.blocks.iter().zip(&other.blocks) { + assert!(x.values.iter().zip(&y.values).all(|(p, q)| p.to_bits() == q.to_bits()) && x.order == y.order, "{threads} threads"); + } + } + let err = residual(&a, &|b: &mut [f64]| one.solve(b).unwrap(), 0.1); + assert!(err < 1e-9, "{err}"); + } + + #[test] + fn auto_selection_and_refinement() { + let (spd, _) = shifted_grid([4, 4, 4], -1.0); + let s = Symmetric::new(&UpperColumns::from_symmetric(&spd), usize::MAX, Kind::Auto).unwrap(); + assert!(matches!(s, Symmetric::Cholesky(_))); + let (ind, below) = shifted_grid([4, 4, 4], 4.4); + let s = Symmetric::new(&UpperColumns::from_symmetric(&ind), usize::MAX, Kind::Auto).unwrap(); + assert!(matches!(s, Symmetric::Indefinite(_))); + assert_eq!(s.inertia().negative, below); + let x: Vec = (0..ind.n).map(|i| (i as f64).cos()).collect(); + let mut b = vec![0.0; ind.n]; + ind.mul(&x, &mut b); + s.solve_refined(&ind, &mut b, 3).unwrap(); + assert!(b.iter().zip(&x).all(|(p, q)| (p - q).abs() < 1e-12)); + assert!(matches!( + Symmetric::new(&UpperColumns::from_symmetric(&ind), usize::MAX, Kind::Definite), + Err(FactorError::NotPositiveDefinite { .. }) + )); + } + + /// The stiffness pattern of a hex8 mesh (27 neighbours, 3 unknowns per node) on an n³ grid of + /// nodes, ordered by nested dissection, diagonally dominant and shifted by σ. + fn hex_pattern(n: usize, sigma: f64) -> Csr { + let id = |i: usize, j: usize, k: usize| i + n * (j + n * k); + let nodes = n * n * n; + let xyz: Vec<[f64; 3]> = (0..nodes).map(|v| [(v % n) as f64, ((v / n) % n) as f64, (v / (n * n)) as f64]).collect(); + let mut adjacency = vec![Vec::new(); nodes]; + for k in 0..n { + for j in 0..n { + for i in 0..n { + for (a, b, c) in (0..27).map(|t| (t % 3, (t / 3) % 3, t / 9)) { + let (x, y, z) = (i + a, j + b, k + c); + if (1..=n).contains(&x) && (1..=n).contains(&y) && (1..=n).contains(&z) { + adjacency[id(i, j, k)].push(id(x - 1, y - 1, z - 1)); + } + } + } + } + } + let order = crate::ordering::nested_dissection(&xyz, &adjacency); + let mut place = vec![0; nodes]; + for (p, v) in order.iter().enumerate() { + place[*v] = p; + } + let mut rows: Vec> = vec![Vec::new(); 3 * nodes]; + for (v, adj) in adjacency.iter().enumerate() { + for w in adj { + for a in 0..3 { + for b in 0..3 { + let value = if v == *w { + if a == b { 30.0 - sigma } else { 0.5 } + } else if a == b { + -1.0 + } else { + 0.01 * ((v + w) % 3) as f64 + }; + rows[3 * place[v] + a].push((3 * place[*w] + b, value)); + } + } + } + } + csr(&rows) + } + + /// Factor time and fill of Cholesky and LDLᵀ on a large stiffness pattern (`--nocapture`). + #[test] + #[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] + fn speed_against_cholesky() { + for n in [24, 36] { + let a = hex_pattern(n, 0.0); + let upper = UpperColumns::from_symmetric(&a); + let t = std::time::Instant::now(); + let chol = Factor::new(&upper, usize::MAX).unwrap(); + let tc = t.elapsed().as_secs_f64(); + let t = std::time::Instant::now(); + let ldlt = LdltFactor::new(&upper, usize::MAX).unwrap(); + let tl = t.elapsed().as_secs_f64(); + assert_eq!((ldlt.entries(), ldlt.inertia().negative), (chol.entries(), 0)); + let shifted = hex_pattern(n, 22.0); + let t = std::time::Instant::now(); + let indefinite = LdltFactor::new(&UpperColumns::from_symmetric(&shifted), usize::MAX).unwrap(); + let ti = t.elapsed().as_secs_f64(); + assert!(residual(&shifted, &|b: &mut [f64]| indefinite.solve(b).unwrap(), 0.3) < 1e-8); + eprintln!( + "{} equations, {} entries: Cholesky {tc:.2} s, LDLᵀ {tl:.2} s (×{:.2}); indefinite (σ = 22: {} negative, {} 2 × 2 pivots) {ti:.2} s", + a.n, + chol.entries(), + tl / tc, + indefinite.inertia().negative, + indefinite.two_by_two() + ); + } + } +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index b7aca2a..36d35a9 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -13,7 +13,9 @@ //! 4. **Solution** ([`cholesky`]): sparse Cholesky; reactions K u − f at the prescribed degrees //! of freedom; stresses at the integration points, extrapolated and averaged at the nodes, and //! von Mises. Frequency analysis ([`eigen`]) runs shift-invert block Lanczos on the same -//! kind of factorization. +//! kind of factorization. Indefinite matrices (tangents past a limit point, strained mixed +//! elements, shifted operators) are factored by the symmetric indefinite LDLᵀ ([`ldlt`]) on +//! the same symbolic factorization, which also counts their negative eigenvalues. //! //! Element matrices and stresses, the factorization's subtrees and dense tiles, the Lanczos //! block solves and matrix products run on threads ([`fieldcraft_parallel`]), but every number @@ -31,6 +33,7 @@ pub mod constraints; pub mod discrete; pub mod eigen; pub mod heat; +pub mod ldlt; pub mod nonlinear; pub mod ordering; pub mod sparse; diff --git a/docs/fem.md b/docs/fem.md index 8b31064..1a5ed6f 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -150,6 +150,29 @@ when confined. A continuous nodal pressure (`DofKind::P`, Taylor–Hood) needs a symmetric indefinite factorization and is planned. +## Symmetric indefinite factorization (LDLᵀ) + +`fieldcraft_solver::ldlt` factors P A Pᵀ = L D Lᵀ (L unit lower triangular, D with 1 × 1 and +2 × 2 pivots) on the Cholesky factor's symbolic factorization: the same nested dissection +ordering, elimination tree, supernodes, fill, factor-size cap and parallel subtree tasks, so it +is deterministic on any number of threads. Pivots are chosen by the bounded Bunch–Kaufman (rook) +search of Ashcraft, Grimes & Lewis (1998) among the columns of the supernode being factored +(Schenk & Gärtner's supernode Bunch–Kaufman, 2006): an interchange inside a supernode changes no +row structure, so nothing is delayed and the static structure holds. Its limits, stated in the +module: L's entries below the supernode are not bounded by the search (`solve_refined` adds +iterative refinement), and a column whose supernode offers no pivot above round-off is a located +`Singular` error (a matrix singular at that column, or a zero diagonal coupled only to later +supernodes, which no FieldCraft system has: hybrid pressures are condensed per element). + +The factor reports its **inertia** (positive and negative eigenvalue counts, Sylvester's law): +a 2 × 2 pivot of the search has a negative determinant (one of each). `Symmetric::new(a, cap, +kind)` picks Cholesky (`Definite`, the default everywhere a matrix is positive definite, so +linear results are bit for bit unchanged), LDLᵀ (`Indefinite`), or Cholesky with LDLᵀ when the +matrix turns out not to be positive definite (`Auto`, one symbolic factorization for both). +On a hex8 stiffness pattern of 140 000 equations (1.5 · 10⁸ entries) the LDLᵀ factor costs 1.1 +times the Cholesky factor, and 1.2 times on a shift with 10 600 negative eigenvalues and 7 700 +2 × 2 pivots (`ldlt::tests::speed_against_cholesky`). + ## Nonlinear static analysis (finite strain) A step `{"type": "nonlinear_static", …}` (`Procedure::NonlinearStatic`, solved by diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 5f4a765..10e6f1e 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1228,3 +1228,21 @@ update are written from (cited again in the module docs): - Thermal strain by the multiplicative split F = F_m F_θ: V. A. Lubarda, *Constitutive theories based on the multiplicative decomposition of deformation gradient: thermoelasticity, elastoplasticity, and biomechanics*, Appl. Mech. Rev. 57(2):95–108, 2004. + +## Part P: Symmetric indefinite factorization, linear buckling and arc-length, as implemented + +- Symmetric indefinite pivoting: J. R. Bunch, L. Kaufman, *Some stable methods for calculating + inertia and solving symmetric linear systems*, Math. Comp. 31(137):163–179, 1977; the bounded + (rook) search, C. Ashcraft, R. G. Grimes, J. G. Lewis, *Accurate symmetric indefinite linear + equation solvers*, SIAM J. Matrix Anal. Appl. 20(2):513–561, 1998; the panel's partial + factorization of LAPACK's `dsytrf`, E. Anderson et al., *LAPACK Users' Guide*, 3rd ed., SIAM + 1999. +- Pivoting within supernodes on a static structure: O. Schenk, K. Gärtner, *On fast + factorization pivoting methods for sparse symmetric indefinite systems*, Electron. Trans. + Numer. Anal. 23:158–179, 2006. Delayed pivots (the multifrontal alternative, not used: it + changes the structure as it goes): I. S. Duff, J. K. Reid, *The multifrontal solution of + indefinite sparse symmetric linear equations*, ACM TOMS 9(3):302–325, 1983; J. W. H. Liu, *The + multifrontal method for sparse matrix solution: theory and practice*, SIAM Review + 34(1):82–109, 1992. +- Iterative refinement: N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., + SIAM 2002, §12.2.