diffsol_la/linear_solver/faer/
lu.rs

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

diffsol_la/linear_solver/faer/
lu.rs
