diffsol_la/linear_solver/faer/
sparse_lu.rs

1use crate::{
2    error::LaError, linear_solver::LinearSolver, linear_solver_error, scalar::IndexType,
3    FaerContext, FaerScalar, FaerSparseMat, FaerVec, LinearOp, Matrix,
4};
5
6use faer::{
7    linalg::solvers::Solve,
8    reborrow::Reborrow,
9    sparse::linalg::{solvers::Lu, solvers::SymbolicLu},
10};
11
12/// A [LinearSolver] that uses the LU decomposition in the [`faer`](https://github.com/sarah-ek/faer-rs) library to solve the linear system.
13pub struct FaerSparseLU<T>
14where
15    T: FaerScalar,
16{
17    lu: Option<Lu<IndexType, T>>,
18    lu_symbolic: Option<SymbolicLu<IndexType>>,
19    matrix: Option<FaerSparseMat<T>>,
20}
21
22impl<T> Default for FaerSparseLU<T>
23where
24    T: FaerScalar,
25{
26    fn default() -> Self {
27        Self {
28            lu: None,
29            matrix: None,
30            lu_symbolic: None,
31        }
32    }
33}
34
35impl<T: FaerScalar> LinearSolver<FaerSparseMat<T>> for FaerSparseLU<T> {
36    fn set_linearisation<
37        C: LinearOp<T = T, V = FaerVec<T>, M = FaerSparseMat<T>, C = FaerContext>,
38    >(
39        &mut self,
40        op: &C,
41    ) {
42        let matrix = self.matrix.as_mut().expect("Matrix not set");
43        op.matrix_inplace(matrix);
44        self.lu = Some(
45            Lu::try_new_with_symbolic(self.lu_symbolic.as_ref().unwrap().clone(), matrix.data.rb())
46                .expect("Failed to factorise matrix"),
47        );
48    }
49
50    fn solve_in_place(&self, x: &mut FaerVec<T>) -> Result<(), LaError> {
51        let lu = self
52            .lu
53            .as_ref()
54            .ok_or_else(|| linear_solver_error!(LuNotInitialized))?;
55        lu.solve_in_place(&mut x.data);
56        Ok(())
57    }
58
59    fn set_sparsity<C: LinearOp<T = T, V = FaerVec<T>, M = FaerSparseMat<T>, C = FaerContext>>(
60        &mut self,
61        op: &C,
62    ) {
63        let ncols = op.ncols();
64        let nrows = op.nrows();
65        let matrix = C::M::new_from_sparsity(nrows, ncols, op.sparsity(), *op.context());
66        self.lu_symbolic = Some(
67            SymbolicLu::try_new(matrix.data.symbolic()).expect("Failed to create symbolic LU"),
68        );
69        self.matrix = Some(matrix);
70    }
71}
72
73#[cfg(test)]
74mod tests {
75    use super::*;
76    use crate::{linear_solver::tests::diagonal_op, Vector};
77
78    #[test]
79    fn test_sparse_lu() {
80        let mut s = FaerSparseLU::<f64>::default();
81        let op = diagonal_op::<FaerSparseMat<f64>>(2.0);
82        s.set_sparsity(&op);
83        s.set_linearisation(&op);
84        let b = FaerVec::from_vec(vec![2.0, 4.0], Default::default());
85        let x = s.solve(&b).unwrap();
86        x.assert_eq_st(
87            &FaerVec::from_vec(vec![1.0, 2.0], Default::default()),
88            1e-10,
89        );
90    }
91}
diffsol_la/linear_solver/faer/sparse_lu.rs

diffsol_la/linear_solver/faer/
sparse_lu.rs
