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
use super::*;
use crate::{norm::Norm, operator::LinearOperator};
use num_traits::One;
use std::iter::*;
pub struct Arnoldi<A, S, F, Ortho>
where
A: Scalar,
S: DataMut<Elem = A>,
F: LinearOperator<Elem = A>,
Ortho: Orthogonalizer<Elem = A>,
{
a: F,
v: ArrayBase<S, Ix1>,
ortho: Ortho,
h: Vec<Array1<A>>,
}
impl<A, S, F, Ortho> Arnoldi<A, S, F, Ortho>
where
A: Scalar + Lapack,
S: DataMut<Elem = A>,
F: LinearOperator<Elem = A>,
Ortho: Orthogonalizer<Elem = A>,
{
pub fn new(a: F, mut v: ArrayBase<S, Ix1>, mut ortho: Ortho) -> Self {
assert_eq!(ortho.len(), 0);
assert!(ortho.tolerance() < One::one());
let norm = v.norm_l2();
azip!((v in &mut v) *v = v.div_real(norm));
ortho.append(v.view());
Arnoldi {
a,
v,
ortho,
h: Vec::new(),
}
}
pub fn dim(&self) -> usize {
self.ortho.len()
}
pub fn complete(mut self) -> (Q<A>, H<A>) {
for _ in &mut self {}
let q = self.ortho.get_q();
let n = self.h.len();
let mut h = Array2::zeros((n, n).f());
for (i, hc) in self.h.iter().enumerate() {
let m = std::cmp::min(n, i + 2);
for j in 0..m {
h[(j, i)] = hc[j];
}
}
(q, h)
}
}
impl<A, S, F, Ortho> Iterator for Arnoldi<A, S, F, Ortho>
where
A: Scalar + Lapack,
S: DataMut<Elem = A>,
F: LinearOperator<Elem = A>,
Ortho: Orthogonalizer<Elem = A>,
{
type Item = Array1<A>;
fn next(&mut self) -> Option<Self::Item> {
self.a.apply_mut(&mut self.v);
let result = self.ortho.div_append(&mut self.v);
let norm = self.v.norm_l2();
azip!((v in &mut self.v) *v = v.div_real(norm));
match result {
AppendResult::Added(coef) => {
self.h.push(coef.clone());
Some(coef)
}
AppendResult::Dependent(coef) => {
self.h.push(coef);
None
}
}
}
}
pub fn arnoldi_householder<A, S>(a: impl LinearOperator<Elem = A>, v: ArrayBase<S, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
where
A: Scalar + Lapack,
S: DataMut<Elem = A>,
{
let householder = Householder::new(v.len(), tol);
Arnoldi::new(a, v, householder).complete()
}
pub fn arnoldi_mgs<A, S>(a: impl LinearOperator<Elem = A>, v: ArrayBase<S, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
where
A: Scalar + Lapack,
S: DataMut<Elem = A>,
{
let mgs = MGS::new(v.len(), tol);
Arnoldi::new(a, v, mgs).complete()
}