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
104
105
106
107
108
109
110
111
112
113
114
115
116

#![cfg(not(feature = "complex"))]

use std::any::Any;
use std::sync::Arc;

use kryst::context::ksp_context::{KspContext, SolverType};
#[cfg(all(
    feature = "backend-faer",
    feature = "dense-direct",
    not(feature = "complex")
))]
use kryst::context::pc_context::PcType;
use kryst::matrix::op::{LinOp, StructureId, ValuesId};
use kryst::parallel::{NoComm, UniverseComm};
use kryst::utils::convergence::ConvergedReason;

#[derive(Clone)]
struct TinyDiagOp {
    diag: [f64; 2],
    comm: UniverseComm,
    sid: StructureId,
    vid: ValuesId,
}

impl TinyDiagOp {
    fn new(diag: [f64; 2], comm: UniverseComm) -> Self {
        Self {
            diag,
            comm,
            sid: StructureId(1),
            vid: ValuesId(1),
        }
    }
}

impl LinOp for TinyDiagOp {
    type S = f64;

    fn dims(&self) -> (usize, usize) {
        (self.diag.len(), self.diag.len())
    }

    fn matvec(&self, x: &[Self::S], y: &mut [Self::S]) {
        for (yi, (&d, &xi)) in y.iter_mut().zip(self.diag.iter().zip(x.iter())) {
            *yi = d * xi;
        }
    }

    fn as_any(&self) -> &dyn Any {
        self
    }

    fn structure_id(&self) -> StructureId {
        self.sid
    }

    fn values_id(&self) -> ValuesId {
        self.vid
    }

    fn comm(&self) -> UniverseComm {
        self.comm.clone()
    }
}

#[test]
fn context_smoke_cg() {
    let op = Arc::new(TinyDiagOp::new([2.0, 4.0], UniverseComm::NoComm(NoComm)));
    let b = vec![2.0, 8.0];
    let mut x = vec![0.0; 2];

    let mut ksp = KspContext::new();
    ksp.rtol = 1e-12;
    ksp.atol = 1e-12;
    ksp.set_type(SolverType::Cg).unwrap();
    ksp.set_operators(op, None);
    ksp.setup().unwrap();

    let stats = ksp.solve(&b, &mut x).unwrap();
    assert!(matches!(
        stats.reason,
        ConvergedReason::ConvergedRtol | ConvergedReason::ConvergedAtol
    ));
    assert!((x[0] - 1.0).abs() < 1e-10);
    assert!((x[1] - 2.0).abs() < 1e-10);
    assert!(stats.final_residual < 1e-10);
}

#[cfg(all(
    feature = "backend-faer",
    feature = "dense-direct",
    not(feature = "complex")
))]
#[test]
fn context_smoke_preonly_lu() {
    use faer::Mat;
    use kryst::matrix::op::DenseOp;

    let mat = Arc::new(Mat::from_fn(2, 2, |i, j| {
        if i == j { 3.0 + i as f64 } else { 0.0 }
    }));
    let op = Arc::new(DenseOp::<f64>::new(mat));
    let b = vec![3.0, 8.0];
    let mut x = vec![0.0; 2];

    let mut ksp = KspContext::new();
    ksp.set_preonly_with_pc(PcType::Lu, None).unwrap();
    ksp.set_operators(op, None);
    ksp.setup().unwrap();

    let stats = ksp.solve(&b, &mut x).unwrap();
    assert!(matches!(stats.reason, ConvergedReason::ConvergedAtol));
    assert!((x[0] - 1.0).abs() < 1e-12);
    assert!((x[1] - 2.0).abs() < 1e-12);
    assert!(stats.final_residual < 1e-12);
}