Merge r8g-fea-sparse-solver (R8/R7 phase 2 round 1; default-off, verified)

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

crates/specialized/rtx-fea/src/analysis/mod.rs

@@ -11,6 +11,7 @@ pub mod flag3d;
 pub mod modal_analysis;
 pub mod nonlinear_analysis;
 pub mod nonlinear_dynamic;
+pub mod sparse_tangent;
 pub mod static_analysis;
 
 use crate::assembly::AssemblyOptions;
@@ -26,6 +27,7 @@ pub use dynamic_analysis::*;
 pub use modal_analysis::*;
 pub use nonlinear_analysis::*;
 pub use nonlinear_dynamic::*;
+pub use sparse_tangent::{TangentSolver, TangentStats};
 pub use static_analysis::*;
 
 /// Base trait for all finite element analyses.

crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs

@@ -45,6 +45,7 @@
 //! change between steps (and between subiterations of one step) through
 //! the stepper.
 
+use super::sparse_tangent::{SparseTangent, TangentSolver, TangentStats};
 use super::{AnalysisConfig, ConvergenceCriteria};
 use crate::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
 use crate::assembly::SparseMatrix;
@@ -56,6 +57,7 @@ use crate::materials::{reduced_constitutive, MaterialDatabase};
 use crate::mesh::{Mesh, NodeId};
 use crate::solvers::{BandedLu, LinearSolver, SolverOptions};
 use nalgebra::{DMatrix, DVector, Vector3};
+use rayon::prelude::*;
 
 /// Time histories and final state of a nonlinear transient run.
 #[derive(Debug, Clone)]
@@ -118,6 +120,9 @@ pub struct NonlinearDynamicAnalysis {
     body_force: Option<Box<dyn Fn(Vector3<f64>) -> Vector3<f64> + Send + Sync>>,
     nodal_forces: Vec<(NodeId, Vector3<f64>)>,
     tracked_nodes: Vec<NodeId>,
+    /// Explicit tangent-solver choice; `None` = `RTX_FEA_TANGENT`, else
+    /// the banded LU.
+    tangent_solver: Option<TangentSolver>,
 }
 
 impl NonlinearDynamicAnalysis {
@@ -145,9 +150,26 @@ impl NonlinearDynamicAnalysis {
             body_force: None,
             nodal_forces: Vec::new(),
             tracked_nodes: Vec::new(),
+            tangent_solver: None,
         }
     }
 
+    /// Choose the linear solver of the Newton tangent (default: the
+    /// banded LU, or the `RTX_FEA_TANGENT` knob when set). See
+    /// [`TangentSolver`].
+    #[must_use]
+    pub fn with_tangent_solver(mut self, solver: TangentSolver) -> Self {
+        self.tangent_solver = Some(solver);
+        self
+    }
+
+    /// The tangent solver a new stepper will use.
+    pub fn tangent_solver(&self) -> TangentSolver {
+        self.tangent_solver
+            .or_else(TangentSolver::from_env)
+            .unwrap_or_default()
+    }
+
     /// Switch to the total-Lagrangian St. Venant–Kirchhoff formulation
     /// (plane strain in 2-D), as on the nonlinear static analysis.
     #[must_use]
@@ -284,6 +306,8 @@ pub struct NonlinearDynamicStepper<'a> {
     /// Steps carried by step subdivision after both the plain loop and
     /// the line search failed.
     rescued_subdivision: usize,
