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

#![cfg(all(feature = "backend-faer", not(feature = "complex")))]
#![cfg(all(feature = "complex", feature = "complex_ilu"))]

use faer::Mat;
use kryst::algebra::prelude::*;
use kryst::preconditioner::PcSide;
use kryst::preconditioner::ilu::{IluBuilder, IluType};
use kryst::preconditioner::legacy::Preconditioner;

fn make_diagonal_matrix(diag: &[S]) -> Mat<S> {
    let n = diag.len();
    Mat::from_fn(n, n, |i, j| if i == j { diag[i] } else { S::zero() })
}

fn mat_vec_mul(a: &Mat<S>, x: &[S]) -> Vec<S> {
    let n = a.nrows();
    let mut out = vec![S::zero(); n];
    for i in 0..n {
        let mut sum = S::zero();
        for j in 0..a.ncols() {
            sum = sum + a[(i, j)] * x[j];
        }
        out[i] = sum;
    }
    out
}

#[test]
fn ilu0_complex_diagonal_matches_inverse() {
    let diag = [
        S::from_parts(3.0, -0.5),
        S::from_parts(2.0, 1.25),
        S::from_parts(4.5, -2.0),
        S::from_parts(5.0, 0.75),
    ];
    let matrix = make_diagonal_matrix(&diag);

    let mut ilu = IluBuilder::new()
        .ilu_type(IluType::ILU0)
        .build::<S>()
        .expect("complex ILU0 build");
    ilu.setup(&matrix).expect("complex ILU0 setup");

    let rhs = [
        S::from_parts(1.0, -2.0),
        S::from_parts(-3.0, 0.5),
        S::from_parts(0.25, 1.5),
        S::from_parts(-0.75, -1.25),
    ];
    let mut z = vec![S::zero(); rhs.len()];
    ilu.apply(PcSide::Left, &rhs, &mut z)
        .expect("complex ILU0 apply");

    for (i, (&expected_diag, &rhs_i)) in diag.iter().zip(rhs.iter()).enumerate() {
        let expected = rhs_i / expected_diag;
        let diff = z[i] - expected;
        assert!(
            diff.abs() < 1e-12,
            "entry {} mismatch: got {:?}, expected {:?}",
            i,
            z[i],
            expected
        );
    }
}

#[test]
fn ilu0_diagonally_dominant_has_finite_solution() {
    let matrix = Mat::from_fn(3, 3, |i, j| {
        if i == j {
            S::from_parts(6.0 + i as f64, (-1.5 + i as f64) * 0.1)
        } else if (i as isize - j as isize).abs() == 1 {
            S::from_parts(0.25 * (1 + i + j) as f64, -0.3 * (i as f64 - j as f64))
        } else {
            S::zero()
        }
    });

    let mut ilu = IluBuilder::new()
        .ilu_type(IluType::ILU0)
        .build::<S>()
        .expect("complex ILU0 build");
    ilu.setup(&matrix)
        .expect("complex ILU0 setup for diagonally dominant matrix");

    assert_eq!(ilu.pivot_stats().num_floors, 0, "unexpected pivot floors");

    let rhs = [
        S::from_parts(1.5, -0.25),
        S::from_parts(-2.25, 0.5),
        S::from_parts(0.75, 1.0),
    ];
    let mut z = vec![S::zero(); rhs.len()];
    ilu.apply(PcSide::Left, &rhs, &mut z)
        .expect("complex ILU0 apply");

    for (idx, value) in z.iter().enumerate() {
        assert!(
            value.real().is_finite() && value.imag().is_finite(),
            "non-finite output at {}: {:?}",
            idx,
            value
        );
    }

    let residual = mat_vec_mul(&matrix, &z)
        .iter()
        .zip(rhs.iter())
        .map(|(&az, &rhs_i)| (az - rhs_i).abs())
        .fold(0.0, f64::max);
    assert!(residual < 1e-10, "residual too large: {residual}");
}