ndarray_linalg/krylov/
arnoldi.rs

1//! Arnoldi iteration
2
3use super::*;
4use crate::{norm::Norm, operator::LinearOperator};
5use num_traits::One;
6use std::iter::*;
7
8/// Execute Arnoldi iteration as Rust iterator
9///
10/// - [Arnoldi iteration - Wikipedia](https://en.wikipedia.org/wiki/Arnoldi_iteration)
11///
12pub struct Arnoldi<A, S, F, Ortho>
13where
14    A: Scalar,
15    S: DataMut<Elem = A>,
16    F: LinearOperator<Elem = A>,
17    Ortho: Orthogonalizer<Elem = A>,
18{
19    a: F,
20    /// Next vector (normalized `|v|=1`)
21    v: ArrayBase<S, Ix1>,
22    /// Orthogonalizer
23    ortho: Ortho,
24    /// Coefficients to be composed into H-matrix
25    h: Vec<Array1<A>>,
26}
27
28impl<A, S, F, Ortho> Arnoldi<A, S, F, Ortho>
29where
30    A: Scalar + Lapack,
31    S: DataMut<Elem = A>,
32    F: LinearOperator<Elem = A>,
33    Ortho: Orthogonalizer<Elem = A>,
34{
35    /// Create an Arnoldi iterator from any linear operator `a`
36    pub fn new(a: F, mut v: ArrayBase<S, Ix1>, mut ortho: Ortho) -> Self {
37        assert_eq!(ortho.len(), 0);
38        assert!(ortho.tolerance() < One::one());
39        // normalize before append because |v| may be smaller than ortho.tolerance()
40        let norm = v.norm_l2();
41        azip!((v in &mut v)  *v = v.div_real(norm));
42        ortho.append(v.view());
43        Arnoldi {
44            a,
45            v,
46            ortho,
47            h: Vec::new(),
48        }
49    }
50
51    /// Dimension of Krylov subspace
52    pub fn dim(&self) -> usize {
53        self.ortho.len()
54    }
55
56    /// Iterate until convergent
57    pub fn complete(mut self) -> (Q<A>, H<A>) {
58        for _ in &mut self {} // execute iteration until convergent
59        let q = self.ortho.get_q();
60        let n = self.h.len();
61        let mut h = Array2::zeros((n, n).f());
62        for (i, hc) in self.h.iter().enumerate() {
63            let m = std::cmp::min(n, i + 2);
64            for j in 0..m {
65                h[(j, i)] = hc[j];
66            }
67        }
68        (q, h)
69    }
70}
71
72impl<A, S, F, Ortho> Iterator for Arnoldi<A, S, F, Ortho>
73where
74    A: Scalar + Lapack,
75    S: DataMut<Elem = A>,
76    F: LinearOperator<Elem = A>,
77    Ortho: Orthogonalizer<Elem = A>,
78{
79    type Item = Array1<A>;
80
81    fn next(&mut self) -> Option<Self::Item> {
82        self.a.apply_mut(&mut self.v);
83        let result = self.ortho.div_append(&mut self.v);
84        let norm = self.v.norm_l2();
85        azip!((v in &mut self.v) *v = v.div_real(norm));
86        match result {
87            AppendResult::Added(coef) => {
88                self.h.push(coef.clone());
89                Some(coef)
90            }
91            AppendResult::Dependent(coef) => {
92                self.h.push(coef);
93                None
94            }
95        }
96    }
97}
98
99/// Utility to execute Arnoldi iteration with Householder reflection
100pub fn arnoldi_householder<A, S>(
101    a: impl LinearOperator<Elem = A>,
102    v: ArrayBase<S, Ix1>,
103    tol: A::Real,
104) -> (Q<A>, H<A>)
105where
106    A: Scalar + Lapack,
107    S: DataMut<Elem = A>,
108{
109    let householder = Householder::new(v.len(), tol);
110    Arnoldi::new(a, v, householder).complete()
111}
112
113/// Utility to execute Arnoldi iteration with modified Gram-Schmit orthogonalizer
114pub fn arnoldi_mgs<A, S>(
115    a: impl LinearOperator<Elem = A>,
116    v: ArrayBase<S, Ix1>,
117    tol: A::Real,
118) -> (Q<A>, H<A>)
119where
120    A: Scalar + Lapack,
121    S: DataMut<Elem = A>,
122{
123    let mgs = MGS::new(v.len(), tol);
124    Arnoldi::new(a, v, mgs).complete()
125}
ndarray_linalg/krylov/arnoldi.rs

ndarray_linalg/krylov/
arnoldi.rs