+    /// The sparse tangent path (`None` = the banded LU path).
+    sparse: Option<SparseTangent>,
 }
 
 impl<'a> NonlinearDynamicStepper<'a> {
@@ -420,6 +444,15 @@ impl<'a> NonlinearDynamicStepper<'a> {
         }
 
         let added_mass = DVector::zeros(total_dofs);
+        let sparse = match analysis.tangent_solver() {
+            TangentSolver::BandedLu => None,
+            TangentSolver::SparseLdlt { reuse } => Some(SparseTangent::new(
+                caches.iter().map(|c| c.dofs.as_slice()),
+                &free_index,
+                num_free,
+                reuse,
+            )?),
+        };
         let mut stepper = Self {
             added_mass,
             analysis,
@@ -435,6 +468,7 @@ impl<'a> NonlinearDynamicStepper<'a> {
             solver_options: SolverOptions::default(),
             rescued_line_search: 0,
             rescued_subdivision: 0,
+            sparse,
         };
         stepper.set_nodal_forces(&analysis.nodal_forces);
         Ok(stepper)
@@ -475,9 +509,16 @@ impl<'a> NonlinearDynamicStepper<'a> {
         // inv_beta_dt2 is only read when assembling the tangent.
         let (f_int0, _) = self.assemble(&u, false, 0.0)?;
         let residual0 = &self.external - &f_int0;
-        let (a0_free, _) = self
+        let (a0_free, _) = match &mut self.sparse {
+            None => self
                 .solver
-            .solve(&self.mass_free, &residual0, &self.solver_options)?;
+                .solve(&self.mass_free, &residual0, &self.solver_options)?,
+            Some(sparse) => {
+                sparse
+                    .mass_solver
+                    .solve(&self.mass_free, &residual0, &self.solver_options)?
+            }
+        };
         let mut a = DVector::zeros(self.total_dofs);
         for (k, &dof) in self.free_dofs.iter().enumerate() {
             a[dof] = a0_free[k];
@@ -602,12 +643,20 @@ impl<'a> NonlinearDynamicStepper<'a> {
         let mut u_iter = u_pred.clone();
         let mut step_converged = false;
         let mut iterations = 0usize;
+        // Sparse path only: the previous iteration's residual norm.
+        let mut previous_norm = f64::INFINITY;
         for _ in 0..criteria.max_iterations {
             let mut a_new = DVector::zeros(self.total_dofs);
             for &dof in &self.free_dofs {
                 a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]);
             }
+            let (f_int, tangent) = if let Some(sparse) = &self.sparse {
+                let with_tangent = sparse.wants_tangent(line_search, inv_beta_dt2);
+                (self.sparse_forces(&u_iter, with_tangent)?, None)
+            } else {
                 let (f_int, tangent) = self.assemble(&u_iter, true, inv_beta_dt2)?;
+                (f_int, Some(tangent))
+            };
             let residual = &self.external - &f_int - self.mass_times(&a_new);
             let residual_norm = residual.norm();
             if line_search && !residual_norm.is_finite() {
@@ -618,9 +667,35 @@ impl<'a> NonlinearDynamicStepper<'a> {
                 break;
             }
             iterations += 1;
-            let (delta, _) = self
+            let delta = match tangent {
-                .solver
+                Some(tangent) => {
+                    let (delta, _) =
+                        self.solver
                             .solve(&tangent, &residual, &self.solver_options)?;
+                    delta
+                }
+                None => {
+                    let rate = residual_norm / previous_norm;
+                    previous_norm = residual_norm;
+                    let sparse = self.sparse.as_ref().expect("sparse path");
+                    if sparse.needs_refresh(line_search, inv_beta_dt2, rate)
+                        && !sparse.tangents_current()
+                    {
+                        self.sparse_tangents(&u_iter)?;
+                    }
+                    let masses: Vec<&DMatrix<f64>> = self.caches.iter().map(|c| &c.mass).collect();
+                    let sparse = self.sparse.as_mut().expect("sparse path");
+                    sparse.solve(
+                        &masses,
+                        &self.added_mass,
+                        &self.free_index,
+                        &residual,
+                        inv_beta_dt2,
+                        line_search,
+                        rate,
+                    )?
+                }
+            };
             let alpha = if line_search {
                 match self.backtrack(&u_iter, &delta, &u_pred, inv_beta_dt2, residual_norm) {
                     Some(alpha) => alpha,
@@ -732,6 +807,63 @@ impl<'a> NonlinearDynamicStepper<'a> {
         (self.rescued_line_search, self.rescued_subdivision)
     }
 
+    /// Counters and timings of the sparse tangent path (`None` on the
+    /// banded path).
+    pub fn tangent_stats(&self) -> Option<TangentStats> {
+        self.sparse.as_ref().map(SparseTangent::stats)
+    }
+
+    /// Sparse path: every element's force and tangent in parallel, the
+    /// forces summed in element order (the banded path's order, so the
+    /// internal force is the same float for float), the tangents kept for
+    /// the factorisation.
+    /// With `with_tangent` false (a modified-Newton iteration that will
+    /// reuse the factor) only the forces are evaluated — the TL force-only
+    /// kernel gives the same forces bit for bit.
+    fn sparse_forces(
+        &mut self,
+        solution: &DVector<f64>,
+        with_tangent: bool,
+    ) -> FeaResult<DVector<f64>> {
+        let started = std::time::Instant::now();
+        let analysis = self.analysis;
+        let results: Vec<(DVector<f64>, Option<DMatrix<f64>>)> = self
+            .caches
+            .par_iter()
+            .map(|cache| {
+                if with_tangent {
+                    element_force_and_tangent(analysis, cache, solution).map(|(f, k)| (f, Some(k)))
+                } else {
+                    element_force(analysis, cache, solution).map(|f| (f, None))
+                }
+            })
+            .collect::<FeaResult<_>>()?;
+        let sparse = self.sparse.as_mut().expect("sparse path");
+        Ok(sparse.accumulate(
+            results,
+            self.caches.iter().map(|c| c.dofs.as_slice()),
+            &self.free_index,
+            started,
+        ))
+    }
+
+    /// Sparse path: the element tangents at `solution` alone (a refresh
+    /// after a force-only assembly).
+    fn sparse_tangents(&mut self, solution: &DVector<f64>) -> FeaResult<()> {
+        let started = std::time::Instant::now();
+        let analysis = self.analysis;
+        let tangents: Vec<DMatrix<f64>> = self
+            .caches
+            .par_iter()
+            .map(|cache| element_force_and_tangent(analysis, cache, solution).map(|(_, k)| k))
+            .collect::<FeaResult<_>>()?;
+        self.sparse
+            .as_mut()
+            .expect("sparse path")
+            .set_tangents(tangents, started);
+        Ok(())
+    }
+
     /// Internal force and (optionally) tangent at a full displacement
     /// vector, reduced to the free DOFs. The tangent includes the Newmark
     /// mass term `M / (β Δt²)` at the caller's `inv_beta_dt2` (the
@@ -742,41 +874,11 @@ impl<'a> NonlinearDynamicStepper<'a> {
         with_tangent: bool,
         inv_beta_dt2: f64,
     ) -> FeaResult<(DVector<f64>, SparseMatrix)> {
-        let dim = self.analysis.mesh.spatial_dimension;
         let num_free = self.free_dofs.len();
         let mut internal = DVector::zeros(num_free);
         let mut tangent = SparseMatrix::new(num_free, num_free);
         for cache in &self.caches {
-            let material = self
+            let (f_int, k_t) = element_force_and_tangent(self.analysis, cache, solution)?;
-                .analysis
-                .materials
-                .get_material(cache.material_id)
-                .unwrap();
-            let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
-            let mut element_displacement = DVector::zeros(cache.dofs.len());
-            for (local, &dof) in cache.dofs.iter().enumerate() {
-                element_displacement[local] = solution[dof];
-            }
-            let (f_int, k_t) = if self.analysis.total_lagrangian {
-                let (lambda, mu) = material.properties().lame_parameters();
-                let constitutive = saint_venant_kirchhoff(lambda, mu, dim);
-                total_lagrangian::internal_force_and_tangent(
-                    &fe,
-                    &cache.coords,
-                    &element_displacement,
-                    constitutive.as_ref(),
-                    None,
-                )?
-            } else {
-                let constitutive = reduced_constitutive(material, dim)?;
-                ElementMatrixComputer::compute_internal_force_and_tangent(
-                    &fe,
-                    &cache.coords,
-                    &element_displacement,
-                    constitutive.as_ref(),
-                    None,
-                )?
-            };
             for (local_row, &dof_row) in cache.dofs.iter().enumerate() {
                 let Some(free_row) = self.free_index[dof_row] else {
                     continue;
@@ -834,3 +936,68 @@ impl<'a> NonlinearDynamicStepper<'a> {
         out
     }
 }
+
+/// One element's internal force and tangent at the full displacement
+/// vector `solution` (total-Lagrangian SVK or small strain, per the
+/// analysis) — shared by the banded and the sparse paths.
+fn element_force_and_tangent(
+    analysis: &NonlinearDynamicAnalysis,
+    cache: &ElementCache,
+    solution: &DVector<f64>,
+) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
+    let dim = analysis.mesh.spatial_dimension;
+    let material = analysis.materials.get_material(cache.material_id).unwrap();
+    let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
+    let mut element_displacement = DVector::zeros(cache.dofs.len());
+    for (local, &dof) in cache.dofs.iter().enumerate() {
+        element_displacement[local] = solution[dof];
+    }
+    if analysis.total_lagrangian {
+        let (lambda, mu) = material.properties().lame_parameters();
+        let constitutive = saint_venant_kirchhoff(lambda, mu, dim);
+        total_lagrangian::internal_force_and_tangent(
+            &fe,
+            &cache.coords,
+            &element_displacement,
+            constitutive.as_ref(),
+            None,
+        )
+    } else {
+        let constitutive = reduced_constitutive(material, dim)?;
+        ElementMatrixComputer::compute_internal_force_and_tangent(
+            &fe,
+            &cache.coords,
+            &element_displacement,
+            constitutive.as_ref(),
+            None,
+        )
+    }
+}
+
+/// One element's internal force alone (the TL force-only kernel; the
+/// small-strain path evaluates the pair and drops the tangent).
+fn element_force(
+    analysis: &NonlinearDynamicAnalysis,
+    cache: &ElementCache,
+    solution: &DVector<f64>,
+) -> FeaResult<DVector<f64>> {
+    if !analysis.total_lagrangian {
+        return element_force_and_tangent(analysis, cache, solution).map(|(f, _)| f);
+    }
+    let dim = analysis.mesh.spatial_dimension;
+    let material = analysis.materials.get_material(cache.material_id).unwrap();
+    let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
+    let mut element_displacement = DVector::zeros(cache.dofs.len());
+    for (local, &dof) in cache.dofs.iter().enumerate() {
+        element_displacement[local] = solution[dof];
+    }
+    let (lambda, mu) = material.properties().lame_parameters();
+    let constitutive = saint_venant_kirchhoff(lambda, mu, dim);
+    total_lagrangian::internal_force(
+        &fe,
+        &cache.coords,
+        &element_displacement,
+        constitutive.as_ref(),
+        None,
+    )
+}

crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs

@@ -0,0 +1,379 @@
+//! The sparse Newton-tangent path of [`super::NonlinearDynamicStepper`]
+//! (R8-g): a fixed-pattern CSR assembled from per-element blocks through
+//! precomputed positions, factorised by [`SparseLdlt`] with the symbolic
+//! analysis done once (the pattern is fixed by the mesh), optionally
+//! reusing the numeric factor for several Newton iterations (modified
+//! Newton).
+//!
+//! Default off: [`TangentSolver::BandedLu`] is the stepper's original
+//! path, float for float. Select with
+//! [`super::NonlinearDynamicAnalysis::with_tangent_solver`] or the
+//! environment knob `RTX_FEA_TANGENT=sparse` (`sparse:K` = reuse each
+//! factor for up to K Newton solves; `banded` = the default).
+//!
+//! What stays identical to the banded path: the element kernels, the
+//! internal force (summed in the same element/local order), the
+//! residual, the convergence tests, the rescue ladder. What differs: the
+//! linear solve (LDLᵀ in a nested-dissection order vs banded LU with
+//! partial pivoting — rounding-level), and the tangent keeps entries the
+//! old `SparseMatrix` assembly drops below 1e-15.
+
+use crate::assembly::SparseMatrix;
+use crate::error::FeaResult;
+use crate::solvers::{BandedLu, LinearSolver, SolverOptions, SparseLdlt};
+use nalgebra::{DMatrix, DVector};
+use rayon::prelude::*;
+use std::time::Instant;
+
+/// CSR entries per parallel scatter task.
+const SCATTER_CHUNK: usize = 16_384;
+
+/// Which linear solver carries the Newton tangent.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum TangentSolver {
+    /// The original path: `SparseMatrix` assembly + banded LU every
+    /// Newton iteration.
+    #[default]
+    BandedLu,
+    /// Fixed-pattern parallel assembly + supernodal sparse LDLᵀ with
+    /// symbolic reuse. `reuse` = how many Newton solves one numeric
+    /// factor may serve (1 = full Newton; more = modified Newton, the
+    /// factor refreshed early when the dt changes, when the residual
+    /// contracts by less than half per iteration, and always in the
+    /// line-search rescue).
+    SparseLdlt {
+        /// Newton solves per numeric factorisation (≥ 1).
+        reuse: usize,
+    },
+}
+
+impl TangentSolver {
+    /// Full-Newton sparse path.
+    pub const SPARSE: Self = Self::SparseLdlt { reuse: 1 };
+
+    /// Parse the `RTX_FEA_TANGENT` knob: `banded`, `sparse`, `sparse:K`.
+    /// Unset or unparsable → `None` (unparsable values are logged).
+    pub fn from_env() -> Option<Self> {
+        let value = std::env::var("RTX_FEA_TANGENT").ok()?;
+        let value = value.trim();
+        let parsed = match value {
+            "banded" => Some(Self::BandedLu),
+            "sparse" => Some(Self::SPARSE),
+            other => other
+                .strip_prefix("sparse:")
+                .and_then(|k| k.parse::<usize>().ok())
+                .filter(|&k| k >= 1)
+                .map(|reuse| Self::SparseLdlt { reuse }),
+        };
+        if parsed.is_none() {
+            tracing::warn!("RTX_FEA_TANGENT={value:?} not understood; using the banded LU");
+        }
+        parsed
+    }
+}
+
+/// Counters and accumulated wall time of the sparse path.
+#[derive(Debug, Clone, Copy, Default, PartialEq)]
+pub struct TangentStats {
+    /// Element evaluations (parallel) + force summation; with the
+    /// tangents under full Newton, forces only on reused iterations.
+    pub assemblies: usize,
+    /// Extra tangent evaluations for a refresh decided after a
+    /// force-only assembly (slow contraction).
+    pub tangent_refreshes: usize,
+    /// Numeric factorisations.
+    pub factorizations: usize,
+    /// Linear solves (Newton corrections).
+    pub solves: usize,
+    /// Symbolic analyses (1 unless the pattern changed).
+    pub analyses: usize,
+    /// Solves that fell back to the banded LU (failed pivot).
+    pub fallbacks: usize,
+    /// Seconds in element evaluation + force summation.
+    pub assembly_seconds: f64,
+    /// Seconds in the tangent scatter (element blocks → CSR).
+    pub scatter_seconds: f64,
+    /// Seconds in the numeric factorisation.
+    pub factor_seconds: f64,
+    /// Seconds in triangular solves.
+    pub solve_seconds: f64,
+    /// Nonzeros of the factor.
+    pub nnz_l: usize,
+}
+
+/// State of the sparse path, owned by the stepper.
+pub(super) struct SparseTangent {
+    num_free: usize,
+    row_ptr: Vec<usize>,
+    col_idx: Vec<usize>,
+    values: Vec<f64>,
+    /// Per CSR entry, its element contributions in element order:
+    /// `contrib[contrib_ptr[p]..contrib_ptr[p + 1]]` = `(element,
+    /// column-major local index)` — the transposed scatter, so entries can
+    /// be summed in parallel in exactly the serial scatter's order.
+    contrib_ptr: Vec<usize>,
+    contrib: Vec<(u32, u32)>,
+    /// CSR index of each free DOF's diagonal.
+    diag_pos: Vec<usize>,
+    /// The element tangents of the latest assembly that computed them.
+    element_tangents: Vec<DMatrix<f64>>,
+    /// Whether `element_tangents` belong to the latest assembly's state.
+    tangents_current: bool,
+    solver: SparseLdlt,
+    /// Separate instance for the rest state's mass solve.
+    pub(super) mass_solver: SparseLdlt,
+    reuse: usize,
+    uses_left: usize,
+    factor_coef: f64,
+    factor_valid: bool,
+    fallback: BandedLu,
+    stats: TangentStats,
+}
+
+impl SparseTangent {
+    /// Build the fixed pattern (every free–free pair sharing an element)
+    /// and the element position maps; analyse it once.
+    pub(super) fn new<'e>(
+        elements: impl Iterator<Item = &'e [usize]> + Clone,
+        free_index: &[Option<usize>],
+        num_free: usize,
+        reuse: usize,
+    ) -> FeaResult<Self> {
+        let mut rows: Vec<Vec<usize>> = vec![Vec::new(); num_free];
+        for dofs in elements.clone() {
+            for &r in dofs {
+                let Some(fr) = free_index[r] else { continue };
+                for &c in dofs {
+                    if let Some(fc) = free_index[c] {
+                        rows[fr].push(fc);
+                    }
+                }
+            }
+        }
+        for (fr, row) in rows.iter_mut().enumerate() {
+            row.push(fr);
+            row.sort_unstable();
+            row.dedup();
+        }
+        let mut row_ptr = Vec::with_capacity(num_free + 1);
+        row_ptr.push(0);
+        let mut col_idx = Vec::new();
+        for row in &rows {
+            col_idx.extend_from_slice(row);
+            row_ptr.push(col_idx.len());
+        }
+        drop(rows);
+        let find = |r: usize, c: usize| -> usize {
+            let span = &col_idx[row_ptr[r]..row_ptr[r + 1]];
+            row_ptr[r]
+                + span
+                    .binary_search(&c)
+                    .expect("pattern covers every element pair")
+        };
+        // Element-major, row-major within the element: the order in
+        // which a serial scatter would add the contributions.
+        let mut flat: Vec<(usize, u32, u32)> = Vec::new();
+        for (e, dofs) in elements.enumerate() {
+            let n_e = dofs.len();
+            for (a, &r) in dofs.iter().enumerate() {
+                let Some(fr) = free_index[r] else { continue };
+                for (b, &c) in dofs.iter().enumerate() {
+                    if let Some(fc) = free_index[c] {
+                        flat.push((find(fr, fc), e as u32, (b * n_e + a) as u32));
+                    }
+                }
+            }
+        }
+        let mut contrib_ptr = vec![0usize; col_idx.len() + 1];
+        for &(p, _, _) in &flat {
+            contrib_ptr[p + 1] += 1;
+        }
+        for p in 0..col_idx.len() {
+            contrib_ptr[p + 1] += contrib_ptr[p];
+        }
+        let mut fill = contrib_ptr.clone();
+        let mut contrib = vec![(0u32, 0u32); flat.len()];
+        for (p, e, l) in flat {
+            contrib[fill[p]] = (e, l);
+            fill[p] += 1;
+        }
+        let diag_pos = (0..num_free).map(|d| find(d, d)).collect();
+        let mut solver = SparseLdlt::new();
+        solver.analyze(num_free, &row_ptr, &col_idx)?;
+        let values = vec![0.0; col_idx.len()];
+        let stats = TangentStats {
+            analyses: 1,
+            nnz_l: solver.stats().map_or(0, |s| s.nnz_l),
+            ..TangentStats::default()
+        };
+        Ok(Self {
+            num_free,
+            row_ptr,
+            col_idx,
+            values,
+            contrib_ptr,
+            contrib,
+            diag_pos,
+            element_tangents: Vec::new(),
+            tangents_current: false,
+            solver,
+            mass_solver: SparseLdlt::new(),
+            reuse: reuse.max(1),
+            uses_left: 0,
+            factor_coef: f64::NAN,
+            factor_valid: false,
+            fallback: BandedLu::new(),
+            stats,
+        })
+    }
+
+    pub(super) fn stats(&self) -> TangentStats {
+        self.stats
+    }
+
+    /// Whether the next Newton solve will factorise for sure (so the
+    /// assembly should produce the element tangents): always under full
+    /// Newton; under modified Newton when the factor is missing, used up,
+    /// for another dt, or `force_fresh`.
+    pub(super) fn wants_tangent(&self, force_fresh: bool, coef: f64) -> bool {
+        #[allow(clippy::float_cmp)]
+        let same_coef = coef == self.factor_coef;
+        self.reuse == 1 || force_fresh || !self.factor_valid || self.uses_left == 0 || !same_coef
+    }
+
+    /// Whether the solve at contraction `rate` will factorise.
+    pub(super) fn needs_refresh(&self, force_fresh: bool, coef: f64, rate: f64) -> bool {
+        self.wants_tangent(force_fresh, coef) || rate > 0.5
+    }
+
+    /// Whether the stored element tangents are those of the latest state.
+    pub(super) fn tangents_current(&self) -> bool {
+        self.tangents_current
+    }
+
+    /// Sum the element internal forces (element order, local order — the
+    /// banded path's order) into the free-DOF vector, and keep the element
+    /// tangents when they were computed.
+    pub(super) fn accumulate(
+        &mut self,
+        element_results: Vec<(DVector<f64>, Option<DMatrix<f64>>)>,
+        element_dofs: impl Iterator<Item = impl AsRef<[usize]>>,
+        free_index: &[Option<usize>],
+        started: Instant,
+    ) -> DVector<f64> {
+        let mut internal = DVector::zeros(self.num_free);
+        let with_tangents = element_results.first().is_some_and(|r| r.1.is_some());
+        if with_tangents {
+            self.element_tangents.clear();
+        }
+        for ((f_e, k_e), dofs) in element_results.into_iter().zip(element_dofs) {
+            for (local, &dof) in dofs.as_ref().iter().enumerate() {
+                if let Some(free) = free_index[dof] {
+                    internal[free] += f_e[local];
+                }
+            }
+            if let Some(k_e) = k_e {
+                self.element_tangents.push(k_e);
+            }
+        }
+        self.tangents_current = with_tangents;
+        self.stats.assemblies += 1;
+        self.stats.assembly_seconds += started.elapsed().as_secs_f64();
+        internal
+    }
+
+    /// Install element tangents computed separately for the latest state
+    /// (a modified-Newton refresh the force-only assembly did not foresee).
+    pub(super) fn set_tangents(&mut self, tangents: Vec<DMatrix<f64>>, started: Instant) {
+        self.element_tangents = tangents;
+        self.tangents_current = true;
+        self.stats.tangent_refreshes += 1;
+        self.stats.assembly_seconds += started.elapsed().as_secs_f64();
+    }
+
+    /// The Newton correction `(K_T + coef·(M + M_added)) δ = r`, with the
+    /// tangent of the latest [`Self::accumulate`]. `force_fresh` demands a
+    /// new factor; `rate` is `‖r_i‖/‖r_{i−1}‖` within the step (0 at the
+    /// first iteration).
+    pub(super) fn solve(
+        &mut self,
+        masses: &[&DMatrix<f64>],
+        added_mass: &DVector<f64>,
+        free_index: &[Option<usize>],
+        residual: &DVector<f64>,
+        coef: f64,
+        force_fresh: bool,
+        rate: f64,
+    ) -> FeaResult<DVector<f64>> {
+        if self.needs_refresh(force_fresh, coef, rate) {
+            if !self.tangents_current {
+                return Err(crate::error::SolverError::FactorizationFailed {
+                    reason: "sparse tangent: refresh without the current element tangents"
+                        .to_string(),
+                }
+                .into());
+            }
+            let started = Instant::now();
+            let tangents = &self.element_tangents;
+            let contrib_ptr = &self.contrib_ptr;
+            let contrib = &self.contrib;
+            self.values
+                .par_chunks_mut(SCATTER_CHUNK)
+                .enumerate()
+                .for_each(|(chunk, values)| {
+                    let base = chunk * SCATTER_CHUNK;
+                    for (offset, value) in values.iter_mut().enumerate() {
+                        let p = base + offset;
+                        let mut sum = 0.0;
+                        for &(e, l) in &contrib[contrib_ptr[p]..contrib_ptr[p + 1]] {
+                            let (e, l) = (e as usize, l as usize);
+                            sum += tangents[e].as_slice()[l] + coef * masses[e].as_slice()[l];
+                        }
+                        *value = sum;
+                    }
+                });
+            for (dof, &m) in added_mass.iter().enumerate() {
+                if m != 0.0 {
+                    if let Some(free) = free_index[dof] {
+                        self.values[self.diag_pos[free]] += coef * m;
+                    }
+                }
+            }
+            let scattered = Instant::now();
+            self.stats.scatter_seconds += (scattered - started).as_secs_f64();
+            let outcome = self.solver.factorize_values(&self.values);
+            self.stats.factor_seconds += scattered.elapsed().as_secs_f64();
+            if let Err(error) = outcome {
+                // A failed pivot: this solve goes to the pivoting banded
+                // LU on the same matrix; the next one tries LDLᵀ again.
+                if self.stats.fallbacks == 0 {
+                    tracing::warn!("sparse tangent: {error}; falling back to the banded LU");
+                }
+                self.stats.fallbacks += 1;
+                self.factor_valid = false;
+                let mut triplets = Vec::with_capacity(self.values.len());
+                for row in 0..self.num_free {
+                    for idx in self.row_ptr[row]..self.row_ptr[row + 1] {
+                        triplets.push((row, self.col_idx[idx], self.values[idx]));
+                    }
+                }
+                let matrix = SparseMatrix::from_triplets(self.num_free, self.num_free, &triplets)?;
+                let (x, _) = self
+                    .fallback
+                    .solve(&matrix, residual, &SolverOptions::default())?;
+                self.stats.solves += 1;
+                return Ok(x);
+            }
+            self.stats.factorizations += 1;
+            self.factor_valid = true;
+            self.factor_coef = coef;
+            self.uses_left = self.reuse;
+        }
+        let started = Instant::now();
+        let x = self.solver.solve_factored(residual)?;
+        self.stats.solve_seconds += started.elapsed().as_secs_f64();
+        self.stats.solves += 1;
+        self.uses_left -= 1;
+        Ok(x)
+    }
+}

crates/specialized/rtx-fea/src/elements/total_lagrangian.rs

@@ -194,3 +194,96 @@ pub fn internal_force_and_tangent(
 
     Ok((internal_force, tangent))
 }
+
+/// The internal force alone of [`internal_force_and_tangent`] — the same
+/// operations in the same order for `f_int` (so the result is identical
+/// bit for bit), without the tangent's `B_Lᵀ C B_L` and geometric terms.
+/// For Newton iterations that reuse a factorised tangent (modified
+/// Newton, R8-g).
+pub fn internal_force(
+    element: &dyn FiniteElement,
+    node_coords: &[Vector3<f64>],
+    element_displacement: &DVector<f64>,
+    constitutive: &dyn Fn(&DVector<f64>) -> FeaResult<(DVector<f64>, DMatrix<f64>)>,
+    quadrature_order: Option<usize>,
+) -> FeaResult<DVector<f64>> {
+    let quad_rule = element.quadrature_rule(quadrature_order)?;
+    let dim = element.spatial_dimension();
+    let num_nodes = element.num_nodes();
+    let total_dofs = num_nodes * dim;
+    if element_displacement.len() != total_dofs {
+        return Err(ElementError::MatrixComputationFailed {
+            reason: format!(
+                "element displacement has {} entries, element has {total_dofs} DOFs",
+                element_displacement.len()
+            ),
+        }
+        .into());
+    }
+    if dim != 2 && dim != 3 {
+        return Err(ElementError::MatrixComputationFailed {
+            reason: format!("total-Lagrangian formulation needs 2-D or 3-D, got {dim}"),
+        }
+        .into());
+    }
+    let n_strain = if dim == 2 { 3 } else { 6 };
+    let pairs: Vec<(usize, usize)> = if dim == 2 {
+        vec![(0, 0), (1, 1), (0, 1)]
+    } else {
+        vec![(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
+    };
+
+    let mut internal_force = DVector::zeros(total_dofs);
+    for point in &quad_rule.points {
+        let shape_eval = element.shape_functions(&point.coords)?;
+        let jacobian_eval = element.jacobian(&point.coords, node_coords)?;
+        if !jacobian_eval.is_valid() {
+            return Err(ElementError::JacobianSingular {
+                det: jacobian_eval.determinant,
+            }
+            .into());
+        }
+        let g = jacobian_eval.transform_derivatives(&shape_eval.derivatives)?;
+        let mut f: DMatrix<f64> = DMatrix::identity(dim, dim);
+        for a in 0..num_nodes {
+            for i in 0..dim {
+                let ui = element_displacement[a * dim + i];
+                for j in 0..dim {
+                    f[(i, j)] += ui * g[(a, j)];
+                }
+            }
+        }
+        let c_right = f.transpose() * &f;
+        let mut e_voigt: DVector<f64> = DVector::zeros(n_strain);
+        for (k, &(i, j)) in pairs.iter().enumerate() {
+            let e_ij = 0.5 * (c_right[(i, j)] - if i == j { 1.0 } else { 0.0 });
+            e_voigt[k] = if i == j { e_ij } else { 2.0 * e_ij };
+        }
+        let (s_voigt, _) = constitutive(&e_voigt)?;
+        if s_voigt.len() != n_strain {
+            return Err(ElementError::MatrixComputationFailed {
+                reason: format!(
+                    "constitutive closure returned stress of {} components, expected {n_strain}",
+                    s_voigt.len()
+                ),
+            }
+            .into());
+        }
+        let mut b_l: DMatrix<f64> = DMatrix::zeros(n_strain, total_dofs);
+        for (k, &(i, j)) in pairs.iter().enumerate() {
+            for a in 0..num_nodes {
+                for m in 0..dim {
+                    let col = a * dim + m;
+                    b_l[(k, col)] = if i == j {
+                        f[(m, i)] * g[(a, i)]
+                    } else {
+                        f[(m, i)] * g[(a, j)] + f[(m, j)] * g[(a, i)]
+                    };
+                }
+            }
+        }
+        let scale = jacobian_eval.determinant().abs() * point.weight;
+        internal_force += b_l.transpose() * &s_voigt * scale;
+    }
+    Ok(internal_force)
+}

crates/specialized/rtx-fea/src/solvers/mod.rs

@@ -15,6 +15,7 @@ pub mod iterative;
 #[cfg(all(target_os = "macos", feature = "metal"))]
 pub mod metal_solvers;
 pub mod nonlinear;
+pub mod sparse_ldlt;
 pub mod time_integration;
 
 #[cfg(disabled)]
@@ -34,6 +35,7 @@ pub use iterative::*;
 #[cfg(all(target_os = "macos", feature = "metal"))]
 pub use metal_solvers::*;
 pub use nonlinear::*;
+pub use sparse_ldlt::{SparseLdlt, SparseLdltStats};
 pub use time_integration::*;
 
 /// Solver configuration and options.

crates/specialized/rtx-fea/src/solvers/sparse_ldlt/mod.rs

@@ -0,0 +1,369 @@
+// Copyright (c) 2024 RustyTorch++ Team
+// Licensed under the Apache License, Version 2.0
+
+//! Sparse supernodal LDLᵀ direct solver with symbolic reuse.
+//!
+//! Built for the Newton tangent of the 3-D flag (R8-g in omni-cortex):
+//! `K_T + M/(β Δt²)` is symmetric, its sparsity pattern is fixed by the
+//! mesh, and the banded LU it replaces does `O(n·b²)` work on a Hex20
+//! bandwidth of one x-slab (~500 DOFs at 35×2×8). Here:
+//!
+//! * [`SparseLdlt::analyze`] — nested-dissection ordering, elimination
+//!   tree, supernodes and every scatter map, once per pattern;
+//! * [`SparseLdlt::factorize_csr`] — the numeric multifrontal
+//!   factorisation on the analysed pattern (re-analysing only when the
+//!   pattern changes), parallel over independent subtrees and
+//!   bit-deterministic at any thread count;
+//! * [`SparseLdlt::solve_factored`] — the two triangular solves.
+//!
+//! No pivoting: the matrix must be symmetric and have nonzero pivots in
+//! the chosen order (true of an SPD matrix, and of the mildly indefinite
+//! Newmark tangents met in practice); a failed pivot is an error the
+//! caller can answer with a pivoting solver. Only the lower triangle is
+//! read — the upper is assumed equal.
+
+mod numeric;
+mod ordering;
+mod symbolic;
+
+use super::{ConvergenceInfo, LinearSolver, SolverCapabilities, SolverOptions};
+use crate::assembly::SparseMatrix;
+use crate::error::{FeaResult, SolverError};
+use nalgebra::DVector;
+use std::time::Instant;
+
+/// Pattern and cost statistics of the current analysis.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct SparseLdltStats {
+    /// Matrix dimension.
+    pub n: usize,
+    /// Supernodes in the assembly tree.
+    pub supernodes: usize,
+    /// Nonzeros of L including the diagonal.
+    pub nnz_l: usize,
+    /// Largest front dimension.
+    pub max_front: usize,
+    /// Estimated multiply–adds of one numeric factorisation.
+    pub work: f64,
+    /// Symbolic analyses performed so far.
+    pub analyses: usize,
+    /// Numeric factorisations performed so far.
+    pub factorizations: usize,
+}
+
+/// Sparse LDLᵀ with symbolic reuse. See the module docs.
+#[derive(Debug, Default)]
+pub struct SparseLdlt {
+    symbolic: Option<symbolic::Symbolic>,
+    factor: numeric::Factor,
+    factorized: bool,
+    analyses: usize,
+    factorizations: usize,
+    work: Vec<f64>,
+}
+
+impl SparseLdlt {
+    /// A solver with no analysis yet.
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    /// Symbolic analysis of an `n × n` CSR pattern (structurally
+    /// symmetric; sorted column indices per row are not required).
+    pub fn analyze(&mut self, n: usize, row_ptr: &[usize], col_idx: &[usize]) -> FeaResult<()> {
+        if row_ptr.len() != n + 1 || row_ptr[n] != col_idx.len() {
+            return Err(SolverError::FactorizationFailed {
+                reason: format!(
+                    "sparse LDLt: CSR pattern inconsistent with n = {n} ({} row pointers, {} columns)",
+                    row_ptr.len(),
+                    col_idx.len()
+                ),
+            }
+            .into());
+        }
+        if col_idx.len() > u32::MAX as usize {
+            return Err(SolverError::FactorizationFailed {
+                reason: "sparse LDLt: more than 2^32 nonzeros".to_string(),
+            }
+            .into());
+        }
+        let sym = symbolic::Symbolic::analyze(n, row_ptr, col_idx);
+        self.factor.resize(&sym);
+        self.symbolic = Some(sym);
+        self.factorized = false;
+        self.analyses += 1;
+        Ok(())
+    }
+
+    /// Numeric factorisation of the CSR matrix `(row_ptr, col_idx,
+    /// values)`. Reuses the symbolic analysis when the pattern is the
+    /// analysed one, re-analyses otherwise.
+    pub fn factorize_csr(
+        &mut self,
+        row_ptr: &[usize],
+        col_idx: &[usize],
+        values: &[f64],
+    ) -> FeaResult<()> {
+        let n = row_ptr.len().saturating_sub(1);
+        let reuse = self
+            .symbolic
+            .as_ref()
+            .is_some_and(|s| s.n == n && s.same_pattern(row_ptr, col_idx));
+        if !reuse {
+            self.analyze(n, row_ptr, col_idx)?;
+        }
+        self.factorize_values(values)
+    }
+
+    /// Numeric factorisation of new values on the analysed pattern
+    /// (the caller guarantees the pattern is unchanged).
+    pub fn factorize_values(&mut self, values: &[f64]) -> FeaResult<()> {
+        let sym = self
+            .symbolic
+            .as_ref()
+            .ok_or(SolverError::FactorizationRequired)?;
+        if values.len() != sym.pattern_col_idx.len() {
+            return Err(SolverError::FactorizationFailed {
+                reason: "sparse LDLt: value count differs from the analysed pattern".to_string(),
+            }
+            .into());
+        }
+        self.factorized = false;
+        // Pivot floor relative to the largest diagonal entry.
+        let mut diag_max = 0.0f64;
+        for row in 0..sym.n {
+            for idx in sym.pattern_row_ptr[row]..sym.pattern_row_ptr[row + 1] {
+                if sym.pattern_col_idx[idx] == row {
+                    diag_max = diag_max.max(values[idx].abs());
+                }
+            }
+        }
+        let floor = 1e-14 * diag_max;
+        numeric::factorize(sym, values, &self.factor, floor).map_err(|fail| {
+            SolverError::FactorizationFailed {
+                reason: format!(
+                    "sparse LDLt: pivot {:.3e} at factor column {} (floor {floor:.3e})",
+                    fail.value, fail.column
+                ),
+            }
+        })?;
+        self.factorized = true;
+        self.factorizations += 1;
+        Ok(())
+    }
+
+    /// Solve with the current factors.
+    pub fn solve_factored(&mut self, rhs: &DVector<f64>) -> FeaResult<DVector<f64>> {
+        let sym = self
+            .symbolic
+            .as_ref()
+            .ok_or(SolverError::FactorizationRequired)?;
+        if !self.factorized {
+            return Err(SolverError::FactorizationRequired.into());
+        }
+        if rhs.len() != sym.n {
+            return Err(SolverError::DimensionMismatch {
+                matrix_rows: sym.n,
+                matrix_cols: sym.n,
+                rhs_rows: rhs.len(),
+                rhs_cols: 1,
+            }
+            .into());
+        }
+        self.work.clear();
+        self.work.extend(sym.perm.iter().map(|&old| rhs[old]));
+        numeric::solve_in_place(sym, &self.factor, &mut self.work);
+        let mut x = DVector::zeros(sym.n);
+        for (new, &old) in sym.perm.iter().enumerate() {
+            x[old] = self.work[new];
+        }
+        Ok(x)
+    }
+
+    /// Whether a numeric factor is available.
+    pub fn is_factorized(&self) -> bool {
+        self.factorized
+    }
+
+    /// Statistics of the current analysis (`None` before the first).
+    pub fn stats(&self) -> Option<SparseLdltStats> {
+        self.symbolic.as_ref().map(|s| SparseLdltStats {
+            n: s.n,
+            supernodes: s.num_supernodes(),
+            nnz_l: s.nnz_l,
+            max_front: s.max_front(),
+            work: s.work,
+            analyses: self.analyses,
+            factorizations: self.factorizations,
+        })
+    }
+}
+
+impl LinearSolver for SparseLdlt {
+    fn solve(
+        &mut self,
+        matrix: &SparseMatrix,
+        rhs: &DVector<f64>,
+        _options: &SolverOptions,
+    ) -> FeaResult<(DVector<f64>, ConvergenceInfo)> {
+        let start_time = Instant::now();
+        if matrix.nrows() != matrix.ncols() || matrix.nrows() != rhs.len() {
+            return Err(SolverError::DimensionMismatch {
+                matrix_rows: matrix.nrows(),
+                matrix_cols: matrix.ncols(),
+                rhs_rows: rhs.len(),
+                rhs_cols: 1,
+            }
+            .into());
+        }
+        let (row_ptr, col_idx) = matrix.structure();
+        if row_ptr.len() != matrix.nrows() + 1 {
+            return Err(SolverError::FactorizationFailed {
+                reason: "sparse LDLt needs a finalized (CSR) matrix".to_string(),
+            }
+            .into());
+        }
+        self.factorize_csr(row_ptr, col_idx, matrix.values())?;
+        let x = self.solve_factored(rhs)?;
+        let mut info = ConvergenceInfo::new();
+        info.set_solve_time(start_time.elapsed());
+        info.set_converged(1, 0.0, 0.0);
+        Ok((x, info))
+    }
+
+    fn name(&self) -> &'static str {
+        "Sparse supernodal LDLt"
+    }
+
+    fn capabilities(&self) -> SolverCapabilities {
+        SolverCapabilities {
+            symmetric: true,
+            positive_definite: false,
+            gpu_acceleration: false,
+            multiple_rhs: true,
+            iterative_refinement: false,
+            memory_efficiency: 5,
+            computational_efficiency: 5,
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::solvers::BandedLu;
+
+    /// A 3-D 7-point Laplacian with `dofs` components per node coupled
+    /// by a small symmetric block, shifted by `shift` (negative shifts
+    /// make it indefinite).
+    fn laplacian_3d(nx: usize, ny: usize, nz: usize, dofs: usize, shift: f64) -> SparseMatrix {
+        let nn = nx * ny * nz;
+        let n = nn * dofs;
+        let id = |i: usize, j: usize, k: usize| (k * ny + j) * nx + i;
+        let mut m = SparseMatrix::new(n, n);
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    let p = id(i, j, k);
+                    let mut nbrs = Vec::new();
+                    if i > 0 {
+                        nbrs.push(id(i - 1, j, k));
+                    }
+                    if i + 1 < nx {
+                        nbrs.push(id(i + 1, j, k));
+                    }
+                    if j > 0 {
+                        nbrs.push(id(i, j - 1, k));
+                    }
+                    if j + 1 < ny {
+                        nbrs.push(id(i, j + 1, k));
+                    }
+                    if k > 0 {
+                        nbrs.push(id(i, j, k - 1));
+                    }
+                    if k + 1 < nz {
+                        nbrs.push(id(i, j, k + 1));
+                    }
+                    for a in 0..dofs {
+                        let row = p * dofs + a;
+                        m.add_entry(row, row, 6.0 + shift + 0.1 * a as f64).unwrap();
+                        for b in 0..dofs {
+                            if a != b {
+                                m.add_entry(row, p * dofs + b, 0.3).unwrap();
+                            }
+                        }
+                        for &q in &nbrs {
+                            m.add_entry(row, q * dofs + a, -1.0).unwrap();
+                            for b in 0..dofs {
+                                if a != b {
+                                    m.add_entry(row, q * dofs + b, 0.05).unwrap();
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        m.finalize().unwrap();
+        m
+    }
+
+    fn check(m: &SparseMatrix) {
+        let n = m.nrows();
+        let x_exact = DVector::from_fn(n, |i, _| ((i as f64) * 0.37).sin() + 0.2);
+        let b = m.multiply_vector(&x_exact).unwrap();
+        let opts = SolverOptions::default();
+        let (x, _) = SparseLdlt::new().solve(m, &b, &opts).unwrap();
+        let rel = (&x - &x_exact).norm() / x_exact.norm();
+        assert!(rel < 1e-11, "manufactured rel err {rel:.3e}");
+        let (xb, _) = BandedLu::new().solve(m, &b, &opts).unwrap();
+        let rel_b = (&x - &xb).norm() / xb.norm();
+        assert!(rel_b < 1e-11, "vs banded LU rel err {rel_b:.3e}");
+    }
+
+    #[test]
+    fn spd_3d_block_laplacian() {
+        check(&laplacian_3d(12, 3, 5, 3, 0.0));
+    }
+
+    #[test]
+    fn indefinite_shifted_laplacian() {
+        check(&laplacian_3d(9, 4, 4, 2, -5.3));
+    }
+
+    #[test]
+    fn tiny_and_diagonal() {
+        check(&laplacian_3d(1, 1, 1, 1, 0.0));
+        check(&laplacian_3d(3, 1, 1, 1, 0.0));
+    }
+
+    #[test]
+    fn symbolic_reuse_across_values() {
+        let m1 = laplacian_3d(10, 3, 4, 3, 0.0);
+        let m2 = laplacian_3d(10, 3, 4, 3, 1.5);
+        let mut solver = SparseLdlt::new();
+        let opts = SolverOptions::default();
+        for m in [&m1, &m2, &m1] {
+            let x_exact = DVector::from_fn(m.nrows(), |i, _| (i as f64 * 0.11).cos());
+            let b = m.multiply_vector(&x_exact).unwrap();
+            let (x, _) = solver.solve(m, &b, &opts).unwrap();
+            assert!((&x - &x_exact).norm() / x_exact.norm() < 1e-11);
+        }
+        let stats = solver.stats().unwrap();
+        assert_eq!(stats.analyses, 1);
+        assert_eq!(stats.factorizations, 3);
+    }
+
+    #[test]
+    fn zero_pivot_is_an_error() {
+        let mut m = SparseMatrix::new(2, 2);
+        m.add_entry(0, 1, 1.0).unwrap();
+        m.add_entry(1, 0, 1.0).unwrap();
+        m.add_entry(1, 1, 1.0).unwrap();
+        m.finalize().unwrap();
+        let b = DVector::from_vec(vec![1.0, 2.0]);
+        // The (0,0) pivot is structurally absent: zero.
+        let result = SparseLdlt::new().solve(&m, &b, &SolverOptions::default());
+        assert!(result.is_err());
+    }
+}

crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs

@@ -0,0 +1,373 @@
+//! Numeric phase: multifrontal supernodal LDLᵀ (no pivoting) over the
+//! symbolic assembly tree, and the triangular solves.
+//!
+//! Each supernode's front is a dense column-major `m × m` block (lower
+//! triangle used): the original entries are scattered in, the children's
+//! update matrices extend-added (in child order), then the `k` pivot
+//! columns are factorised by a blocked right-looking LDLᵀ whose trailing
+//! updates are `matrixmultiply` gemms (through nalgebra views). Sibling
+//! subtrees run on rayon when their estimated work is large; the
+//! arithmetic of every supernode is independent of the thread schedule,
+//! so results are bit-identical at any thread count.
+
+use super::symbolic::Symbolic;
+use nalgebra::{DMatrixView, DMatrixViewMut};
+use rayon::prelude::*;
+use std::sync::Mutex;
+
+/// Panel width of the blocked partial factorisation.
+const PANEL: usize = 32;
+/// Column block of the lower-triangular trailing update.
+const TRAIL_BLOCK: usize = 64;
+/// Subtrees below this many multiply–adds are factorised serially.
+const PARALLEL_WORK: f64 = 2.0e6;
+/// Trailing updates below this many multiply–adds (×2) run serially.
+const PARALLEL_FRONT_WORK: f64 = 4.0e6;
+
+/// A zero, non-finite or relatively tiny pivot.
+#[derive(Debug, Clone)]
+pub(super) struct PivotFailure {
+    pub column: usize,
+    pub value: f64,
+}
+
+/// Per-supernode factor blocks: `m × k` column-major, unit-lower L with
+/// D on the diagonal. Allocations persist across factorisations.
+#[derive(Debug, Default)]
+pub(super) struct Factor {
+    pub blocks: Vec<Mutex<Vec<f64>>>,
+}
+
+impl Factor {
+    pub fn resize(&mut self, sym: &Symbolic) {
+        let ns = sym.num_supernodes();
+        if self.blocks.len() != ns {
+            self.blocks = (0..ns).map(|_| Mutex::new(Vec::new())).collect();
+        }
+    }
+}
+
+/// Factorise the matrix whose CSR values are `values` (on the analysed
+/// pattern). `pivot_floor` rejects |d| below it.
+pub(super) fn factorize(
+    sym: &Symbolic,
+    values: &[f64],
+    factor: &Factor,
+    pivot_floor: f64,
+) -> Result<(), PivotFailure> {
+    let ctx = Ctx {
+        sym,
+        values,
+        factor,
+        pivot_floor,
+    };
+    let results: Vec<Result<Vec<f64>, PivotFailure>> = if sym.roots.len() > 1 {
+        sym.roots.par_iter().map(|&r| ctx.subtree(r)).collect()
+    } else {
+        sym.roots.iter().map(|&r| ctx.subtree(r)).collect()
+    };
+    for r in results {
+        r?;
+    }
+    Ok(())
+}
+
+struct Ctx<'a> {
+    sym: &'a Symbolic,
+    values: &'a [f64],
+    factor: &'a Factor,
+    pivot_floor: f64,
+}
+
+impl Ctx<'_> {
+    /// Factorise the subtree rooted at `s`; returns its update matrix
+    /// (`(m − k)²`, column-major).
+    fn subtree(&self, s: usize) -> Result<Vec<f64>, PivotFailure> {
+        let sym = self.sym;
+        let kids = sym.node_children(s);
+        let updates: Vec<Result<Vec<f64>, PivotFailure>> =
+            if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_WORK {
+                kids.par_iter().map(|&c| self.subtree(c)).collect()
+            } else {
+                kids.iter().map(|&c| self.subtree(c)).collect()
+            };
+
+        let rows = sym.front_rows(s);
+        let m = rows.len();
+        let f = sym.start[s];
+        let k = sym.start[s + 1] - f;
+        let mut front = vec![0.0f64; m * m];
+        for &(idx, off) in &sym.amap[sym.amap_ptr[s]..sym.amap_ptr[s + 1]] {
+            front[off as usize] += self.values[idx as usize];
+        }
+        for (&c, update) in kids.iter().zip(updates) {
+            let update = update?;
+            let rel = &sym.relind[sym.relind_ptr[c]..sym.relind_ptr[c + 1]];
+            let u = rel.len();
+            for b in 0..u {
+                let col = rel[b] as usize * m;
+                let src = &update[b * u..(b + 1) * u];
+                for a in b..u {
+                    front[rel[a] as usize + col] += src[a];
+                }
+            }
+        }
+        partial_ldlt(&mut front, m, k, self.pivot_floor).map_err(|(j, value)| PivotFailure {
+            column: f + j,
+            value,
+        })?;
+        {
+            let mut block = self.factor.blocks[s].lock().unwrap();
+            block.clear();
+            block.extend_from_slice(&front[..m * k]);
+        }
+        let u = m - k;
+        let mut update = vec![0.0f64; u * u];
+        for b in 0..u {
+            let src = (k + b) * m + k;
+            update[b * u + b..(b + 1) * u].copy_from_slice(&front[src + b..src + u]);
+        }
+        Ok(update)
+    }
+}
+
+/// Blocked right-looking LDLᵀ of the first `k` columns of the `m × m`
+/// column-major front, updating the trailing `(m−k)²` lower block.
+/// On exit columns `0..k` hold L (unit, strictly below the diagonal) and
+/// D (on the diagonal).
+fn partial_ldlt(front: &mut [f64], m: usize, k: usize, floor: f64) -> Result<(), (usize, f64)> {
+    let mut w = Vec::new();
+    let mut p0 = 0;
+    while p0 < k {
+        let p1 = (p0 + PANEL).min(k);
+        // Unblocked factorisation of the panel's columns over rows p0..m.
+        for j in p0..p1 {
+            let d = front[j + j * m];
+            if !d.is_finite() || d.abs() <= floor {
+                return Err((j, d));
+            }
+            let inv_d = 1.0 / d;
+            for c in (j + 1)..p1 {
+                let lcj_d = front[c + j * m] * inv_d;
+                if lcj_d != 0.0 {
+                    let (left, right) = front.split_at_mut(c * m);
+                    let src = &left[j * m + c..j * m + m];
+                    let dst = &mut right[c..m];
+                    for (x, &y) in dst.iter_mut().zip(src) {
+                        *x -= lcj_d * y;
+                    }
+                }
+            }
+            for x in &mut front[j * m + j + 1..j * m + m] {
+                *x *= inv_d;
+            }
+        }
+        // Trailing update of columns p1..m (lower): F -= L_P D_P L_Pᵀ.
+        if p1 < m {
+            let nb = p1 - p0;
+            let r = m - p1;
+            // W = L_P · D_P, rows p1..m, column-major r × nb.
+            w.clear();
+            w.resize(r * nb, 0.0);
+            for (t, j) in (p0..p1).enumerate() {
+                let d = front[j + j * m];
+                for (x, &l) in w[t * r..(t + 1) * r]
+                    .iter_mut()
+                    .zip(&front[j * m + p1..j * m + m])
+                {
+                    *x = l * d;
+                }
+            }
+            let (left, right) = front.split_at_mut(p1 * m);
+            let left: &[f64] = left;
+            let w: &[f64] = &w;
+            // Each column block of the trailing matrix is an independent
+            // gemm; large fronts run them on rayon (the arithmetic of a
+            // block does not depend on the schedule).
+            let block_update = |(t, chunk): (usize, &mut [f64])| {
+                let c0 = t * TRAIL_BLOCK;
+                let c1 = (c0 + TRAIL_BLOCK).min(r);
+                let rows = r - c0;
+                // L_P rows (p1 + c0)..m: a view into `left` (columns p0..p1).
+                let lp = DMatrixView::from_slice_with_strides_generic(
+                    &left[p0 * m + p1 + c0..],
+                    nalgebra::Dyn(rows),
+                    nalgebra::Dyn(nb),
+                    nalgebra::Dyn(1),
+                    nalgebra::Dyn(m),
+                );
+                // W rows c0..c1 (the columns being updated), transposed.
+                let wb = DMatrixView::from_slice_with_strides_generic(
+                    &w[c0..],
+                    nalgebra::Dyn(c1 - c0),
+                    nalgebra::Dyn(nb),
+                    nalgebra::Dyn(1),
+                    nalgebra::Dyn(r),
+                );
+                let mut target = DMatrixViewMut::from_slice_with_strides_generic(
+                    &mut chunk[p1 + c0..],
+                    nalgebra::Dyn(rows),
+                    nalgebra::Dyn(c1 - c0),
+                    nalgebra::Dyn(1),
+                    nalgebra::Dyn(m),
+                );
+                target.gemm(-1.0, &lp, &wb.transpose(), 1.0);
+            };
+            if (r * r * nb) as f64 > PARALLEL_FRONT_WORK {
+                right
+                    .par_chunks_mut(TRAIL_BLOCK * m)
+                    .enumerate()
+                    .for_each(block_update);
+            } else {
+                right
+                    .chunks_mut(TRAIL_BLOCK * m)
+                    .enumerate()
+                    .for_each(block_update);
+            }
+        }
+        p0 = p1;
+    }
+    Ok(())
+}
+
+/// Subtrees below this factorisation work solve serially.
+const PARALLEL_SOLVE_WORK: f64 = 2.0e7;
+
+/// Solve `L D Lᵀ x = b` in the factor numbering, in place, over the
+/// assembly tree: the forward solve is multifrontal (each supernode
+/// solves its dense unit-triangular block and hands its parent an update
+/// vector, children in order), the backward solve top-down (each
+/// supernode gathers its off-diagonal unknowns from its parent's front
+/// vector); independent subtrees run on rayon, and every supernode's
+/// arithmetic is independent of the schedule.
+pub(super) fn solve_in_place(sym: &Symbolic, factor: &Factor, x: &mut [f64]) {
+    let blocks: Vec<std::sync::MutexGuard<'_, Vec<f64>>> =
+        factor.blocks.iter().map(|b| b.lock().unwrap()).collect();
+    let blocks: Vec<&[f64]> = blocks.iter().map(|b| b.as_slice()).collect();
+    // One disjoint output slice per supernode (its own columns).
+    let mut slots: Vec<Mutex<&mut [f64]>> = Vec::with_capacity(sym.num_supernodes());
+    {
+        let mut rest: &mut [f64] = x;
+        for s in 0..sym.num_supernodes() {
+            let k = sym.start[s + 1] - sym.start[s];
+            let (own, tail) = rest.split_at_mut(k);
+            slots.push(Mutex::new(own));
+            rest = tail;
+        }
+    }
+    let ctx = SolveCtx {
+        sym,
+        blocks: &blocks,
+        slots: &slots,
+    };
+    let roots = &sym.roots;
+    let run_forward = |&r: &usize| {
+        ctx.forward(r);
+    };
+    if roots.len() > 1 {
+        roots.par_iter().for_each(run_forward);
+    } else {
+        roots.iter().for_each(run_forward);
+    }
+    let run_backward = |&r: &usize| ctx.backward(r, &[]);
+    if roots.len() > 1 {
+        roots.par_iter().for_each(run_backward);
+    } else {
+        roots.iter().for_each(run_backward);
+    }
+}
+
+struct SolveCtx<'a, 'x> {
+    sym: &'a Symbolic,
+    blocks: &'a [&'a [f64]],
+    slots: &'a [Mutex<&'x mut [f64]>],
+}
+
+impl SolveCtx<'_, '_> {
+    /// Forward (`L y = b`, then `z = D⁻¹ y` on the own columns) for the
+    /// subtree at `s`; returns the update vector for the parent's front.
+    fn forward(&self, s: usize) -> Vec<f64> {
+        let sym = self.sym;
+        let kids = sym.node_children(s);
+        let updates: Vec<Vec<f64>> = if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_SOLVE_WORK
+        {
+            kids.par_iter().map(|&c| self.forward(c)).collect()
+        } else {
+            kids.iter().map(|&c| self.forward(c)).collect()
+        };
+        let m = sym.rows_ptr[s + 1] - sym.rows_ptr[s];
+        let k = sym.start[s + 1] - sym.start[s];
+        let block = self.blocks[s];
+        let mut v = vec![0.0f64; m];
+        {
+            let own = self.slots[s].lock().unwrap();
+            v[..k].copy_from_slice(&own);
+        }
+        for (&c, update) in kids.iter().zip(&updates) {
+            let rel = &sym.relind[sym.relind_ptr[c]..sym.relind_ptr[c + 1]];
+            for (&r, &u) in rel.iter().zip(update) {
+                v[r as usize] += u;
+            }
+        }
+        let (head, tail) = v.split_at_mut(k);
+        for j in 0..k {
+            let xj = head[j];
+            if xj != 0.0 {
+                let col = &block[j * m..(j + 1) * m];
+                for (y, &l) in head[j + 1..].iter_mut().zip(&col[j + 1..k]) {
+                    *y -= l * xj;
+                }
+                for (y, &l) in tail.iter_mut().zip(&col[k..]) {
+                    *y -= l * xj;
+                }
+            }
+        }
+        {
+            let mut own = self.slots[s].lock().unwrap();
+            for j in 0..k {
+                own[j] = head[j] / block[j + j * m];
+            }
+        }
+        v.drain(..k);
+        v
+    }
+
+    /// Backward (`Lᵀ x = z`) for the subtree at `s`, given the parent's
+    /// front vector (the solution at the parent's front rows).
+    fn backward(&self, s: usize, parent_front: &[f64]) {
+        let sym = self.sym;
+        let m = sym.rows_ptr[s + 1] - sym.rows_ptr[s];
+        let k = sym.start[s + 1] - sym.start[s];
+        let block = self.blocks[s];
+        let mut front = vec![0.0f64; m];
+        let rel = &sym.relind[sym.relind_ptr[s]..sym.relind_ptr[s + 1]];
+        for (slot, &r) in front[k..].iter_mut().zip(rel) {
+            *slot = parent_front[r as usize];
+        }
+        {
+            let mut own = self.slots[s].lock().unwrap();
+            let (head, tail) = front.split_at_mut(k);
+            head.copy_from_slice(&own);
+            for j in (0..k).rev() {
+                let col = &block[j * m..(j + 1) * m];
+                let mut acc = 0.0;
+                for (&l, &y) in col[k..].iter().zip(tail.iter()) {
+                    acc += l * y;
+                }
+                for i in (j + 1)..k {
+                    acc += col[i] * head[i];
+                }
+                head[j] -= acc;
+            }
+            own.copy_from_slice(head);
+        }
+        let kids = sym.node_children(s);
+        if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_SOLVE_WORK {
+            kids.par_iter().for_each(|&c| self.backward(c, &front));
+        } else {
+            for &c in kids {
+                self.backward(c, &front);
+            }
+        }
+    }
+}

crates/specialized/rtx-fea/src/solvers/sparse_ldlt/ordering.rs

@@ -0,0 +1,336 @@
+//! Fill-reducing ordering: nested dissection on the matrix graph.
+//!
+//! The graph is first compressed to supervariables (DOFs with identical
+//! closed neighbourhoods — the 2 or 3 displacement components of one FE
+//! node), then dissected recursively: a pseudo-peripheral vertex's BFS
+//! level structure supplies the candidate separators (one level each,
+//! trimmed to the vertices that actually touch the far side), the
+//! lightest one that keeps both halves at least 30 % of the weight is
+//! taken, and the ordering is `A, B, S`. On a long, thin FE body (the
+//! flag) the level sets are cross-sections, so the top separators are the
+//! minimal ones; small subgraphs are left in natural order.
+//!
+//! Deterministic: no hashing order, no randomness, no threads.
+
+use std::collections::HashMap;
+
+/// Subgraphs at or below this weight (DOFs) are not dissected further.
+const LEAF_WEIGHT: usize = 96;
+
+/// Compressed adjacency: vertex `v`'s neighbours are
+/// `adj[ptr[v]..ptr[v + 1]]` (no self loops), its weight the number of
+/// DOFs it stands for.
+struct Graph {
+    ptr: Vec<usize>,
+    adj: Vec<u32>,
+    weight: Vec<usize>,
+}
+
+impl Graph {
+    fn neighbours(&self, v: u32) -> &[u32] {
+        &self.adj[self.ptr[v as usize]..self.ptr[v as usize + 1]]
+    }
+}
+
+/// A nested-dissection permutation of the `n × n` symmetric pattern
+/// (`perm[new] = old`). Only the pattern's structure is read; a
+/// nonsymmetric pattern is symmetrised.
+pub(super) fn nested_dissection(n: usize, row_ptr: &[usize], col_idx: &[usize]) -> Vec<usize> {
+    if n == 0 {
+        return Vec::new();
+    }
+    // Symmetrised adjacency lists over DOFs, self loops excluded.
+    let mut lists: Vec<Vec<u32>> = vec![Vec::new(); n];
+    for row in 0..n {
+        for &col in &col_idx[row_ptr[row]..row_ptr[row + 1]] {
+            if col != row {
+                lists[row].push(col as u32);
+                lists[col].push(row as u32);
+            }
+        }
+    }
+    for list in &mut lists {
+        list.sort_unstable();
+        list.dedup();
+    }
+
+    // Supervariables: identical closed neighbourhoods.
+    let mut group_of = vec![0u32; n];
+    let mut members: Vec<Vec<u32>> = Vec::new();
+    {
+        let mut index: HashMap<Vec<u32>, u32> = HashMap::new();
+        for v in 0..n {
+            let mut closed = lists[v].clone();
+            let at = closed.partition_point(|&x| x < v as u32);
+            closed.insert(at, v as u32);
+            let next = members.len() as u32;
+            let g = *index.entry(closed).or_insert(next);
+            if g == next {
+                members.push(Vec::new());
+            }
+            members[g as usize].push(v as u32);
+            group_of[v] = g;
+        }
+    }
+    let ng = members.len();
+    let mut ptr = Vec::with_capacity(ng + 1);
+    let mut adj = Vec::new();
+    ptr.push(0);
+    for group in &members {
+        let rep = group[0] as usize;
+        let mut nbrs: Vec<u32> = lists[rep]
+            .iter()
+            .map(|&u| group_of[u as usize])
+            .filter(|&g| g != group_of[rep])
+            .collect();
+        nbrs.sort_unstable();
+        nbrs.dedup();
+        adj.extend_from_slice(&nbrs);
+        ptr.push(adj.len());
+    }
+    drop(lists);
+    let graph = Graph {
+        ptr,
+        adj,
+        weight: members.iter().map(Vec::len).collect(),
+    };
+
+    let mut state = Dissector {
+        graph: &graph,
+        stamp: vec![0u32; ng],
+        level: vec![u32::MAX; ng],
+        next_stamp: 0,
+        out: Vec::with_capacity(ng),
+    };
+    let all: Vec<u32> = (0..ng as u32).collect();
+    state.dissect(all);
+
+    let mut perm = Vec::with_capacity(n);
+    for &g in &state.out {
+        for &v in &members[g as usize] {
+            perm.push(v as usize);
+        }
+    }
+    debug_assert_eq!(perm.len(), n);
+    perm
+}
+
+struct Dissector<'g> {
+    graph: &'g Graph,
+    /// `stamp[v] == id` marks membership of the subgraph being processed.
+    stamp: Vec<u32>,
+    level: Vec<u32>,
+    next_stamp: u32,
+    out: Vec<u32>,
+}
+
+impl Dissector<'_> {
+    fn mark(&mut self, verts: &[u32]) -> u32 {
+        self.next_stamp += 1;
+        let id = self.next_stamp;
+        for &v in verts {
+            self.stamp[v as usize] = id;
+        }
+        id
+    }
+
+    fn weight(&self, verts: &[u32]) -> usize {
+        verts.iter().map(|&v| self.graph.weight[v as usize]).sum()
+    }
+
+    /// BFS inside the stamped subgraph from `start`; fills `level` for
+    /// the reached vertices and returns them grouped by level.
+    fn bfs(&mut self, id: u32, start: u32) -> Vec<Vec<u32>> {
+        let mut levels: Vec<Vec<u32>> = vec![vec![start]];
+        self.level[start as usize] = 0;
+        let mut visited = vec![start];
+        loop {
+            let depth = levels.len() as u32;
+            let mut next = Vec::new();
+            for &v in levels.last().unwrap() {
+                for &u in self.graph.neighbours(v) {
+                    if self.stamp[u as usize] == id && self.level[u as usize] == u32::MAX {
+                        self.level[u as usize] = depth;
+                        next.push(u);
+                        visited.push(u);
+                    }
+                }
+            }
+            if next.is_empty() {
+                break;
+            }
+            next.sort_unstable();
+            levels.push(next);
+        }
+        for &v in &visited {
+            self.level[v as usize] = u32::MAX;
+        }
+        levels
+    }
+
+    fn dissect(&mut self, verts: Vec<u32>) {
+        if self.weight(&verts) <= LEAF_WEIGHT || verts.len() < 4 {
+            self.out.extend_from_slice(&verts);
+            return;
+        }
+        let id = self.mark(&verts);
+        // Connected components first.
+        let first = self.bfs(id, verts[0]);
+        let reached: usize = first.iter().map(Vec::len).sum();
+        if reached < verts.len() {
+            let mut component: Vec<u32> = first.into_iter().flatten().collect();
+            component.sort_unstable();
+            let comp_id = self.mark(&component);
+            let rest: Vec<u32> = verts
+                .iter()
+                .copied()
+                .filter(|&v| self.stamp[v as usize] != comp_id)
+                .collect();
+            self.dissect(component);
+            self.dissect(rest);
+            return;
+        }
+
+        // Pseudo-peripheral start: min-degree vertex of the last level,
+        // repeated while the eccentricity grows.
+        let degree = |g: &Graph, v: u32| g.ptr[v as usize + 1] - g.ptr[v as usize];
+        let mut start = *verts
+            .iter()
+            .min_by_key(|&&v| (degree(self.graph, v), v))
+            .unwrap();
+        let mut levels = self.bfs(id, start);
+        for _ in 0..8 {
+            let far = *levels
+                .last()
+                .unwrap()
+                .iter()
+                .min_by_key(|&&v| (degree(self.graph, v), v))
+                .unwrap();
+            let trial = self.bfs(id, far);
+            if trial.len() > levels.len() {
+                start = far;
+                levels = trial;
+            } else {
+                break;
+            }
+        }
+        let _ = start;
+        if levels.len() < 3 {
+            self.out.extend_from_slice(&verts);
+            return;
+        }
+
+        for (depth, level) in levels.iter().enumerate() {
+            for &v in level {
+                self.level[v as usize] = depth as u32;
+            }
+        }
+        let total = self.weight(&verts);
+        let level_weight: Vec<usize> = levels.iter().map(|l| self.weight(l)).collect();
+        // Candidate separators around the weighted median level.
+        let mut best: Option<(usize, usize)> = None; // (separator weight, level)
+        let mut before = 0usize;
+        for l in 0..levels.len() {
+            let after = total - before - level_weight[l];
+            if l >= 1 && l + 1 < levels.len() {
+                let sep: usize = levels[l]
+                    .iter()
+                    .filter(|&&v| {
+                        self.graph.neighbours(v).iter().any(|&u| {
+                            self.stamp[u as usize] == id && self.level[u as usize] == l as u32 + 1
+                        })
+                    })
+                    .map(|&v| self.graph.weight[v as usize])
+                    .sum();
+                let side_a = before + level_weight[l] - sep;
+                let balanced = 10 * side_a.min(after) >= 3 * total;
+                if balanced && best.is_none_or(|(w, _)| sep < w) {
+                    best = Some((sep, l));
+                }
+            }
+            before += level_weight[l];
+        }
+        let cut = match best {
+            Some((_, l)) => l,
+            None => {
+                // No balanced level: take the weighted median.
+                let mut acc = 0usize;
+                let mut pick = levels.len() / 2;
+                for (l, w) in level_weight.iter().enumerate() {
+                    acc += w;
+                    if 2 * acc >= total {
+                        pick = l;
+                        break;
+                    }
+                }
+                pick.clamp(1, levels.len() - 2)
+            }
+        };
+        let mut part_a = Vec::new();
+        let mut part_b = Vec::new();
+        let mut separator = Vec::new();
+        for &v in &verts {
+            let l = self.level[v as usize] as usize;
+            if l < cut {
+                part_a.push(v);
+            } else if l > cut {
+                part_b.push(v);
+            } else if self
+                .graph
+                .neighbours(v)
+                .iter()
+                .any(|&u| self.stamp[u as usize] == id && self.level[u as usize] == cut as u32 + 1)
+            {
+                separator.push(v);
+            } else {
+                part_a.push(v);
+            }
+        }
+        for &v in &verts {
+            self.level[v as usize] = u32::MAX;
+        }
+        self.dissect(part_a);
+        self.dissect(part_b);
+        self.out.extend_from_slice(&separator);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn grid_pattern(nx: usize, ny: usize) -> (usize, Vec<usize>, Vec<usize>) {
+        let n = nx * ny;
+        let mut row_ptr = vec![0];
+        let mut col_idx = Vec::new();
+        for j in 0..ny {
+            for i in 0..nx {
+                let mut cols = Vec::new();
+                for (di, dj) in [(-1i64, 0i64), (1, 0), (0, -1), (0, 1), (0, 0)] {
+                    let (a, b) = (i as i64 + di, j as i64 + dj);
+                    if a >= 0 && b >= 0 && (a as usize) < nx && (b as usize) < ny {
+                        cols.push(b as usize * nx + a as usize);
+                    }
+                }
+                cols.sort_unstable();
+                col_idx.extend(cols);
+                row_ptr.push(col_idx.len());
+            }
+        }
+        (n, row_ptr, col_idx)
+    }
+
+    #[test]
+    fn nested_dissection_is_a_permutation() {
+        let (n, rp, ci) = grid_pattern(40, 7);
+        let perm = nested_dissection(n, &rp, &ci);
+        let mut seen = vec![false; n];
+        for &p in &perm {
+            assert!(!seen[p]);
+            seen[p] = true;
+        }
+        assert_eq!(perm.len(), n);
+        // Deterministic.
+        assert_eq!(perm, nested_dissection(n, &rp, &ci));
+    }
+}

crates/specialized/rtx-fea/src/solvers/sparse_ldlt/symbolic.rs

@@ -0,0 +1,375 @@
+//! Symbolic analysis: ordering, elimination tree, postorder, column
+//! structures, fundamental supernodes, the supernodal assembly tree, and
+//! the maps the numeric phase scatters through. Computed once per
+//! sparsity pattern and reused for every numeric factorisation.
+
+use super::ordering::nested_dissection;
+
+/// Everything the numeric factorisation needs that depends only on the
+/// pattern. Column indices below are in the *factor* numbering
+/// (`perm[new] = old`).
+#[derive(Debug)]
+pub(super) struct Symbolic {
+    pub n: usize,
+    pub perm: Vec<usize>,
+    /// Supernode `s` owns columns `start[s]..start[s + 1]`.
+    pub start: Vec<usize>,
+    /// Front row indices of supernode `s` (its own columns first, then the
+    /// sorted off-diagonal structure): `rows[rows_ptr[s]..rows_ptr[s+1]]`.
+    pub rows_ptr: Vec<usize>,
+    pub rows: Vec<usize>,
+    /// Children in the assembly tree, increasing (postorder) order.
+    pub children_ptr: Vec<usize>,
+    pub children: Vec<usize>,
+    /// Roots of the assembly forest.
+    pub roots: Vec<usize>,
+    /// For a non-root supernode: position of each of its update rows
+    /// (`rows[k..m]`) in its parent's front — `relind[relind_ptr[s]..]`.
+    pub relind_ptr: Vec<usize>,
+    pub relind: Vec<u32>,
+    /// Per supernode, the lower-triangle entries of A it assembles:
+    /// `(CSR index, column-major offset in the front)`.
+    pub amap_ptr: Vec<usize>,
+    pub amap: Vec<(u32, u32)>,
+    /// The pattern this analysis was built for (for reuse checks).
+    pub pattern_row_ptr: Vec<usize>,
+    pub pattern_col_idx: Vec<usize>,
+    /// Estimated multiply–adds of one numeric factorisation, per subtree
+    /// (used to decide where to fork threads).
+    pub subtree_work: Vec<f64>,
+    /// Nonzeros in L (including the diagonal) and total work.
+    pub nnz_l: usize,
+    pub work: f64,
+}
+
+impl Symbolic {
+    pub fn num_supernodes(&self) -> usize {
+        self.start.len() - 1
+    }
+
+    pub fn front_rows(&self, s: usize) -> &[usize] {
+        &self.rows[self.rows_ptr[s]..self.rows_ptr[s + 1]]
+    }
+
+    pub fn node_children(&self, s: usize) -> &[usize] {
+        &self.children[self.children_ptr[s]..self.children_ptr[s + 1]]
+    }
+
+    pub fn max_front(&self) -> usize {
+        (0..self.num_supernodes())
+            .map(|s| self.rows_ptr[s + 1] - self.rows_ptr[s])
+            .max()
+            .unwrap_or(0)
+    }
+
+    /// Full analysis of an `n × n` CSR pattern (symmetric structure
+    /// expected; only the lower triangle is factorised).
+    pub fn analyze(n: usize, row_ptr: &[usize], col_idx: &[usize]) -> Self {
+        let nd = nested_dissection(n, row_ptr, col_idx);
+        // Elimination tree under the ND order, then postorder it.
+        let iperm_nd = inverse(&nd);
+        let rowpat = lower_rows(n, row_ptr, col_idx, &iperm_nd);
+        let parent_nd = etree(n, &rowpat);
+        let post = postorder(&parent_nd);
+        let perm: Vec<usize> = post.iter().map(|&k| nd[k]).collect();
+        drop(rowpat);
+
+        let iperm = inverse(&perm);
+        let rowpat = lower_rows(n, row_ptr, col_idx, &iperm);
+        let parent = etree(n, &rowpat);
+        // Column patterns of the strictly lower triangle of PAPᵀ.
+        let mut colpat: Vec<Vec<usize>> = vec![Vec::new(); n];
+        for (i, cols) in rowpat.iter().enumerate() {
+            for &j in cols {
+                colpat[j].push(i);
+            }
+        }
+        drop(rowpat);
+        let mut nchild = vec![0usize; n];
+        let mut child_lists: Vec<Vec<usize>> = vec![Vec::new(); n];
+        for j in 0..n {
+            if parent[j] != NONE {
+                nchild[parent[j]] += 1;
+                child_lists[parent[j]].push(j);
+            }
+        }
+
+        // Column structures by merging children (postorder: children of j
+        // precede j). A structure is kept only while it may still be
+        // needed: until its parent is processed, and for supernode heads.
+        let mut structs: Vec<Option<Vec<usize>>> = vec![None; n];
+        let mut counts = vec![0usize; n];
+        let mut start = vec![0usize];
+        let mut heads_struct: Vec<Vec<usize>> = Vec::new();
+        let mut mark = vec![NONE; n];
+        for j in 0..n {
+            let mut s: Vec<usize> = Vec::new();
+            mark[j] = j;
+            for &i in &colpat[j] {
+                if mark[i] != j {
+                    mark[i] = j;
+                    s.push(i);
+                }
+            }
+            for &c in &child_lists[j] {
+                if let Some(cs) = structs[c].take() {
+                    for i in cs {
+                        if i != j && mark[i] != j {
+                            mark[i] = j;
+                            s.push(i);
+                        }
+                    }
+                }
+            }
+            s.sort_unstable();
+            counts[j] = s.len();
+            // Does j extend the previous supernode?
+            let joins =
+                j > 0 && parent[j - 1] == j && nchild[j] == 1 && counts[j - 1] == counts[j] + 1;
+            if !joins {
+                if j > 0 {
+                    start.push(j);
+                }
+                heads_struct.push(s.clone());
+            }
+            structs[j] = Some(s);
+        }
+        start.push(n);
+        drop(structs);
+        drop(colpat);
+        let ns = start.len() - 1;
+
+        let mut col_to_sn = vec![0usize; n];
+        for s in 0..ns {
+            for c in start[s]..start[s + 1] {
+                col_to_sn[c] = s;
+            }
+        }
+        let mut rows_ptr = vec![0usize];
+        let mut rows = Vec::new();
+        for s in 0..ns {
+            rows.push(start[s]);
+            rows.extend_from_slice(&heads_struct[s]);
+            rows_ptr.push(rows.len());
+        }
+        drop(heads_struct);
+
+        let mut sn_parent = vec![NONE; ns];
+        for s in 0..ns {
+            let last = start[s + 1] - 1;
+            if parent[last] != NONE {
+                sn_parent[s] = col_to_sn[parent[last]];
+            }
+        }
+        let mut children_lists: Vec<Vec<usize>> = vec![Vec::new(); ns];
+        let mut roots = Vec::new();
+        for s in 0..ns {
+            if sn_parent[s] == NONE {
+                roots.push(s);
+            } else {
+                children_lists[sn_parent[s]].push(s);
+            }
+        }
+        let mut children_ptr = vec![0usize];
+        let mut children = Vec::new();
+        for list in &children_lists {
+            children.extend_from_slice(list);
+            children_ptr.push(children.len());
+        }
+
+        // Relative indices of each update row in the parent front.
+        let mut relind_ptr = vec![0usize];
+        let mut relind = Vec::new();
+        for s in 0..ns {
+            let k = start[s + 1] - start[s];
+            let own = &rows[rows_ptr[s]..rows_ptr[s + 1]];
+            if sn_parent[s] != NONE {
+                let p = sn_parent[s];
+                let prow = &rows[rows_ptr[p]..rows_ptr[p + 1]];
+                let mut q = 0usize;
+                for &r in &own[k..] {
+                    while prow[q] != r {
+                        q += 1;
+                    }
+                    relind.push(q as u32);
+                }
+            }
+            relind_ptr.push(relind.len());
+        }
+
+        let mut sym = Self {
+            n,
+            perm,
+            start,
+            rows_ptr,
+            rows,
+            children_ptr,
+            children,
+            roots,
+            relind_ptr,
+            relind,
+            amap_ptr: Vec::new(),
+            amap: Vec::new(),
+            pattern_row_ptr: Vec::new(),
+            pattern_col_idx: Vec::new(),
+            subtree_work: Vec::new(),
+            nnz_l: 0,
+            work: 0.0,
+        };
+        sym.build_amap(row_ptr, col_idx, &iperm, &col_to_sn);
+        sym.estimate_work();
+        sym
+    }
+
+    fn build_amap(
+        &mut self,
+        row_ptr: &[usize],
+        col_idx: &[usize],
+        iperm: &[usize],
+        col_to_sn: &[usize],
+    ) {
+        let ns = self.num_supernodes();
+        let mut per_sn: Vec<Vec<(u32, u32)>> = vec![Vec::new(); ns];
+        for old_row in 0..self.n {
+            let i = iperm[old_row];
+            for idx in row_ptr[old_row]..row_ptr[old_row + 1] {
+                let j = iperm[col_idx[idx]];
+                if i < j {
+                    continue;
+                }
+                let s = col_to_sn[j];
+                let front = self.front_rows(s);
+                let m = front.len();
+                let a = front
+                    .binary_search(&i)
+                    .expect("A entry outside the symbolic structure");
+                let b = j - self.start[s];
+                per_sn[s].push((idx as u32, (a + b * m) as u32));
+            }
+        }
+        self.amap_ptr = vec![0];
+        for list in per_sn {
+            self.amap.extend(list);
+            self.amap_ptr.push(self.amap.len());
+        }
+        self.pattern_row_ptr = row_ptr.to_vec();
+        self.pattern_col_idx = col_idx.to_vec();
+    }
+
+    fn estimate_work(&mut self) {
+        let ns = self.num_supernodes();
+        let mut own = vec![0.0f64; ns];
+        let mut nnz = 0usize;
+        let mut total = 0.0;
+        for (s, slot) in own.iter_mut().enumerate() {
+            let k = (self.start[s + 1] - self.start[s]) as f64;
+            let m = (self.rows_ptr[s + 1] - self.rows_ptr[s]) as f64;
+            let ku = k as usize;
+            let mu = m as usize;
+            nnz += ku * mu - ku * (ku - 1) / 2;
+            // Partial LDLᵀ of k pivots in an m-front (plus the extend-add).
+            let w = k * (m * m - m * k + k * k / 3.0) + (m - k) * (m - k);
+            *slot = w;
+            total += w;
+        }
+        // Subtree sums in postorder (children precede parents).
+        let mut sub = own;
+        for s in 0..ns {
+            for &c in self.node_children(s).to_vec().iter() {
+                sub[s] += sub[c];
+            }
+        }
+        self.subtree_work = sub;
+        self.nnz_l = nnz;
+        self.work = total;
+    }
+
+    /// Whether `(row_ptr, col_idx)` is exactly the analysed pattern.
+    pub fn same_pattern(&self, row_ptr: &[usize], col_idx: &[usize]) -> bool {
+        self.pattern_row_ptr == row_ptr && self.pattern_col_idx == col_idx
+    }
+}
+
+const NONE: usize = usize::MAX;
+
+fn inverse(perm: &[usize]) -> Vec<usize> {
+    let mut inv = vec![0usize; perm.len()];
+    for (new, &old) in perm.iter().enumerate() {
+        inv[old] = new;
+    }
+    inv
+}
+
+/// Row `i`'s strictly-lower columns `{j < i}` of the permuted,
+/// symmetrised pattern.
+fn lower_rows(n: usize, row_ptr: &[usize], col_idx: &[usize], iperm: &[usize]) -> Vec<Vec<usize>> {
+    let mut rows: Vec<Vec<usize>> = vec![Vec::new(); n];
+    for old_row in 0..n {
+        let a = iperm[old_row];
+        for &old_col in &col_idx[row_ptr[old_row]..row_ptr[old_row + 1]] {
+            let b = iperm[old_col];
+            if a > b {
+                rows[a].push(b);
+            } else if b > a {
+                rows[b].push(a);
+            }
+        }
+    }
+    for r in &mut rows {
+        r.sort_unstable();
+        r.dedup();
+    }
+    rows
+}
+
+/// Liu's elimination tree with path compression.
+fn etree(n: usize, rowpat: &[Vec<usize>]) -> Vec<usize> {
+    let mut parent = vec![NONE; n];
+    let mut ancestor = vec![NONE; n];
+    for i in 0..n {
+        for &j in &rowpat[i] {
+            let mut k = j;
+            while ancestor[k] != NONE && ancestor[k] != i {
+                let next = ancestor[k];
+                ancestor[k] = i;
+                k = next;
+            }
+            if ancestor[k] == NONE {
+                ancestor[k] = i;
+                parent[k] = i;
+            }
+        }
+    }
+    parent
+}
+
+/// Postorder of a forest (`post[k]` = the k-th node visited), children
+/// in increasing order, roots in increasing order; iterative DFS.
+fn postorder(parent: &[usize]) -> Vec<usize> {
+    let n = parent.len();
+    let mut kids: Vec<Vec<usize>> = vec![Vec::new(); n];
+    let mut roots = Vec::new();
+    for j in 0..n {
+        if parent[j] == NONE {
+            roots.push(j);
+        } else {
+            kids[parent[j]].push(j);
+        }
+    }
+    let mut post = Vec::with_capacity(n);
+    let mut stack: Vec<(usize, usize)> = Vec::new();
+    for &r in &roots {
+        stack.push((r, 0));
+        while let Some(&mut (v, ref mut next)) = stack.last_mut() {
+            if *next < kids[v].len() {
+                let c = kids[v][*next];
+                *next += 1;
+                stack.push((c, 0));
+            } else {
+                post.push(v);
+                stack.pop();
+            }
+        }
+    }
+    post
+}

crates/specialized/rtx-fea/tests/sparse_tangent.rs

@@ -0,0 +1,592 @@
+//! R8-g: the sparse Newton-tangent path of the nonlinear Newmark stepper
+//! ([`TangentSolver::SparseLdlt`]) against the default banded LU.
+//!
+//! Suite (fast, run by default):
+//!
+//! 1. `sparse_matches_banded_csm3_2d` — the 2-D 35×2 Quad8 CSM3 march
+//!    (the FSI2 harness's flag), 60 steps: every step's displacement
+//!    within 1e-10 of the banded path's (relative to the peak).
+//! 2. `sparse_matches_banded_csm1_3d_plane_strain` — CSM1 static on the
+//!    3-D 35×2×1 plane-strain flag.
+//! 3. `sparse_matches_banded_csm3_3d_free` — a free-edge 3-D flag
+//!    (12×2×2), 40 CSM3 steps.
+//! 4. `modified_newton_reuses_the_factor` — `reuse = 4`: fewer
+//!    factorisations than solves, same march to the Newton tolerance.
+//! 5. `force_only_kernel_is_bit_identical` — the TL force-only kernel
+//!    (modified Newton's reused iterations) = the pair's force, bit for
+//!    bit, on Hex20 and Quad8.
+//! 6. `tangent_knob_parses` — the solver enum's defaults.
+//!
+//! Instruments (`#[ignore]`, env-driven, write under `R8G_OUT`):
+//!
+//! * `r8g_g1_equivalence` — CSM1 static and a CSM3 march per
+//!   configuration with both solvers side by side; per-step CSV.
+//! * `r8g_g2_cost` — wall time per Newton iteration and per step, per
+//!   solver mode and configuration, with the sparse path's breakdown.
+
+use std::io::Write as _;
+use std::time::Instant;
+
+use nalgebra::{DVector, Vector3};
+use rtx_fea::analysis::flag3d::{Flag3d, Flag3dSpec, LateralFaces};
+use rtx_fea::analysis::{
+    ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis, TangentSolver,
+};
+use rtx_fea::assembly::dof_mapping::DofComponent;
+use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use rtx_fea::elements::{ElementMatrixComputer, StandardFiniteElement};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+
+const E_MOD: f64 = 1.4e6;
+const NU: f64 = 0.4;
+const RHO_CSM: f64 = 1000.0;
+const G: f64 = 2.0;
+
+fn env_str(name: &str, default: &str) -> String {
+    std::env::var(name).unwrap_or_else(|_| default.to_string())
+}
+
+fn env_num(name: &str, default: f64) -> f64 {
+    std::env::var(name)
+        .map(|v| v.parse().expect(name))
+        .unwrap_or(default)
+}
+
+/// The 2-D flag, as the FSI2 harness builds it (copied from
+/// `flag3d_structure.rs`).
+fn quad8_flag(nx: usize, ny: usize) -> Mesh {
+    let (x0, x1, y0, y1) = (0.25, 0.6, 0.19, 0.21);
+    let mut mesh = Mesh::new(2).unwrap();
+    let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
+    let mut grid = vec![vec![None; ly]; lx];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            if i % 2 == 1 && j % 2 == 1 {
+                continue;
+            }
+            let x = x0 + (x1 - x0) * i as f64 / (2 * nx) as f64;
+            let y = y0 + (y1 - y0) * j as f64 / (2 * ny) as f64;
+            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
+        }
+    }
+    for i in 0..nx {
+        for j in 0..ny {
+            let (a, b) = (2 * i, 2 * j);
+            let nodes = vec![
+                grid[a][b].unwrap(),
+                grid[a + 2][b].unwrap(),
+                grid[a + 2][b + 2].unwrap(),
+                grid[a][b + 2].unwrap(),
+                grid[a + 1][b].unwrap(),
+                grid[a + 2][b + 1].unwrap(),
+                grid[a + 1][b + 2].unwrap(),
+                grid[a][b + 1].unwrap(),
+            ];
+            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    mesh
+}
+
+fn clamp_2d(mesh: &Mesh) -> BoundaryConditionSet {
+    let clamped: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| (node.position().x - 0.25).abs() < 1e-12)
+        .map(|(&id, _)| id)
+        .collect();
+    let mut set = BoundaryConditionSet::new();
+    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
+        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
+            nodes: clamped.clone(),
+            components: vec![component],
+            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
+            time_range: None,
+            ramping_factor: 1.0,
+            gradual_enforcement: false,
+        }));
+    }
+    set
+}
+
+/// Consistent gravity nodal forces `∫ N_a ρ g dV`.
+fn gravity_forces(mesh: &Mesh, rho: f64, g: f64) -> Vec<(NodeId, Vector3<f64>)> {
+    let dim = mesh.spatial_dimension;
+    let mut acc: std::collections::BTreeMap<NodeId, Vector3<f64>> = Default::default();
+    for element in mesh.elements.values() {
+        let coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+        let f = ElementMatrixComputer::compute_body_force_vector(
+            &fe,
+            &coords,
+            &|_| Vector3::new(0.0, -rho * g, 0.0),
+            None,
+        )
+        .unwrap();
+        for (a, id) in element.nodes.iter().enumerate() {
+            let e = acc.entry(*id).or_insert_with(Vector3::zeros);
+            for c in 0..dim {
+                e[c] += f[a * dim + c];
+            }
+        }
+    }
+    acc.into_iter().collect()
+}
+
+fn static_criteria() -> ConvergenceCriteria {
+    ConvergenceCriteria {
+        force_tolerance: 1e-12,
+        displacement_tolerance: 1e-14,
+        max_iterations: 60,
+        ..ConvergenceCriteria::default()
+    }
+}
+
+/// The CSM3 analysis (gravity from rest) on the 2-D 35×2 flag.
+fn csm3_2d(dt: f64, solver: TangentSolver) -> NonlinearDynamicAnalysis {
+    let mesh = quad8_flag(35, 2);
+    let mut analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        Flag3d::materials(E_MOD, NU, RHO_CSM),
+        clamp_2d(&mesh),
+        dt,
+        1,
+        Default::default(),
+    )
+    .with_total_lagrangian()
+    .with_tangent_solver(solver);
+    analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0));
+    analysis
+}
+
+/// The CSM3 analysis on a 3-D flag (R8-b's instrument settings).
+fn csm3_3d(
+    flag: &Flag3d,
+    lateral: LateralFaces,
+    dt: f64,
+    solver: TangentSolver,
+) -> NonlinearDynamicAnalysis {
+    let mut analysis = flag
+        .dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, dt, 1, 0.5)
+        .with_tangent_solver(solver);
+    analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0));
+    analysis
+}
+
+/// `max_i |a_i − b_i| / max_i |b_i|`.
+fn rel_max(a: &DVector<f64>, b: &DVector<f64>) -> f64 {
+    let scale = b.amax().max(1e-300);
+    (a - b).amax() / scale
+}
+
+/// March both analyses `steps` steps side by side; per step the
+/// relative displacement difference. Returns (max rel, per-step rows).
+fn march_pair(
+    banded: &NonlinearDynamicAnalysis,
+    sparse: &NonlinearDynamicAnalysis,
+    steps: usize,
+) -> (f64, Vec<(f64, f64, usize, usize)>) {
+    let mut sb = banded.stepper().unwrap();
+    let mut ss = sparse.stepper().unwrap();
+    let mut stb = sb.rest_state().unwrap();
+    let mut sts = ss.rest_state().unwrap();
+    let mut worst: f64 = rel_max(&sts.acceleration, &stb.acceleration);
+    let mut rows = Vec::new();
+    for _ in 0..steps {
+        let t0 = Instant::now();
+        let (nb, ib) = sb.step(&stb).unwrap();
+        let tb = t0.elapsed().as_secs_f64();
+        let t1 = Instant::now();
+        let (ns, is) = ss.step(&sts).unwrap();
+        let ts = t1.elapsed().as_secs_f64();
+        stb = nb;
+        sts = ns;
+        let r = rel_max(&sts.displacement, &stb.displacement);
+        worst = worst.max(r);
+        rows.push((r, tb / ts.max(1e-12), ib, is));
+    }
+    (worst, rows)
+}
+
+#[test]
+fn sparse_matches_banded_csm3_2d() {
+    let dt = 0.005;
+    let (worst, rows) = march_pair(
+        &csm3_2d(dt, TangentSolver::BandedLu),
+        &csm3_2d(dt, TangentSolver::SPARSE),
+        60,
+    );
+    let same_newton = rows.iter().all(|r| r.2 == r.3);
+    println!("CSM3 2-D 35x2, 60 steps: max rel |Δu| {worst:.3e}, same Newton counts {same_newton}");
+    assert!(worst < 1e-10, "sparse departs from banded: {worst:.3e}");
+    assert!(same_newton);
+}
+
+#[test]
+fn sparse_matches_banded_csm1_3d_plane_strain() {
+    let flag = Flag3d::build(Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1)).unwrap();
+    let forces = gravity_forces(&flag.mesh, RHO_CSM, G);
+    let solve = |solver: TangentSolver| {
+        let analysis = flag
+            .dynamic_analysis(E_MOD, NU, RHO_CSM, LateralFaces::PlaneStrain, 1e4, 1, 0.5)
+            .with_convergence_criteria(static_criteria())
+            .with_tangent_solver(solver);
+        let mut stepper = analysis.stepper().unwrap();
+        let n = stepper.rest_state().unwrap().displacement.len();
+        let mut u = DVector::zeros(n);
+        for s in 1..=5 {
+            let scale = s as f64 / 5.0;
+            let scaled: Vec<_> = forces.iter().map(|(id, f)| (*id, f * scale)).collect();
+            stepper.set_nodal_forces(&scaled);
+            let state = DynamicState {
+                displacement: u.clone(),
+                velocity: DVector::zeros(n),
+                acceleration: DVector::zeros(n),
+            };
+            u = stepper.step(&state).unwrap().0.displacement;
+        }
+        let d = stepper.node_dofs(flag.point_a());
+        (u.clone(), u[d[1]])
+    };
+    let (ub, ay_b) = solve(TangentSolver::BandedLu);
+    let (us, ay_s) = solve(TangentSolver::SPARSE);
+    let r = rel_max(&us, &ub);
+    println!(
+        "CSM1 3-D 35x2x1 ps: uy(A) banded {ay_b:.12e} sparse {ay_s:.12e}, max rel |Δu| {r:.3e}"
+    );
+    assert!(r < 1e-10, "CSM1 departs: {r:.3e}");
+    assert!((ay_b + 65.1406e-3).abs() < 1e-6, "CSM1 uy(A) moved: {ay_b}");
+}
+
+#[test]
+fn sparse_matches_banded_csm3_3d_free() {
+    let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 12, 2, 2)).unwrap();
+    let dt = 0.005;
+    let (worst, _) = march_pair(
+        &csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::BandedLu),
+        &csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::SPARSE),
+        40,
+    );
+    println!("CSM3 3-D 12x2x2 free, 40 steps: max rel |Δu| {worst:.3e}");
+    assert!(worst < 1e-10, "sparse departs from banded: {worst:.3e}");
+}
+
+#[test]
+fn modified_newton_reuses_the_factor() {
+    let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 12, 2, 2)).unwrap();
+    let dt = 0.005;
+    let full = csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::SPARSE);
+    let modified = csm3_3d(
+        &flag,
+        LateralFaces::Free,
+        dt,
+        TangentSolver::SparseLdlt { reuse: 4 },
+    );
+    let mut sf = full.stepper().unwrap();
+    let mut sm = modified.stepper().unwrap();
+    let mut stf = sf.rest_state().unwrap();
+    let mut stm = sm.rest_state().unwrap();
+    for _ in 0..40 {
+        stf = sf.step(&stf).unwrap().0;
+        stm = sm.step(&stm).unwrap().0;
+    }
+    let r = rel_max(&stm.displacement, &stf.displacement);
+    let stats_f = sf.tangent_stats().unwrap();
+    let stats_m = sm.tangent_stats().unwrap();
+    println!(
+        "modified Newton (reuse 4) vs full after 40 steps: rel {r:.3e}; full {stats_f:?}; modified {stats_m:?}"
+    );
+    assert_eq!(stats_f.factorizations, stats_f.solves);
+    assert!(stats_m.factorizations < stats_m.solves);
+    assert_eq!(stats_m.analyses, 1);
+    // Newton's own displacement tolerance is 1e-6 (relative).
+    assert!(r < 1e-4, "modified Newton drifts: {r:.3e}");
+}
+
+#[test]
+fn force_only_kernel_is_bit_identical() {
+    use rtx_fea::elements::total_lagrangian::{
+        internal_force, internal_force_and_tangent, saint_venant_kirchhoff,
+    };
+    let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 3, 1, 1)).unwrap();
+    let quad = quad8_flag(3, 1);
+    for mesh in [&flag.mesh, &quad] {
+        let dim = mesh.spatial_dimension;
+        let constitutive = saint_venant_kirchhoff(8.0e5, 5.0e5, dim);
+        for (e, element) in mesh.elements.values().enumerate() {
+            let coords: Vec<Vector3<f64>> = element
+                .nodes
+                .iter()
+                .map(|id| mesh.get_node(*id).unwrap().position())
+                .collect();
+            let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+            let n = coords.len() * dim;
+            let u = DVector::from_fn(n, |i, _| 1e-3 * ((i + 7 * e) as f64 * 0.73).sin());
+            let (f_pair, _) =
+                internal_force_and_tangent(&fe, &coords, &u, constitutive.as_ref(), None).unwrap();
+            let f_only = internal_force(&fe, &coords, &u, constitutive.as_ref(), None).unwrap();
+            assert!(f_pair.norm() > 0.0);
+            assert_eq!(f_pair, f_only, "force-only kernel differs");
+        }
+    }
+}
+
+#[test]
+fn tangent_knob_parses() {
+    // Only the parser (the knob is read when a stepper is built).
+    assert_eq!(TangentSolver::default(), TangentSolver::BandedLu);
+    assert_eq!(
+        TangentSolver::SPARSE,
+        TangentSolver::SparseLdlt { reuse: 1 }
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Instruments
+// ---------------------------------------------------------------------------
+
+/// `NXxNYxNZ:span:free|ps` entries (as in `flag3d_structure.rs`); `2d`
+/// = the 2-D 35×2 Quad8 flag.
+fn parse_configs(spec: &str) -> Vec<Option<(usize, usize, usize, f64, LateralFaces)>> {
+    spec.split(',')
+        .map(|entry| {
+            if entry.trim() == "2d" {
+                return None;
+            }
+            let mut parts = entry.trim().split(':');
+            let mesh = parts.next().unwrap();
+            let span: f64 = parts.next().unwrap().parse().unwrap();
+            let lateral = match parts.next().unwrap() {
+                "free" => LateralFaces::Free,
+                "ps" => LateralFaces::PlaneStrain,
+                other => panic!("lateral {other}"),
+            };
+            let n: Vec<usize> = mesh.split('x').map(|t| t.parse().unwrap()).collect();
+            Some((n[0], n[1], n[2], span, lateral))
+        })
+        .collect()
+}
+
+fn tag(cfg: &Option<(usize, usize, usize, f64, LateralFaces)>) -> String {
+    match cfg {
+        None => "2d_35x2".to_string(),
+        Some((nx, ny, nz, span, lateral)) => {
+            let l = if *lateral == LateralFaces::Free {
+                "free"
+            } else {
+                "ps"
+            };
+            format!("{nx}x{ny}x{nz}_s{span}_{l}")
+        }
+    }
+}
+
+fn build_flag(cfg: &(usize, usize, usize, f64, LateralFaces)) -> Flag3d {
+    let (nx, ny, nz, span, _) = *cfg;
+    Flag3d::build(Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz)).unwrap()
+}
+
+fn analysis_for(
+    cfg: &Option<(usize, usize, usize, f64, LateralFaces)>,
+    flag: Option<&Flag3d>,
+    dt: f64,
+    solver: TangentSolver,
+) -> NonlinearDynamicAnalysis {
+    match cfg {
+        None => csm3_2d(dt, solver),
+        Some(c) => csm3_3d(flag.unwrap(), c.4, dt, solver),
+    }
+}
+
+#[test]
+#[ignore = "instrument: G1 — sparse vs banded, CSM1 static and a CSM3 march"]
+fn r8g_g1_equivalence() {
+    let out = env_str("R8G_OUT", ".");
+    let configs = parse_configs(&env_str("R8G_CONFIGS", "35x2x1:0.05:ps,35x2x8:0.41:free"));
+    let steps = env_num("R8G_STEPS", 200.0) as usize;
+    let dt = env_num("R8G_DT", 0.005);
+    let do_csm1 = env_str("R8G_CSM1", "1") == "1";
+    let mut table = std::fs::OpenOptions::new()
+        .create(true)
+        .append(true)
+        .open(format!("{out}/g1_table.txt"))
+        .unwrap();
+    for cfg in &configs {
+        let name = tag(cfg);
+        let flag = cfg.as_ref().map(build_flag);
+        if do_csm1 {
+            if let (Some(c), Some(flag)) = (cfg, flag.as_ref()) {
+                let forces = gravity_forces(&flag.mesh, RHO_CSM, G);
+                let mut results = Vec::new();
+                for solver in [TangentSolver::BandedLu, TangentSolver::SPARSE] {
+                    let start = Instant::now();
+                    let analysis = flag
+                        .dynamic_analysis(E_MOD, NU, RHO_CSM, c.4, 1e4, 1, 0.5)
+                        .with_convergence_criteria(static_criteria())
+                        .with_tangent_solver(solver);
+                    let mut stepper = analysis.stepper().unwrap();
+                    let n = stepper.rest_state().unwrap().displacement.len();
+                    let mut u = DVector::zeros(n);
+                    let mut newton = 0;
+                    for s in 1..=5 {
+                        let scale = s as f64 / 5.0;
+                        let scaled: Vec<_> =
+                            forces.iter().map(|(id, f)| (*id, f * scale)).collect();
+                        stepper.set_nodal_forces(&scaled);
+                        let state = DynamicState {
+                            displacement: u.clone(),
+                            velocity: DVector::zeros(n),
+                            acceleration: DVector::zeros(n),
+                        };
+                        let (next, it) = stepper.step(&state).unwrap();
+                        newton += it;
+                        u = next.displacement;
+                    }
+                    let d = stepper.node_dofs(flag.point_a());
+                    results.push((
+                        u.clone(),
+                        u[d[0]],
+                        u[d[1]],
+                        newton,
+                        start.elapsed().as_secs_f64(),
+                    ));
+                }
+                let r = rel_max(&results[1].0, &results[0].0);
+                let line = format!(
+                    "G1 CSM1 {name}: A banded ux {:.9e} uy {:.9e} [{} Newton, {:.1} s] | sparse ux \
+                     {:.9e} uy {:.9e} [{} Newton, {:.1} s] | max rel |Δu| {r:.3e}",
+                    results[0].1,
+                    results[0].2,
+                    results[0].3,
+                    results[0].4,
+                    results[1].1,
+                    results[1].2,
+                    results[1].3,
+                    results[1].4
+                );
+                println!("{line}");
+                writeln!(table, "{line}").unwrap();
+            }
+        }
+        let banded = analysis_for(cfg, flag.as_ref(), dt, TangentSolver::BandedLu);
+        let sparse = analysis_for(cfg, flag.as_ref(), dt, TangentSolver::SPARSE);
+        let start = Instant::now();
+        let (worst, rows) = march_pair(&banded, &sparse, steps);
+        let path = format!("{out}/g1_csm3_{name}_dt{dt}.csv");
+        let mut csv = std::fs::File::create(&path).unwrap();
+        writeln!(csv, "step,rel_u,speedup,newton_banded,newton_sparse").unwrap();
+        for (k, (r, speed, ib, is)) in rows.iter().enumerate() {
+            writeln!(csv, "{},{r:.6e},{speed:.3},{ib},{is}", k + 1).unwrap();
+        }
+        let differ = rows.iter().filter(|r| r.2 != r.3).count();
+        let line = format!(
+            "G1 CSM3 {name} dt {dt}: {steps} steps, max rel |Δu| {worst:.3e}, Newton-count \
+             differences {differ}, {:.0} s → {path}",
+            start.elapsed().as_secs_f64()
+        );
+        println!("{line}");
+        writeln!(table, "{line}").unwrap();
+    }
+}
+
+#[test]
+#[ignore = "instrument: G2 — cost per Newton iteration and per step"]
+fn r8g_g2_cost() {
+    let out = env_str("R8G_OUT", ".");
+    let configs = parse_configs(&env_str(
+        "R8G_CONFIGS",
+        "2d,35x2x8:0.41:free,35x2x16:0.41:free",
+    ));
+    let modes: Vec<String> = env_str("R8G_MODES", "banded,sparse,sparse:3")
+        .split(',')
+        .map(str::to_string)
+        .collect();
+    let steps = env_num("R8G_STEPS", 40.0) as usize;
+    let warm = env_num("R8G_WARM", 20.0) as usize;
+    let dt = env_num("R8G_DT", 0.005);
+    let mut table = std::fs::OpenOptions::new()
+        .create(true)
+        .append(true)
+        .open(format!("{out}/g2_table.txt"))
+        .unwrap();
+    for cfg in &configs {
+        let name = tag(cfg);
+        let flag = cfg.as_ref().map(build_flag);
+        let mut reference: Option<DVector<f64>> = None;
+        for mode in &modes {
+            let solver = match mode.as_str() {
+                "banded" => TangentSolver::BandedLu,
+                "sparse" => TangentSolver::SPARSE,
+                other => TangentSolver::SparseLdlt {
+                    reuse: other.strip_prefix("sparse:").unwrap().parse().unwrap(),
+                },
+            };
+            let analysis = analysis_for(cfg, flag.as_ref(), dt, solver);
+            let t_build = Instant::now();
+            let mut stepper = analysis.stepper().unwrap();
+            let mut state = stepper.rest_state().unwrap();
+            let build = t_build.elapsed().as_secs_f64();
+            // Warm-up steps (off the clock) move the flag off rest.
+            for _ in 0..warm {
+                state = stepper.step(&state).unwrap().0;
+            }
+            let before = stepper.tangent_stats();
+            let mut times = Vec::with_capacity(steps);
+            let mut newton = 0usize;
+            for _ in 0..steps {
+                let t = Instant::now();
+                let (next, it) = stepper.step(&state).unwrap();
+                times.push(t.elapsed().as_secs_f64());
+                newton += it;
+                state = next;
+            }
+            let total: f64 = times.iter().sum();
+            let mut sorted = times.clone();
+            sorted.sort_by(f64::total_cmp);
+            let median = sorted[sorted.len() / 2];
+            let drift = reference
+                .as_ref()
+                .map_or(f64::NAN, |r| rel_max(&state.displacement, r));
+            if reference.is_none() {
+                reference = Some(state.displacement.clone());
+            }
+            let breakdown = match (before, stepper.tangent_stats()) {
+                (Some(b), Some(a)) => format!(
+                    " | assemblies {} ({:.4} s each), factorisations {} (scatter {:.4} + \
+                     numeric {:.4} s each), solves {} ({:.4} s each), fallbacks {}, nnz(L) {}",
+                    a.assemblies - b.assemblies,
+                    (a.assembly_seconds - b.assembly_seconds)
+                        / (a.assemblies - b.assemblies).max(1) as f64,
+                    a.factorizations - b.factorizations,
+                    (a.scatter_seconds - b.scatter_seconds)
+                        / (a.factorizations - b.factorizations).max(1) as f64,
+                    (a.factor_seconds - b.factor_seconds)
+                        / (a.factorizations - b.factorizations).max(1) as f64,
+                    a.solves - b.solves,
+                    (a.solve_seconds - b.solve_seconds) / (a.solves - b.solves).max(1) as f64,
+                    a.fallbacks,
+                    a.nnz_l
+                ),
+                _ => String::new(),
+            };
+            let line = format!(
+                "G2 {name} {mode} dt {dt}: dofs {} | build+rest {build:.3} s | {steps} steps after \
+                 {warm} warm: {:.4} s/step mean, {median:.4} median, {:.2} Newton/step, {:.4} \
+                 s/Newton | final rel to first mode {drift:.2e} | threads {}{breakdown}",
+                state.displacement.len(),
+                total / steps as f64,
+                newton as f64 / steps as f64,
+                total / newton.max(1) as f64,
+                rayon::current_num_threads()
+            );
+            println!("{line}");
+            writeln!(table, "{line}").unwrap();
+        }
+    }
+}

crates/specialized/rtx-fea/tests/sparse_tangent_bench.rs

@@ -0,0 +1,89 @@
+//! R8-g instrument: the sparse LDLᵀ against the banded LU on the dumped
+//! 3-D flag operators (R8-b's `k_*.coo` / `m_*.coo`, the TL tangent at
+//! u = 0 and the consistent mass on the free DOFs). The Newmark tangent
+//! `K + M/(β Δt²)` at CSM3's Δt = 0.005 is solved by both; the timing and
+//! the agreement are printed.
+//!
+//! `R8G_COO_DIR` (default the R8-b study dir), `R8G_TAGS` (comma list).
+
+use nalgebra::DVector;
+use rtx_fea::assembly::SparseMatrix;
+use rtx_fea::solvers::{BandedLu, LinearSolver, SolverOptions, SparseLdlt};
+use std::time::Instant;
+
+fn read_coo(path: &str) -> Vec<(usize, usize, f64)> {
+    let text = std::fs::read_to_string(path).unwrap_or_else(|e| panic!("{path}: {e}"));
+    text.lines()
+        .filter(|l| !l.trim().is_empty())
+        .map(|l| {
+            let mut it = l.split_whitespace();
+            let i: usize = it.next().unwrap().parse().unwrap();
+            let j: usize = it.next().unwrap().parse().unwrap();
+            let v: f64 = it.next().unwrap().parse().unwrap();
+            (i, j, v)
+        })
+        .collect()
+}
+
+#[test]
+#[ignore = "instrument: sparse LDLt vs banded LU on the dumped flag tangents"]
+fn ldlt_vs_banded_on_dumped_tangents() {
+    let dir = std::env::var("R8G_COO_DIR").unwrap_or_else(|_| {
+        "/home/osobh/projects/omni-cortex-data/fsi_studies/p1_regress/r8b".to_string()
+    });
+    let tags = std::env::var("R8G_TAGS")
+        .unwrap_or_else(|_| "35x2x1_s0.05_ps,35x2x8_s0.41_free".to_string());
+    let skip_banded = std::env::var("R8G_SKIP_BANDED").is_ok();
+    let coef = 1.0 / (0.25 * 0.005 * 0.005);
+    for tag in tags.split(',') {
+        let k = read_coo(&format!("{dir}/k_{tag}.coo"));
+        let m = read_coo(&format!("{dir}/m_{tag}.coo"));
+        let n = k.iter().map(|t| t.0.max(t.1)).max().unwrap() + 1;
+        let mut triplets = k.clone();
+        triplets.extend(m.iter().map(|&(i, j, v)| (i, j, coef * v)));
+        let a = SparseMatrix::from_triplets(n, n, &triplets).unwrap();
+        let x_exact = DVector::from_fn(n, |i, _| ((i as f64) * 0.37).sin() + 0.2);
+        let b = a.multiply_vector(&x_exact).unwrap();
+        let opts = SolverOptions::default();
+
+        let mut ldlt = SparseLdlt::new();
+        let t0 = Instant::now();
+        let (row_ptr, col_idx) = a.structure();
+        ldlt.analyze(n, row_ptr, col_idx).unwrap();
+        let t_sym = t0.elapsed().as_secs_f64();
+        let mut t_num = f64::MAX;
+        let mut t_solve = f64::MAX;
+        let mut x = DVector::zeros(n);
+        for _ in 0..5 {
+            let t1 = Instant::now();
+            ldlt.factorize_values(a.values()).unwrap();
+            t_num = t_num.min(t1.elapsed().as_secs_f64());
+            let t2 = Instant::now();
+            x = ldlt.solve_factored(&b).unwrap();
+            t_solve = t_solve.min(t2.elapsed().as_secs_f64());
+        }
+        let stats = ldlt.stats().unwrap();
+        let rel_exact = (&x - &x_exact).norm() / x_exact.norm();
+        let residual = (a.multiply_vector(&x).unwrap() - &b).norm() / b.norm();
+        let (t_band, rel_band) = if skip_banded {
+            (f64::NAN, f64::NAN)
+        } else {
+            let t3 = Instant::now();
+            let (xb, _) = BandedLu::new().solve(&a, &b, &opts).unwrap();
+            let t = t3.elapsed().as_secs_f64();
+            (t, (&x - &xb).norm() / xb.norm())
+        };
+        println!(
+            "{tag}: n {n} nnz(A) {} | LDLt symbolic {t_sym:.3} s, numeric {t_num:.4} s, solve \
+             {t_solve:.4} s | nnz(L) {} supernodes {} max front {} work {:.3e} | banded LU \
+             {t_band:.3} s | rel vs exact {rel_exact:.2e}, vs banded {rel_band:.2e}, residual \
+             {residual:.2e} | threads {}",
+            a.nnz(),
+            stats.nnz_l,
+            stats.supernodes,
+            stats.max_front,
+            stats.work,
+            rayon::current_num_threads()
+        );
+    }
+}