diffsol_la/linear_solver/
mod.rs

1use crate::{error::LaError, LinearOp, Matrix};
2
3#[cfg(feature = "nalgebra")]
4pub mod nalgebra;
5
6#[cfg(feature = "faer")]
7pub mod faer;
8
9#[cfg(feature = "suitesparse")]
10pub mod suitesparse;
11
12#[cfg(feature = "cuda")]
13pub mod cuda;
14
15pub use faer::lu::LU as FaerLU;
16pub use nalgebra::lu::LU as NalgebraLU;
17
18/// A solver for the linear problem `Ax = b`, where `A` is a [LinearOp].
19pub trait LinearSolver<M: Matrix>: Default {
20    /// Set the point at which the linear operator `A` is evaluated and factorise it.
21    ///
22    /// The operator is assumed to have the same sparsity as that given to
23    /// [Self::set_sparsity].
24    fn set_linearisation<C: LinearOp<V = M::V, T = M::T, M = M, C = M::C>>(&mut self, op: &C);
25
26    /// Set the sparsity of the problem to be solved, any previous problem is discarded.
27    ///
28    /// Any internal state of the solver is reset. This function will normally set
29    /// the sparsity pattern of the matrix to be solved.
30    fn set_sparsity<C: LinearOp<V = M::V, T = M::T, M = M, C = M::C>>(&mut self, op: &C);
31
32    /// Solve the problem `Ax = b` and return the solution `x`.
33    ///
34    /// Panics if [Self::set_linearisation] has not been called previously.
35    fn solve(&self, b: &M::V) -> Result<M::V, LaError> {
36        let mut b = b.clone();
37        self.solve_in_place(&mut b)?;
38        Ok(b)
39    }
40
41    fn solve_in_place(&self, b: &mut M::V) -> Result<(), LaError>;
42}
43
44#[cfg(test)]
45pub(crate) mod tests {
46    use crate::{IndexType, LinearOp, Matrix, Vector};
47
48    /// A simple diagonal [LinearOp] used for testing the linear solvers.
49    pub struct DiagonalOp<M: Matrix> {
50        matrix: M,
51    }
52
53    /// Create a 2x2 diagonal operator `A = diag(value, value)`.
54    pub fn diagonal_op<M: Matrix>(value: f64) -> DiagonalOp<M> {
55        use num_traits::FromPrimitive;
56        let v = M::T::from_f64(value).unwrap();
57        let diag = M::V::from_vec(vec![v, v], Default::default());
58        DiagonalOp {
59            matrix: M::from_diagonal(&diag),
60        }
61    }
62
63    impl<M: Matrix> LinearOp for DiagonalOp<M> {
64        type T = M::T;
65        type V = M::V;
66        type M = M;
67        type C = M::C;
68
69        fn nrows(&self) -> IndexType {
70            self.matrix.nrows()
71        }
72        fn ncols(&self) -> IndexType {
73            self.matrix.ncols()
74        }
75        fn context(&self) -> &Self::C {
76            self.matrix.context()
77        }
78        fn matrix_inplace(&self, y: &mut Self::M) {
79            y.copy_from(&self.matrix);
80        }
81        fn sparsity(&self) -> Option<<Self::M as Matrix>::Sparsity> {
82            self.matrix.sparsity().map(|s| {
83                use crate::matrix::sparsity::MatrixSparsityRef;
84                s.to_owned()
85            })
86        }
87    }
88}
diffsol_la/linear_solver/mod.rs

diffsol_la/linear_solver/
mod.rs
