Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Arnoldi iterator #155

Merged
merged 15 commits into from
Jun 15, 2019
Prev Previous commit
Next Next commit
Add test
  • Loading branch information
termoshtt committed Jun 14, 2019
commit fe4efae2e2f7f796544d63fed9c508420689b06f
56 changes: 48 additions & 8 deletions src/krylov/arnoldi.rs
Original file line number Diff line number Diff line change
Expand Up @@ -54,10 +54,10 @@ where
pub fn complete(mut self) -> (Q<A>, H<A>) {
for _ in &mut self {} // execute iteration until convergent
let q = self.ortho.get_q();
let n = self.dim();
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::max(n, i + 1);
let m = std::cmp::min(n, i + 2);
for j in 0..m {
h[(j, i)] = hc[j];
}
Expand All @@ -79,12 +79,17 @@ where
(self.a)(&mut self.v);
let result = self.ortho.div_append(&mut self.v);
azip!(mut v(&mut self.v) in { *v = v.div_real(result.residual_norm()) });
if result.is_dependent() {
None
} else {
let coef = result.into_coeff();
self.h.push(coef.clone());
Some(coef)
match result {
AppendResult::Added(coef) => {
dbg!(&coef);
self.h.push(coef.clone());
Some(coef)
}
AppendResult::Dependent(coef) => {
dbg!(&coef);
self.h.push(coef);
None
}
}
}
}
Expand Down Expand Up @@ -126,3 +131,38 @@ where
let mgs = MGS::new(v.len(), tol);
Arnoldi::new(mul_mat(a), v, mgs).complete()
}

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

#[test]
fn aq_qh() {
let a = array![[1.0, 2.0], [3.0, 5.0]];
let mut v = Array::zeros(2);
v[0] = 1.0;
let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
println!("A = \n{:?}", &a);
println!("Q = \n{:?}", &q);
println!("H = \n{:?}", &h);
println!("AQ = \n{:?}", a.dot(&q));
println!("QH = \n{:?}", q.dot(&h));
panic!()
}

#[test]
fn aq_qh_random() {
let a: Array2<f64> = random((5, 5));
let mut v = Array::zeros(5);
v[0] = 1.0;
let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
println!("A = \n{:?}", &a);
println!("Q = \n{:?}", &q);
println!("H = \n{:?}", &h);
println!("AQ = \n{:?}", a.dot(&q));
println!("QH = \n{:?}", q.dot(&h));
panic!()
}

}