tiny_solver/linear/
sparse_cholesky.rs

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

tiny_solver/linear/
sparse_cholesky.rs
