ndarray_linalg/krylov/
mgs.rs1use super::*;
4use crate::{generate::*, inner::*, norm::Norm};
5
6#[derive(Debug, Clone)]
8pub struct MGS<A: Scalar> {
9 dim: usize,
11
12 q: Vec<Array1<A>>,
14
15 tol: A::Real,
17}
18
19impl<A: Scalar + Lapack> MGS<A> {
20 pub fn new(dim: usize, tol: A::Real) -> Self {
22 Self {
23 dim,
24 q: Vec::new(),
25 tol,
26 }
27 }
28}
29
30impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
31 type Elem = A;
32
33 fn dim(&self) -> usize {
34 self.dim
35 }
36
37 fn len(&self) -> usize {
38 self.q.len()
39 }
40
41 fn tolerance(&self) -> A::Real {
42 self.tol
43 }
44
45 fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
46 where
47 S: DataMut<Elem = A>,
48 {
49 assert_eq!(a.len(), self.dim());
50 let mut coef = Array1::zeros(self.len() + 1);
51 for i in 0..self.len() {
52 let q = &self.q[i];
53 let c = q.inner(a);
54 azip!((a in &mut *a, &q in q) *a -= c * q);
55 coef[i] = c;
56 }
57 let nrm = a.norm_l2();
58 coef[self.len()] = A::from_real(nrm);
59 coef
60 }
61
62 fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
63 where
64 A: Lapack,
65 S: Data<Elem = A>,
66 {
67 let mut a = a.into_owned();
68 self.decompose(&mut a)
69 }
70
71 fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
72 where
73 A: Lapack,
74 S: Data<Elem = A>,
75 {
76 let mut a = a.into_owned();
77 self.div_append(&mut a)
78 }
79
80 fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
81 where
82 A: Lapack,
83 S: DataMut<Elem = A>,
84 {
85 let coef = self.decompose(a);
86 let nrm = coef[coef.len() - 1].re();
87 if nrm < self.tol {
88 return AppendResult::Dependent(coef);
90 }
91 azip!((a in &mut *a) *a /= A::from_real(nrm));
92 self.q.push(a.to_owned());
93 AppendResult::Added(coef)
94 }
95
96 fn get_q(&self) -> Q<A> {
97 hstack(&self.q).unwrap()
98 }
99}
100
101pub fn mgs<A, S>(
103 iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
104 dim: usize,
105 rtol: A::Real,
106 strategy: Strategy,
107) -> (Q<A>, R<A>)
108where
109 A: Scalar + Lapack,
110 S: Data<Elem = A>,
111{
112 let mgs = MGS::new(dim, rtol);
113 qr(iter, mgs, strategy)
114}