kryst 3.2.1

Krylov subspace and preconditioned iterative solvers for dense and sparse linear systems, with shared and distributed memory parallelism.
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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

use crate::algebra::prelude::*;
use crate::error::KError;
use crate::matrix::op::LinOp;
use crate::preconditioner::{PcSide, Preconditioner};
use faer::Mat;

#[cfg(feature = "complex")]
use crate::algebra::bridge::BridgeScratch;
#[cfg(feature = "complex")]
use crate::ops::kpc::KPreconditioner;
#[cfg(feature = "complex")]
use crate::preconditioner::bridge::apply_pc_s;
#[cfg(feature = "mpi")]
use crate::parallel::Comm;

pub struct LuPc;

impl LuPc {
    pub fn new() -> Self {
        Self
    }
}

impl Preconditioner for LuPc {
    fn setup(&mut self, pmat: &dyn LinOp<S = S>) -> Result<(), KError> {
        if pmat.comm().size() > 1 {
            return Err(KError::Unsupported(
                "LU PC is local-only under MPI; use SuperLU_DIST (feature superlu_dist) or a distributed preconditioner."
                    .into(),
            ));
        }
        pmat.as_any()
            .downcast_ref::<Mat<f64>>()
            .ok_or_else(|| KError::InvalidInput("LU PC requires faer::Mat<f64>".into()))?;
        Ok(())
    }

    fn apply(&self, _side: PcSide, r: &[S], z: &mut [S]) -> Result<(), KError> {
        let _ = (r, z);
        Err(KError::Unsupported(
            "LuPc::apply is PREONLY-only; use SolverType::Preonly or call direct_solve",
        ))
    }

    fn direct_solve(
        &mut self,
        pmat: &dyn LinOp<S = S>,
        b: &[S],
        x: &mut [S],
    ) -> Result<(), KError> {
        let a = pmat
            .as_any()
            .downcast_ref::<Mat<f64>>()
            .ok_or_else(|| KError::InvalidInput("LU PC requires faer::Mat<f64>".into()))?;
        #[cfg(not(feature = "complex"))]
        {
            crate::solver::dense_lu::solve(a, b, x)
        }
        #[cfg(feature = "complex")]
        {
            let n = b.len();
            let mut b_re = vec![0.0; n];
            let mut b_im = vec![0.0; n];
            for (i, bi) in b.iter().enumerate() {
                b_re[i] = bi.real();
                b_im[i] = bi.imag();
            }
            let mut x_re = vec![0.0; x.len()];
            let mut x_im = vec![0.0; x.len()];
            crate::solver::dense_lu::solve(a, &b_re, &mut x_re)?;
            crate::solver::dense_lu::solve(a, &b_im, &mut x_im)?;
            for (xi, (&re, &im)) in x.iter_mut().zip(x_re.iter().zip(x_im.iter())) {
                *xi = S::from_parts(re, im);
            }
            Ok(())
        }
    }

    fn required_format(&self) -> crate::matrix::format::OpFormat {
        crate::matrix::format::OpFormat::Dense
    }
}

#[cfg(feature = "complex")]
impl KPreconditioner for LuPc {
    type Scalar = S;

    #[inline]
    fn dims(&self) -> (usize, usize) {
        <Self as Preconditioner>::dims(self)
    }

    #[inline]
    fn apply_s(
        &self,
        side: PcSide,
        x: &[S],
        y: &mut [S],
        scratch: &mut BridgeScratch,
    ) -> Result<(), KError> {
        apply_pc_s(self, side, x, y, scratch)
    }
}