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

#[allow(unused_imports)]
use crate::algebra::blas::{dot_conj, nrm2};
use crate::algebra::prelude::*;

/// Build a complex Givens rotation that zeroes the second entry of the vector [a; b].
/// Returns (c, s) where `c` is real and non-negative, and `s` is the complex sine.
#[inline]
pub fn build_complex_givens(a: S, b: S) -> (R, S) {
    let absa = a.abs();
    let absb = b.abs();
    if absb <= R::default() {
        return (1.0, S::zero());
    }
    if absa <= R::default() {
        return (0.0, b.conj() / S::from_real(absb));
    }

    let r = absa.hypot(absb);
    let c = absa / r;
    let s = (a * b.conj()) / S::from_real(r * absa);
    (c, s)
}

/// Apply a complex Givens rotation to the in-place vector [h_ik; h_ip1k].
#[inline]
pub fn apply_complex_givens(h_ik: &mut S, h_ip1k: &mut S, c: R, s: S) {
    let c_s = S::from_real(c);
    let t = c_s * *h_ik + s * *h_ip1k;
    let t2 = -s.conj() * *h_ik + c_s * *h_ip1k;
    *h_ik = t;
    *h_ip1k = t2;
}

/// Apply the previously accumulated Givens rotations to a Hessenberg column.
#[inline]
pub fn apply_prev_givens_to_col(hcol: &mut [S], k: usize, cs: &[R], sn: &[S]) {
    for i in 0..k {
        let (head, tail) = hcol.split_at_mut(i + 1);
        let hij = &mut head[i];
        let hij1 = &mut tail[0];
        apply_complex_givens(hij, hij1, cs[i], sn[i]);
    }
}

/// Apply a freshly built Givens rotation to the Hessenberg column and update g.
#[inline]
pub fn apply_new_givens_and_update_g(
    hcol: &mut [S],
    k: usize,
    cs: &mut [R],
    sn: &mut [S],
    g: &mut [S],
) {
    let (c, s) = build_complex_givens(hcol[k], hcol[k + 1]);
    cs[k] = c;
    sn[k] = s;
    let (head, tail) = hcol.split_at_mut(k + 1);
    let hik = &mut head[k];
    let hip1k = &mut tail[0];
    apply_complex_givens(hik, hip1k, c, s);

    let c_s = S::from_real(c);
    let gk = g[k];
    let gk1 = g[k + 1];
    g[k] = c_s * gk + s * gk1;
    g[k + 1] = -s.conj() * gk + c_s * gk1;
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn complex_givens_zeroes_second_entry() {
        let a: S = S::from_real(1.2);
        #[cfg(feature = "complex")]
        let b: S = num_complex::Complex64::new(-0.7, 0.5);
        #[cfg(not(feature = "complex"))]
        let b: S = S::from_real(-0.7);

        let (c, s) = build_complex_givens(a, b);
        let mut h0 = a;
        let mut h1 = b;
        apply_complex_givens(&mut h0, &mut h1, c, s);

        let scale = (a.abs() + b.abs()).max(1.0);
        assert!(h1.abs() <= 1e-12 * scale);
        assert!(c >= 0.0);
    }

    #[test]
    fn givens_handles_zero_head_entry() {
        let a = S::zero();
        #[cfg(feature = "complex")]
        let b: S = num_complex::Complex64::new(0.3, -0.4);
        #[cfg(not(feature = "complex"))]
        let b: S = S::from_real(-0.4);

        let (c, s) = build_complex_givens(a, b);
        let mut h0 = a;
        let mut h1 = b;
        apply_complex_givens(&mut h0, &mut h1, c, s);

        assert!(c.abs() <= 1e-14);
        assert!((h0.abs() - b.abs()).abs() <= 1e-12 * b.abs().max(1.0));
        assert!(h1.abs() <= 1e-12 * b.abs().max(1.0));
    }
}