tiny_solver/linear/
sparse_cholesky.rs

1use std::fmt::Debug;
2use std::ops::Mul;
3
4use faer::linalg::solvers::Solve;
5use faer::sparse::linalg::solvers;
6
7use super::sparse::SparseLinearSolver;
8
9// #[pyclass]
10#[derive(Debug, Clone)]
11pub struct SparseCholeskySolver {
12    symbolic_pattern: Option<solvers::SymbolicLlt<usize>>,
13}
14
15impl SparseCholeskySolver {
16    pub fn new() -> Self {
17        SparseCholeskySolver {
18            symbolic_pattern: None,
19        }
20    }
21}
22impl Default for SparseCholeskySolver {
23    fn default() -> Self {
24        Self::new()
25    }
26}
27impl SparseLinearSolver for SparseCholeskySolver {
28    fn solve(
29        &mut self,
30        residuals: &faer::Mat<f64>,
31        jacobians: &faer::sparse::SparseColMat<usize, f64>,
32    ) -> Option<faer::Mat<f64>> {
33        let jtj = jacobians
34            .as_ref()
35            .transpose()
36            .to_col_major()
37            .unwrap()
38            .mul(jacobians.as_ref());
39        let jtr = jacobians.as_ref().transpose().mul(-residuals);
40
41        self.solve_jtj(&jtr, &jtj)
42    }
43
44    fn solve_jtj(
45        &mut self,
46        jtr: &faer::Mat<f64>,
47        jtj: &faer::sparse::SparseColMat<usize, f64>,
48    ) -> Option<faer::Mat<f64>> {
49        // initialize the pattern
50        if self.symbolic_pattern.is_none() {
51            self.symbolic_pattern =
52                Some(solvers::SymbolicLlt::try_new(jtj.symbolic(), faer::Side::Lower).unwrap());
53        }
54
55        let sym = self.symbolic_pattern.as_ref().unwrap();
56        if let Ok(cholesky) =
57            solvers::Llt::try_new_with_symbolic(sym.clone(), jtj.as_ref(), faer::Side::Lower)
58        {
59            let dx = cholesky.solve(jtr);
60
61            Some(dx)
62        } else {
63            None
64        }
65    }
66}
tiny_solver/linear/sparse_cholesky.rs

tiny_solver/linear/
sparse_cholesky.rs
