tiny_solver/
linear.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
use std::ops::Mul;

use faer::prelude::SpSolver;
use faer::sparse::linalg::solvers;

// #[pyclass]
#[derive(Default, Clone)]
pub enum LinearSolver {
    #[default]
    SparseCholesky,
    SparseQR,
}

pub fn sparse_cholesky(
    residuals: &faer::Mat<f64>,
    jacobians: &faer::sparse::SparseColMat<usize, f64>,
    symbolic_pattern: &mut Option<solvers::SymbolicCholesky<usize>>,
) -> Option<faer::Mat<f64>> {
    let hessian = jacobians
        .as_ref()
        .transpose()
        .to_col_major()
        .unwrap()
        .mul(jacobians.as_ref());
    let b = jacobians.as_ref().transpose().mul(-residuals);
    let sym = if symbolic_pattern.is_some() {
        symbolic_pattern.as_ref().unwrap()
    } else {
        // initialize the pattern
        *symbolic_pattern = Some(
            solvers::SymbolicCholesky::try_new(hessian.symbolic(), faer::Side::Lower).unwrap(),
        );
        symbolic_pattern.as_ref().unwrap()
    };
    if let Ok(cholesky) =
        solvers::Cholesky::try_new_with_symbolic(sym.clone(), hessian.as_ref(), faer::Side::Lower)
    {
        let dx = cholesky.solve(b);

        Some(dx)
    } else {
        None
    }
}

pub fn sparse_qr(
    residuals: &faer::Mat<f64>,
    jacobians: &faer::sparse::SparseColMat<usize, f64>,
) -> faer::Mat<f64> {
    jacobians.sp_qr().unwrap().solve(-residuals)
}