faer-precond 0.1.0

Numerical preconditioners (Jacobi, block-Jacobi, ILU(0), IC(0)) for iterative solvers, built on faer.
Documentation 
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

//! Triangular solve dispatch for IC(0).

use faer::sparse::linalg::triangular_solve;
use faer::{Conj, MatMut, Par};
use faer_traits::{ComplexField, Index};

use super::numeric::Ic0;

/// Apply the IC(0) preconditioner in place.
///
/// `conj = Conj::No`  produces `out = M^{-1} rhs`.
/// `conj = Conj::Yes` produces `out = conj(M)^{-1} rhs = M^{-T} rhs`
///                    (these coincide for Hermitian `M`).
pub(crate) fn solve_in_place<I, T>(ic: &Ic0<I, T>, conj: Conj, rhs: MatMut<'_, T>, par: Par)
where
    I: Index,
    T: ComplexField,
{
    let l = ic.l_view();
    let mut rhs = rhs;
    match conj {
        Conj::No => {
            // M^{-1} = L^{-H} L^{-1}
            triangular_solve::solve_lower_triangular_in_place(l, Conj::No, rhs.as_mut(), par);
            triangular_solve::solve_lower_triangular_transpose_in_place(
                l,
                Conj::Yes,
                rhs.as_mut(),
                par,
            );
        }
        Conj::Yes => {
            // conj(M)^{-1} = L^{-T} conj(L)^{-1}
            triangular_solve::solve_lower_triangular_in_place(l, Conj::Yes, rhs.as_mut(), par);
            triangular_solve::solve_lower_triangular_transpose_in_place(
                l,
                Conj::No,
                rhs.as_mut(),
                par,
            );
        }
    }
}