conspire/math/sparse/solver/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    SparseError,
6    factor::{CscLdl, CscLu},
7    matrix::CscMatrix,
8};
9use crate::math::{Scalar, Vector};
10use std::{
11    cell::{Ref, RefCell},
12    rc::Rc,
13};
14
15/// A sparse direct solver for repeated solves on a fixed sparsity pattern.
16///
17/// Cloning shares the cached factorization, whose pivot order and fill pattern
18/// are reused across solves until a pivot degrades. Symmetric values use an
19/// LDLᵀ factorization, falling back to LU otherwise. Symmetry is a caller-supplied
20/// guarantee about the source values, not detected at runtime — a source whose
21/// symmetry can vary between solves (e.g. a tangent that is only symmetric near
22/// a particular configuration) must be declared asymmetric.
23#[derive(Clone)]
24pub struct SparseSolver {
25    matrix: RefCell<CscMatrix>,
26    ldl: Rc<RefCell<Option<CscLdl>>>,
27    lu: Rc<RefCell<Option<CscLu>>>,
28}
29
30impl SparseSolver {
31    pub fn from_pattern(num: usize, pattern: Vec<(usize, usize)>, symmetric: bool) -> Self {
32        let matrix = CscMatrix::from_pattern(num, num, pattern);
33        let ldl = if symmetric {
34            matrix.ldl_symbolic().ok()
35        } else {
36            None
37        };
38        let lu = if ldl.is_none() {
39            matrix.lu_symbolic().ok()
40        } else {
41            None
42        };
43        Self {
44            matrix: RefCell::new(matrix),
45            ldl: Rc::new(RefCell::new(ldl)),
46            lu: Rc::new(RefCell::new(lu)),
47        }
48    }
49    /// The nonzero (row, column) positions this structure was built from.
50    pub fn pattern(&self) -> Ref<'_, [(usize, usize)]> {
51        Ref::map(self.matrix.borrow(), |matrix| matrix.pattern())
52    }
53    /// Solve a system of linear equations with values from a source,
54    /// refactoring the cached factorization when possible.
55    pub fn solve(
56        &self,
57        source: impl FnMut(usize, usize) -> Scalar,
58        b: &Vector,
59    ) -> Result<Vector, SparseError> {
60        let mut matrix = self.matrix.borrow_mut();
61        matrix.fill(source);
62        let mut ldl = self.ldl.borrow_mut();
63        if ldl.is_some() {
64            if let Some(cached) = ldl.as_mut()
65                && cached.refactor(&matrix).is_ok()
66            {
67                return Ok(cached.solve(b));
68            }
69            *ldl = None;
70        }
71        drop(ldl);
72        let mut lu = self.lu.borrow_mut();
73        if let Some(cached) = lu.as_mut()
74            && cached.refactor(&matrix).is_ok()
75        {
76            return Ok(cached.solve(b));
77        }
78        let fresh = match matrix.lu_symbolic() {
79            Ok(mut symbolic) => match symbolic.refactor(&matrix) {
80                Ok(()) => symbolic,
81                Err(_) => matrix.lu_amd()?,
82            },
83            Err(_) => matrix.lu_amd()?,
84        };
85        let solution = fresh.solve(b);
86        *lu = Some(fresh);
87        Ok(solution)
88    }
89}
conspire/math/sparse/solver/mod.rs

conspire/math/sparse/solver/
mod.rs
