diffsol/linear_solver/faer/
lu.rs

1use crate::FaerContext;
2use crate::{error::LinearSolverError, linear_solver_error};
3
4use crate::{
5    error::DiffsolError, linear_solver::LinearSolver, FaerMat, FaerVec, Matrix,
6    NonLinearOpJacobian, Scalar,
7};
8
9use faer::{linalg::solvers::FullPivLu, linalg::solvers::Solve};
10/// A [LinearSolver] that uses the LU decomposition in the [`faer`](https://github.com/sarah-ek/faer-rs) library to solve the linear system.
11pub struct LU<T>
12where
13    T: Scalar,
14{
15    lu: Option<FullPivLu<T>>,
16    matrix: Option<FaerMat<T>>,
17}
18
19impl<T> Default for LU<T>
20where
21    T: Scalar,
22{
23    fn default() -> Self {
24        Self {
25            lu: None,
26            matrix: None,
27        }
28    }
29}
30
31impl<T: Scalar> LinearSolver<FaerMat<T>> for LU<T> {
32    fn set_linearisation<C: NonLinearOpJacobian<T = T, V = FaerVec<T>, M = FaerMat<T>>>(
33        &mut self,
34        op: &C,
35        x: &FaerVec<T>,
36        t: T,
37    ) {
38        let matrix = self.matrix.as_mut().expect("Matrix not set");
39        op.jacobian_inplace(x, t, matrix);
40        self.lu = Some(matrix.data.full_piv_lu());
41    }
42
43    fn solve_in_place(&self, x: &mut FaerVec<T>) -> Result<(), DiffsolError> {
44        if self.lu.is_none() {
45            return Err(linear_solver_error!(LuNotInitialized))?;
46        }
47        let lu = self.lu.as_ref().unwrap();
48        lu.solve_in_place(x.data.as_mut());
49        Ok(())
50    }
51
52    fn set_problem<
53        C: NonLinearOpJacobian<T = T, V = FaerVec<T>, M = FaerMat<T>, C = FaerContext>,
54    >(
55        &mut self,
56        op: &C,
57    ) {
58        let ncols = op.nstates();
59        let nrows = op.nout();
60        let matrix =
61            C::M::new_from_sparsity(nrows, ncols, op.jacobian_sparsity(), op.context().clone());
62        self.matrix = Some(matrix);
63    }
64}
diffsol/linear_solver/faer/lu.rs

diffsol/linear_solver/faer/
lu.rs
