carapace_core/analysis/
solver.rs

1use std::sync::Mutex;
2
3use faer::prelude::*;
4use faer::sparse::linalg::solvers::{Lu, SymbolicLu};
5use faer::sparse::{SymbolicSparseColMat, SymbolicSparseColMatRef};
6use nalgebra::DVector;
7
8use crate::model::SparseMatrix;
9
10use super::AnalysisError;
11
12/// The one linear solver Carapace ships — a concrete struct, not a
13/// trait or enum, since there is exactly one implementation to commit to.
14///
15/// `faer`'s
16/// sparse LU (with built-in COLAMD/AMD fill-reducing ordering), not a dense
17/// `nalgebra` solve. A dense O(n^3) LU doesn't scale to real models. It builds and
18/// solves correctly under wasm32.
19pub struct SparseSolver {
20    // Keep factor/solve callable through &self and preserve Send + Sync.
21    // The lock is released before numeric factorization and solving.
22    symbolic_cache: Mutex<Option<CachedSymbolic>>,
23    #[cfg(test)]
24    symbolic_analyses: std::sync::atomic::AtomicUsize,
25}
26
27struct CachedSymbolic {
28    pattern: SymbolicSparseColMat<usize>,
29    factorization: SymbolicLu<usize>,
30}
31
32impl CachedSymbolic {
33    fn matches(&self, pattern: SymbolicSparseColMatRef<'_, usize>) -> bool {
34        self.pattern.nrows() == pattern.nrows()
35            && self.pattern.ncols() == pattern.ncols()
36            && self.pattern.col_ptr() == pattern.col_ptr()
37            && self.pattern.col_nnz() == pattern.col_nnz()
38            && self.pattern.row_idx() == pattern.row_idx()
39    }
40}
41
42impl SparseSolver {
43    pub fn new() -> Self {
44        Self {
45            symbolic_cache: Mutex::new(None),
46            #[cfg(test)]
47            symbolic_analyses: std::sync::atomic::AtomicUsize::new(0),
48        }
49    }
50
51    /// Factor `k`, without solving — the half of a solve that
52    /// `TangentStrategy::ReuseAtStepStart`/`Initial` (`docs/algorithms.md`
53    /// §3.1) let a caller skip repeating across iterations. `faer`'s `Lu`
54    /// is a fully owned value (no lifetime tied back to `k`), so the result
55    /// outlives the matrix it was formed from — see `SparseFactorization`.
56    /// Symbolic analysis is reused while the full CSC pattern matches;
57    /// values (including zeros) never participate in cache invalidation.
58    /// Numeric factorization is always recomputed here, independently of
59    /// callers' numeric tangent-reuse caches.
60    pub fn factor(&self, k: &SparseMatrix) -> Result<SparseFactorization, AnalysisError> {
61        let symbolic = {
62            let mut cache = self
63                .symbolic_cache
64                .lock()
65                .expect("symbolic LU cache poisoned");
66            let pattern = k.symbolic();
67            if !cache.as_ref().is_some_and(|cached| cached.matches(pattern)) {
68                let factorization =
69                    SymbolicLu::try_new(pattern).map_err(|_| AnalysisError::SingularSystem)?;
70                let pattern = pattern
71                    .to_owned()
72                    .map_err(|_| AnalysisError::SingularSystem)?;
73                *cache = Some(CachedSymbolic {
74                    pattern,
75                    factorization,
76                });
77                #[cfg(test)]
78                self.symbolic_analyses
79                    .fetch_add(1, std::sync::atomic::Ordering::Relaxed);
80            }
81            cache
82                .as_ref()
83                .expect("symbolic cache populated above")
84                .factorization
85                .clone()
86        };
87        Ok(SparseFactorization(
88            Lu::try_new_with_symbolic(symbolic, k.as_ref())
89                .map_err(|_| AnalysisError::SingularSystem)?,
90        ))
91    }
92
93    /// Factor-then-solve in one call — for `Algorithm::Linear` and anywhere
94    /// else a single ad-hoc solve is enough. Equivalent to `self.factor(k)?
95    /// .solve(r)`, kept as its own method so a one-off caller doesn't have
96    /// to name `SparseFactorization` just to discard it immediately.
97    pub fn solve(&self, k: &SparseMatrix, r: &DVector<f64>) -> Result<DVector<f64>, AnalysisError> {
98        Ok(self.factor(k)?.solve(r))
99    }
100}
101
102impl Default for SparseSolver {
103    fn default() -> Self {
104        Self::new()
105    }
106}
107
108/// An LU factorization held independently of the matrix it was formed
109/// from — cached across Newton iterations by `TangentStrategy::
110/// ReuseAtStepStart`/`Initial` so re-solving against the same tangent
111/// doesn't re-factor it (`docs/algorithms.md` §3.1's real prerequisite for
112/// those strategies to actually save anything).
113pub struct SparseFactorization(Lu<usize, f64>);
114
115impl SparseFactorization {
116    pub fn solve(&self, r: &DVector<f64>) -> DVector<f64> {
117        let n = r.len();
118        let rhs = faer::col::Col::from_fn(n, |i| r[i]);
119        let x = self.0.solve(&rhs);
120        DVector::from_fn(n, |i, _| x[i])
121    }
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127    use faer::sparse::Triplet;
128    use std::sync::atomic::Ordering;
129
130    fn matrix(n: usize, entries: &[(usize, usize, f64)]) -> SparseMatrix {
131        let triplets: Vec<_> = entries
132            .iter()
133            .map(|&(i, j, v)| Triplet::new(i, j, v))
134            .collect();
135        SparseMatrix::try_new_from_triplets(n, n, &triplets).unwrap()
136    }
137
138    fn assert_solution(factor: &SparseFactorization, rhs: &[f64], expected: &[f64]) {
139        let solution = factor.solve(&DVector::from_column_slice(rhs));
140        assert!((solution - DVector::from_column_slice(expected)).norm() < 1e-12);
141    }
142
143    #[test]
144    fn changing_values_and_pivots_reuses_symbolic_but_not_numeric_factors() {
145        let solver = SparseSolver::new();
146        let first = matrix(2, &[(0, 0, 4.0), (0, 1, 2.0), (1, 0, 3.0), (1, 1, 5.0)]);
147        let factor = solver.factor(&first).unwrap();
148        assert_solution(&factor, &[8.0, 13.0], &[1.0, 2.0]);
149
150        // A structural entry becomes zero and numeric LU must change pivots.
151        let second = matrix(2, &[(0, 0, 0.0), (0, 1, 2.0), (1, 0, 3.0), (1, 1, 4.0)]);
152        let updated = solver.factor(&second).unwrap();
153        assert_solution(&updated, &[4.0, 11.0], &[1.0, 2.0]);
154        assert_eq!(solver.symbolic_analyses.load(Ordering::Relaxed), 1);
155        // Previously returned numeric factors remain independently owned.
156        assert_solution(&factor, &[8.0, 13.0], &[1.0, 2.0]);
157    }
158
159    #[test]
160    fn changed_pattern_with_same_size_and_nnz_and_changed_dimensions_rebuild_cache() {
161        let solver = SparseSolver::new();
162        let upper = matrix(2, &[(0, 0, 4.0), (0, 1, 1.0), (1, 1, 3.0)]);
163        assert_solution(&solver.factor(&upper).unwrap(), &[6.0, 6.0], &[1.0, 2.0]);
164        let lower = matrix(2, &[(0, 0, 4.0), (1, 0, 1.0), (1, 1, 3.0)]);
165        assert_solution(&solver.factor(&lower).unwrap(), &[4.0, 7.0], &[1.0, 2.0]);
166        assert_eq!(solver.symbolic_analyses.load(Ordering::Relaxed), 2);
167
168        let larger = matrix(3, &[(0, 0, 2.0), (1, 1, 3.0), (2, 2, 4.0)]);
169        assert_solution(
170            &solver.factor(&larger).unwrap(),
171            &[2.0, 6.0, 12.0],
172            &[1.0, 2.0, 3.0],
173        );
174        assert_eq!(solver.symbolic_analyses.load(Ordering::Relaxed), 3);
175        assert_solution(&solver.factor(&lower).unwrap(), &[4.0, 7.0], &[1.0, 2.0]);
176        assert_eq!(solver.symbolic_analyses.load(Ordering::Relaxed), 4);
177    }
178
179    #[test]
180    fn singular_pattern_does_not_poison_cache_or_invalidate_owned_factors() {
181        let solver = SparseSolver::new();
182        let nonsingular = matrix(2, &[(0, 0, 2.0), (1, 1, 3.0)]);
183        let factor = solver.factor(&nonsingular).unwrap();
184        let singular = matrix(2, &[(1, 1, 3.0)]);
185        assert!(matches!(
186            solver.factor(&singular),
187            Err(AnalysisError::SingularSystem)
188        ));
189        assert_solution(&factor, &[2.0, 6.0], &[1.0, 2.0]);
190        assert_solution(
191            &solver.factor(&nonsingular).unwrap(),
192            &[2.0, 6.0],
193            &[1.0, 2.0],
194        );
195    }
196}
carapace_core/analysis/solver.rs

carapace_core/analysis/
solver.rs
