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

use crate::algebra::prelude::*;
use crate::context::ksp_context::Workspace;
use crate::parallel::{NoComm, UniverseComm};
use crate::preconditioner::PcSide;
use crate::solver::LinearSolver;
use crate::solver::block::{BlockKrylovOptions, bicgstab::BlockBicgstabSolver, gmres::BlockGmresSolver};
use crate::utils::convergence::ConvergedReason;

#[test]
fn block_gmres_converges_on_block_rhs() {
    let mut opts = BlockKrylovOptions::default();
    opts.block_size = 2;
    opts.restart_blocks = 4;
    opts.max_iters = 50;
    let mut solver = BlockGmresSolver::new(opts);
    let a = super::util::spd_poisson2d(4);
    let n = a.nrows();
    let mut b: Vec<R> = vec![R::default(); n * 2];
    for col in 0..2 {
        for row in 0..n {
            b[col * n + row] = R::from(1.0 + (col + row) as f64);
        }
    }
    let mut x: Vec<R> = vec![R::default(); n * 2];
    let mut ws = Workspace::default();
    let comm = UniverseComm::NoComm(NoComm);
    let op: &dyn crate::matrix::op::LinOp<S = f64> = &a;
    let stats = solver
        .solve(
            op,
            None,
            &b,
            &mut x,
            PcSide::Left,
            &comm,
            None,
            Some(&mut ws),
        )
        .unwrap();
    assert!(matches!(
        stats.reason,
        ConvergedReason::ConvergedAtol | ConvergedReason::ConvergedRtol
    ));
}

#[test]
fn block_bicgstab_converges_on_block_rhs() {
    let mut opts = BlockKrylovOptions::default();
    opts.block_size = 2;
    opts.max_iters = 100;
    let mut solver = BlockBicgstabSolver::new(opts);
    let a = super::util::spd_poisson2d(4);
    let n = a.nrows();
    let mut b: Vec<R> = vec![R::default(); n * 2];
    for col in 0..2 {
        for row in 0..n {
            b[col * n + row] = R::from(1.0 + (col + row) as f64);
        }
    }
    let mut x: Vec<R> = vec![R::default(); n * 2];
    let mut ws = Workspace::default();
    let comm = UniverseComm::NoComm(NoComm);
    let op: &dyn crate::matrix::op::LinOp<S = f64> = &a;
    let stats = solver
        .solve(
            op,
            None,
            &b,
            &mut x,
            PcSide::Left,
            &comm,
            None,
            Some(&mut ws),
        )
        .unwrap();
    assert!(matches!(
        stats.reason,
        ConvergedReason::ConvergedAtol | ConvergedReason::ConvergedRtol
    ));
}