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
134
135
136
137

use crate::error::KError;

/// How to fix "too-small" pivots during ILU numerics.
#[derive(Clone, Debug)]
pub struct PivotPolicy {
    pub mode: PivotMode,
    pub tau: f64,
    pub scale: PivotScale,
    pub sign: PivotSignPolicy,
    pub deterministic: bool,
}

impl Default for PivotPolicy {
    fn default() -> Self {
        Self {
            mode: PivotMode::AdditiveFloor,
            tau: 1e-12,
            scale: PivotScale::RowGershgorin,
            sign: PivotSignPolicy::Positive,
            deterministic: true,
        }
    }
}

/// The action to take when |u_ii| is below the floor.
#[derive(Clone, Copy, Debug)]
pub enum PivotMode {
    /// Error out (legacy Strict).
    Strict,
    /// Add an additive shift to reach the floor (recommended).
    AdditiveFloor,
    /// Clamp toward the floor by minimum additive change (alias of AdditiveFloor).
    ThresholdFloor,
    /// (Optional later) ILUTP pivoting path will hook here.
    PivotingAllowed,
}

/// Scaling for the floor: floor_i = tau * scale_i
#[derive(Clone, Copy, Debug)]
pub enum PivotScale {
    /// s_i = max_j |A_jj| over the input (precomputed once).
    MaxDiagA,
    /// s_i = |A_ii| (local diag magnitude).
    LocalDiagA,
    /// s_i = ||A_{i,:}||_inf  (row inf-norm, precompute once).
    RowInfA,
    /// s_i = |A_ii| + sum_{j!=i} |A_ij|  (Gershgorin-based).
    RowGershgorin,
    /// s_i = running max of |U_kk| for k <= i (fill-growth aware).
    RunningMaxU,
}

/// Sign handling
#[derive(Clone, Copy, Debug)]
pub enum PivotSignPolicy {
    /// Keep current sign(u_ii); if zero, treat as +.
    Preserve,
    /// Force positive (good default for SPD or when using CG).
    Positive,
}

/// Statistics gathered during pivot handling.
#[derive(Clone, Debug)]
pub struct PivotStats {
    pub num_floors: usize,
    pub max_abs_shift: f64,
    pub sum_abs_shift: f64,
    pub num_strict_fail: usize,
    pub last_floor_value: f64,
}

impl Default for PivotStats {
    fn default() -> Self {
        Self {
            num_floors: 0,
            max_abs_shift: 0.0,
            sum_abs_shift: 0.0,
            num_strict_fail: 0,
            last_floor_value: 0.0,
        }
    }
}

fn record_pivot_shift(stats: &mut PivotStats, old: f64, new: f64) {
    let sh: f64 = (new - old).abs();
    stats.num_floors += 1;
    if sh > stats.max_abs_shift {
        stats.max_abs_shift = sh;
    }
    stats.sum_abs_shift += sh;
    stats.last_floor_value = new.abs();
}

/// Apply minimal additive shift to stabilize pivot.
pub fn stabilize_pivot_in_place(
    u_ii: &mut f64,
    s_i: f64,
    tau: f64,
    sign_policy: PivotSignPolicy,
    mode: PivotMode,
    stats: &mut PivotStats,
    row: usize,
) -> Result<(), KError> {
    let floor = tau * s_i;
    let abs: f64 = (*u_ii).abs();

    match mode {
        PivotMode::Strict => {
            if abs < floor {
                stats.num_strict_fail += 1;
                stats.last_floor_value = floor;
                return Err(KError::ZeroPivot(row));
            }
        }
        PivotMode::AdditiveFloor | PivotMode::ThresholdFloor => {
            if abs >= floor {
                return Ok(());
            }

            let magnitude = floor;
            let new_value = match sign_policy {
                PivotSignPolicy::Preserve => {
                    let sign: f64 = if *u_ii >= 0.0 { 1.0 } else { -1.0 };
                    sign * magnitude
                }
                PivotSignPolicy::Positive => magnitude,
            };
            let old = *u_ii;
            *u_ii = new_value;
            record_pivot_shift(stats, old, *u_ii);
        }
        PivotMode::PivotingAllowed => {
            unimplemented!("Pivoting path not yet implemented");
        }
    }
    Ok(())
}