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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

use super::*;

fn simple_identity(n: usize) -> CsrMatrix<f64> {
    let mut row_ptr = Vec::with_capacity(n + 1);
    let mut col_idx = Vec::with_capacity(n);
    let mut vals = Vec::with_capacity(n);
    row_ptr.push(0);
    for i in 0..n {
        col_idx.push(i);
        vals.push(1.0);
        row_ptr.push(col_idx.len());
    }
    CsrMatrix::from_csr(n, n, row_ptr, col_idx, vals)
}

#[test]
fn powdered_columns_retry() {
    let a = simple_identity(2);
    let diag_inv = vec![1.0, 1.0];
    let tp = TentativeP {
        agg_of: vec![0, 1],
        n_coarse: 2,
        num_functions: 1,
        nns: None,
        comp_of: None,
    };
    let mut p_csr = Pcsr {
        m: 2,
        n: 2,
        row_ptr: vec![0, 1, 2],
        col_idx: vec![0, 1],
        vals: vec![1e-16, 1e-16],
    };
    let ctx = LevelPostContext {
        r: 1,
        agg_of: &tp.agg_of,
        nns: None,
        a: Some(&a),
        d_inv: Some(&diag_inv),
    };
    let p0 = CsrMatrix::from_csr(
        p_csr.m,
        p_csr.n,
        p_csr.row_ptr.clone(),
        p_csr.col_idx.clone(),
        p_csr.vals.clone(),
    );
    let mut cfg = AMGConfig::default();
    cfg.rank_min_col_norm = 1e-12;
    let diag0 = check_p_rank_fast(&p0, &cfg).unwrap();
    assert!(diag0.suspect, "expected suspect rank before fix");

    match try_fix_rank(0, &a, &diag_inv, &tp, &ctx, &mut p_csr, &mut cfg).unwrap() {
        RankFixOutcome::Fixed => {
            let p_fixed = CsrMatrix::from_csr(
                p_csr.m,
                p_csr.n,
                p_csr.row_ptr.clone(),
                p_csr.col_idx.clone(),
                p_csr.vals.clone(),
            );
            let diag1 = check_p_rank_fast(&p_fixed, &cfg).unwrap();
            assert!(diag1.min_col_norm > cfg.rank_min_col_norm);
            assert!(!diag1.suspect);
        }
        RankFixOutcome::Unfixed => panic!("expected retry to succeed"),
    }
}

#[test]
fn sketch_condition_detects_colinear_columns() {
    let row_ptr = vec![0, 2, 4];
    let col_idx = vec![0, 1, 0, 1];
    let vals = vec![1.0, 1.0, 2.0, 2.0];
    let p = CsrMatrix::from_csr(2, 2, row_ptr.clone(), col_idx.clone(), vals);
    let mut cfg = AMGConfig::default();
    cfg.rank_cond_threshold = 10.0;
    let diag = check_p_rank_fast(&p, &cfg).unwrap();
    assert!(diag.suspect, "colinear columns should trigger suspect flag");
    assert!(!diag.cond_estimate.is_finite() || diag.cond_estimate > cfg.rank_cond_threshold);
}

#[test]
fn galerkin_identity_matches() {
    let a = CsrMatrix::from_csr(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![4.0, -1.0, -1.0, 4.0],
    );
    let p = simple_identity(2);
    let r = simple_identity(2);
    let a_c = a.clone();
    let (ok, worst) = galerkin_sample_check(&a, &p, &r, &a_c, 3, 1e-12, 0xBEEF).unwrap();
    assert!(
        ok,
        "Galekin identity should hold for identical hierarchy: worst={worst:e}"
    );
}

#[test]
fn galerkin_identity_detects_mismatch() {
    let a = CsrMatrix::from_csr(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![4.0, -1.0, -1.0, 4.0],
    );
    let p = simple_identity(2);
    let r = simple_identity(2);
    let a_c = CsrMatrix::from_csr(
        2,
        2,
        vec![0, 2, 4],
        vec![0, 1, 0, 1],
        vec![4.0, -2.0, -1.0, 4.0],
    );
    let (ok, worst) = galerkin_sample_check(&a, &p, &r, &a_c, 2, 1e-12, 0xC0FF).unwrap();
    assert!(!ok);
    assert!(worst > 0.0);
}

#[test]
fn rank_condition_deterministic() {
    let p = CsrMatrix::from_csr(3, 3, vec![0, 1, 2, 3], vec![0, 1, 2], vec![1.0, 0.5, 0.25]);
    let mut cfg = AMGConfig::default();
    cfg.rank_sketch_cols = 3;
    let (_, c1) = rank_condition_estimate(&p, cfg.rank_sketch_cols, 0xABC).unwrap();
    let (_, c2) = rank_condition_estimate(&p, cfg.rank_sketch_cols, 0xABC).unwrap();
    assert_eq!(c1, c2);
}