solver: symmetric indefinite LDLᵀ factorization with inertia

A supernodal LDLᵀ on the Cholesky factor's symbolic factorization (ordering, supernodes, fill,
cap and parallel subtree tasks shared through a Symbolic and a generic driver; the Cholesky
numerics are unchanged bit for bit). Pivots by the bounded Bunch–Kaufman (rook) search of
Ashcraft, Grimes & Lewis among each supernode's columns (Schenk & Gärtner's supernode
Bunch–Kaufman), 1 × 1 and 2 × 2, right-looking panels whose trailing update reuses the tiled
Cholesky kernels with L D formed on the fly. Reports the inertia; a column without a pivot above
round-off is a located Singular error; NaN is refused. Symmetric::new picks Cholesky, LDLᵀ, or
Cholesky with an LDLᵀ fallback; solve_refined adds iterative refinement. 1.1–1.2× the Cholesky
time on a 140 000-equation hex8 pattern.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

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)"
+            ))
+        }
     }
 }
 

crates/solver/src/cholesky.rs

@@ -211,6 +211,10 @@ fn postorder(parent: &[usize]) -> Vec<usize> {
 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<usize> {
+        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<usize>,
     /// 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<usize> {
+    pub(crate) fn cols(&self, s: usize) -> std::ops::Range<usize> {
         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> {
+pub(crate) struct Panel<'a> {
-    m: usize,
+    pub(crate) m: usize,
-    data: &'a [f64],
+    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
+/// The pivots of a symmetric indefinite factor's block diagonal D, per column of a supernode
-/// 0..mu and columns c in 0..wu (meaningful for i ≥ c), each sum in column order from zero.
+/// ([`crate::ldlt`]): a 1 × 1 pivot, or the first or second column of a 2 × 2 one. D's entries
-/// Returned by row tiles: tile t holds rows t·TILE.. as `tile[c * TILE + i]`.
+/// sit on the trapezoid's diagonal, and a 2 × 2 pivot's off-diagonal entry below its first
-fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec<Vec<f64>> {
+/// column's diagonal (where L's unit block has a zero).
-    // Rows p..p+wu of the kw columns, transposed so that column c's coefficients are contiguous.
+#[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<usize>, p: usize, mu: usize, wu: usize, pivots: Option<&[Pivot]>) -> Vec<Vec<f64>> {
+    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<f64> {
         let i0 = t * TILE;
@@ -305,7 +372,7 @@ fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec<Vec<f64
         let mut x = vec![0.0; kw * TILE];
         for k in 0..kw {
             if let Some(dst) = x.get_mut(k * TILE..k * TILE + h) {
-                dst.copy_from_slice(l.column(k, p + i0, h));
+                dst.copy_from_slice(l.column(k0 + k, p + i0, h));
             }
         }
         let cols = wu.min(i0 + h);
@@ -328,7 +395,7 @@ fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec<Vec<f64
 }
 
 /// own(rel(i), rel(c)) −= S(i, c) for i ≥ c, with S in the tiles of [`products`].
-fn subtract(own: &mut [f64], m: usize, tiles: &[Vec<f64>], mu: usize, wu: usize, rel: impl Fn(usize) -> usize) {
+pub(crate) fn subtract(own: &mut [f64], m: usize, tiles: &[Vec<f64>], 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<Vec<f64>> {
 }
 
 /// 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<usize>,
     idx: Vec<usize>,
     val: Vec<f64>,
     diag: Vec<f64>,
+    /// The largest magnitude in each row (= column) of the whole symmetric matrix.
+    pub(crate) scale: Vec<f64>,
 }
 
 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> {
+pub(crate) struct Context<'a> {
-    a: &'a Lower,
+    pub(crate) a: &'a Lower,
-    s: &'a Structure,
+    pub(crate) s: &'a Structure,
     updates: &'a [Vec<Update>],
 }
 
+/// 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<f64> {
+    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
+    /// Supernode `j` assembled into `own`: its columns of A minus the update L_K Y_Kᵀ of every
-    /// the first column (postorder) whose pivot is not positive.
+    /// descendant K that reaches it (Y = L, or L D), in supernode order; the descendants are
-    fn supernode(&self, j: usize, done: &[Vec<f64>], base: usize, own: &mut [f64]) -> std::result::Result<(), FactorError> {
+    /// `done[k - base]`.
-        let internal = |what: &str| FactorError::Internal(format!("supernodal factorization: {what}"));
+    pub(crate) fn assemble<T: Supernode>(&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<f64>], 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.
+/// The symbolic factorization shared by the Cholesky and the symmetric indefinite (LDLᵀ,
-#[derive(Debug, Clone)]
+/// [`crate::ldlt`]) factors: the postorder, the supernodes and their row structures, the updates
-pub struct Factor {
+/// each supernode receives, the matrix in that order, and the subtree tasks of the numeric
-    pub n: usize,
+/// phase.
+pub(crate) struct Symbolic {
+    pub(crate) n: usize,
     /// `perm[k]`: the original column that is column k of the factor.
-    perm: Vec<usize>,
+    pub(crate) perm: Vec<usize>,
-    structure: Structure,
+    pub(crate) structure: Structure,
-    /// The values of each supernode (a lower trapezoid by columns).
+    a: Lower,
-    values: Vec<Vec<f64>>,
+    updates: Vec<Vec<Update>>,
-    entries: usize,
+    /// Subtrees factored as parallel tasks, and the supernodes above them, in order.
+    tasks: Vec<std::ops::Range<usize>>,
+    top: Vec<usize>,
+    /// Entries of L.
+    pub(crate) entries: usize,
 }
 
-impl Factor {
+impl Symbolic {
-    pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result<Factor, FactorError> {
+    /// 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<Symbolic, FactorError> {
         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<f64>> = 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<usize> = (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<Vec<Vec<f64>>, FactorError> {
+        let ns = self.supernodes();
+        let mut values: Vec<Vec<f64>> = 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<T: Supernode>(
+        &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<usize>, &mut [Vec<f64>], Option<usize>)> = Vec::with_capacity(tasks.len());
+        let mut pieces: Vec<(std::ops::Range<usize>, &mut [T], Option<FactorError>)> = Vec::with_capacity(self.tasks.len());
-        let mut rest: &mut [Vec<f64>] = &mut values;
+        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::<FactorError>);
-        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 }) => {
+                    Err(e) if e.column().is_some() => {
-                        *singular = Some(column);
+                        *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);
+        let named = |e: FactorError| match e.column() {
-        for &j in &top {
+            Some(c) => e.at(self.original(c)),
-            if earliest.is_some_and(|c| structure.cols(j).start > 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) {
+            supernode(&context, j, done, 0, own).map_err(named)?;
-                Ok(()) => {}
-                Err(FactorError::NotPositiveDefinite { column }) => return Err(FactorError::NotPositiveDefinite { column: original(column) }),
-                Err(e) => return Err(e),
-            }
         }
-        if let Some(column) = earliest {
+        match earliest {
-            return Err(FactorError::NotPositiveDefinite { column: original(column) });
+            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<usize>,
+    structure: Structure,
+    /// The values of each supernode (a lower trapezoid by columns).
+    values: Vec<Vec<f64>>,
+    entries: usize,
+}
+
+impl Factor {
+    pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result<Factor, FactorError> {
+        Factor::from_symbolic(&Symbolic::new(a, max_entries)?)
+    }
+
+    pub(crate) fn from_symbolic(symbolic: &Symbolic) -> std::result::Result<Factor, FactorError> {
+        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.

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"

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<f64>,
+    pivots: Vec<Pivot>,
+    /// `order[c]`: the supernode's own column (0..w, in the symbolic order) at pivot position c.
+    order: Vec<usize>,
+    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<usize>, skip: Option<usize>) -> (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<f64>,
+    b: Vec<f64>,
+    /// A row of L D over the panel's factored columns.
+    w: Vec<f64>,
+}
+
+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<f64>) {
+        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<usize, FactorError> {
+        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<usize>,
+    structure: Structure,
+    blocks: Vec<Block>,
+    entries: usize,
+    inertia: Inertia,
+    two_by_two: usize,
+}
+
+impl LdltFactor {
+    pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result<LdltFactor, FactorError> {
+        LdltFactor::from_symbolic(&Symbolic::new(a, max_entries)?)
+    }
+
+    pub(crate) fn from_symbolic(symbolic: &Symbolic) -> std::result::Result<LdltFactor, FactorError> {
+        let mut blocks = symbolic.allocate()?.into_iter().map(|values| Block { values, ..Block::default() }).collect::<Vec<_>>();
+        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<f64> = 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<Symmetric, FactorError> {
+        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<f64>| -> Vec<f64> {
+            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<f64> = 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<f64> = (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<(usize, f64)>> = 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<(usize, f64)>> = 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<usize> = (nk..n).chain(0..nk).collect();
+        let mut permuted: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
+        let inverse: Vec<usize> = (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<Vec<(usize, f64)>> = (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<f64> = (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<(usize, f64)>> = 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()
+            );
+        }
+    }
+}

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;

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

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.