From 675e1391c23748b2e80bff892c0a1d2f6ad3525b Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 19:30:13 +0900 Subject: [PATCH 01/30] Test implementation of iterative QR decomposition --- src/arnoldi.rs | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + 2 files changed, 100 insertions(+) create mode 100644 src/arnoldi.rs diff --git a/src/arnoldi.rs b/src/arnoldi.rs new file mode 100644 index 00000000..76676e4f --- /dev/null +++ b/src/arnoldi.rs @@ -0,0 +1,99 @@ +use crate::{generate::*, norm::Norm, types::*}; +use ndarray::*; + +#[derive(Debug, Clone)] +pub struct MGS { + dim: usize, + q: Vec>, + r: Vec>, +} + +impl MGS { + pub fn new(dim: usize) -> Self { + Self { + dim, + q: Vec::new(), + r: Vec::new(), + } + } + + pub fn dim(&self) -> usize { + self.dim + } + + pub fn len(&self) -> usize { + self.q.len() + } + + pub fn append(&mut self, a: ArrayBase) -> A::Real + where + S: Data, + { + assert_eq!(a.len(), self.dim()); + let mut a = a.into_owned(); + let mut coef = Array1::zeros(self.len() + 1); + for i in 0..self.len() { + let q = &self.q[i]; + let c = a.dot(q); + azip!(mut a, q (q) in { *a = *a - c * q } ); + coef[i] = c; + } + let nrm = a.norm_l2(); + coef[self.len()] = A::from_real(nrm); + self.r.push(coef); + azip!(mut a in { *a = *a / A::from_real(nrm) }); + self.q.push(a); + nrm + } + + pub fn get_q(&self) -> Array2 { + hstack(&self.q).unwrap() + } + + pub fn get_r(&self) -> Array2 { + let len = self.len(); + let mut r = Array2::zeros((len, len)); + for i in 0..len { + for j in 0..=i { + r[(j, i)] = self.r[i][j]; + } + } + r + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + + const N: usize = 5; + + #[test] + fn new() { + let mgs: MGS = MGS::new(N); + assert_eq!(mgs.dim(), N); + assert_eq!(mgs.len(), 0); + } + + #[test] + fn append_random() { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = random((N, 3)); + dbg!(&a); + for col in a.axis_iter(Axis(1)) { + let res = mgs.append(col); + dbg!(res); + } + let q = mgs.get_q(); + dbg!(&q); + let r = mgs.get_r(); + dbg!(&r); + + dbg!(q.dot(&r)); + close_l2(&q.dot(&r), &a, 1e-9).unwrap(); + + dbg!(q.t().dot(&q)); + close_l2(&q.t().dot(&q), &Array2::eye(3), 1e-9).unwrap(); + } +} diff --git a/src/lib.rs b/src/lib.rs index 1a55a082..b1cc8b23 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,6 +36,7 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) +pub mod arnoldi; pub mod assert; pub mod cholesky; pub mod convert; From c26b20c69e09e3002e899a243d6105a8e14b45bd Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 19:42:18 +0900 Subject: [PATCH 02/30] Add test for complex --- src/arnoldi.rs | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 76676e4f..ecd30111 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -46,10 +46,12 @@ impl MGS { nrm } + /// Get orthogonal basis as Q matrix pub fn get_q(&self) -> Array2 { hstack(&self.q).unwrap() } + /// Get each vector norm and coefficients as R matrix pub fn get_r(&self) -> Array2 { let len = self.len(); let mut r = Array2::zeros((len, len)); @@ -66,6 +68,7 @@ impl MGS { mod tests { use super::*; use crate::assert::*; + use rand::{distributions::Standard, prelude::*}; const N: usize = 5; @@ -76,10 +79,12 @@ mod tests { assert_eq!(mgs.len(), 0); } - #[test] - fn append_random() { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = random((N, 3)); + fn test(rtol: A::Real) + where + Standard: Distribution, + { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); for col in a.axis_iter(Axis(1)) { let res = mgs.append(col); @@ -91,9 +96,19 @@ mod tests { dbg!(&r); dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a, 1e-9).unwrap(); + close_l2(&q.dot(&r), &a, rtol).unwrap(); dbg!(q.t().dot(&q)); - close_l2(&q.t().dot(&q), &Array2::eye(3), 1e-9).unwrap(); + close_l2(&q.t().dot(&q), &Array2::eye(3), rtol).unwrap(); + } + + #[test] + fn test_f32() { + test::(1e-5); + } + + #[test] + fn test_c32() { + test::(1e-5); } } From 78ad4e784ae7b1ae5c5a86874314abcb95ef042f Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 20:26:20 +0900 Subject: [PATCH 03/30] Drop rand restriction (#143) --- src/arnoldi.rs | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index ecd30111..d1cbc2a3 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -68,7 +68,6 @@ impl MGS { mod tests { use super::*; use crate::assert::*; - use rand::{distributions::Standard, prelude::*}; const N: usize = 5; @@ -79,10 +78,7 @@ mod tests { assert_eq!(mgs.len(), 0); } - fn test(rtol: A::Real) - where - Standard: Distribution, - { + fn test(rtol: A::Real) { let mut mgs: MGS = MGS::new(N); let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); From 28617dbf0f4f89fd6bcb673f5b518a762aa3f262 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 21:00:05 +0900 Subject: [PATCH 04/30] Use conjugate matrix instead of transpose for test --- src/arnoldi.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index d1cbc2a3..e880198f 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -94,8 +94,9 @@ mod tests { dbg!(q.dot(&r)); close_l2(&q.dot(&r), &a, rtol).unwrap(); - dbg!(q.t().dot(&q)); - close_l2(&q.t().dot(&q), &Array2::eye(3), rtol).unwrap(); + let qt: Array2<_> = conjugate(&q); + dbg!(qt.dot(&q)); + close_l2(&qt.dot(&q), &Array2::eye(3), rtol).unwrap(); } #[test] From 22cb2d286fc6564833622c164dc235d4a7700cc0 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 29 Apr 2019 23:22:35 +0900 Subject: [PATCH 05/30] Add inner product --- src/inner.rs | 21 +++++++++++++++++++++ src/lib.rs | 1 + 2 files changed, 22 insertions(+) create mode 100644 src/inner.rs diff --git a/src/inner.rs b/src/inner.rs new file mode 100644 index 00000000..299e9669 --- /dev/null +++ b/src/inner.rs @@ -0,0 +1,21 @@ +use crate::types::*; +use ndarray::*; + +pub trait Inner { + /// Inner product `(self.conjugate, rhs)` + fn inner(&self, rhs: &ArrayBase) -> S::Elem; +} + +impl Inner for ArrayBase +where + A: Scalar, + S1: Data, + S2: Data, +{ + fn inner(&self, rhs: &ArrayBase) -> A { + Zip::from(self) + .and(rhs) + .fold_while(A::zero(), |acc, s, r| FoldWhile::Continue(acc + s.conj() * *r)) + .into_inner() + } +} diff --git a/src/lib.rs b/src/lib.rs index b1cc8b23..68367706 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -44,6 +44,7 @@ pub mod diagonal; pub mod eigh; pub mod error; pub mod generate; +pub mod inner; pub mod lapack; pub mod layout; pub mod norm; From f40f1300998951ae30e9bb709a5aae412ba5f353 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 29 Apr 2019 23:22:48 +0900 Subject: [PATCH 06/30] Use inner product --- src/arnoldi.rs | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index e880198f..0b7e2051 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,4 +1,4 @@ -use crate::{generate::*, norm::Norm, types::*}; +use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; #[derive(Debug, Clone)] @@ -34,7 +34,7 @@ impl MGS { let mut coef = Array1::zeros(self.len() + 1); for i in 0..self.len() { let q = &self.q[i]; - let c = a.dot(q); + let c = q.inner(&a); azip!(mut a, q (q) in { *a = *a - c * q } ); coef[i] = c; } @@ -108,4 +108,14 @@ mod tests { fn test_c32() { test::(1e-5); } + + #[test] + fn test_f64() { + test::(1e-9); + } + + #[test] + fn test_c64() { + test::(1e-9); + } } From edcb8b4f5bc7011b7d794335a37e1bfe0d52e504 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 00:10:00 +0900 Subject: [PATCH 07/30] append_if --- src/arnoldi.rs | 65 ++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 13 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 0b7e2051..3638eab7 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,5 +1,6 @@ use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; +use num_traits::Zero; #[derive(Debug, Clone)] pub struct MGS { @@ -8,6 +9,11 @@ pub struct MGS { r: Vec>, } +pub enum Dependence { + Orthogonal(A::Real), + Coefficient(Array1), +} + impl MGS { pub fn new(dim: usize) -> Self { Self { @@ -25,7 +31,16 @@ impl MGS { self.q.len() } + /// Add new vector, return residual norm pub fn append(&mut self, a: ArrayBase) -> A::Real + where + S: Data, + { + self.append_if(a, A::Real::zero()).unwrap() + } + + /// Add new vector if the residual is larger than relative tolerance + pub fn append_if(&mut self, a: ArrayBase, rtol: A::Real) -> Option where S: Data, { @@ -39,11 +54,14 @@ impl MGS { coef[i] = c; } let nrm = a.norm_l2(); + if nrm < rtol { + return None; + } coef[self.len()] = A::from_real(nrm); self.r.push(coef); azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - nrm + Some(nrm) } /// Get orthogonal basis as Q matrix @@ -78,7 +96,7 @@ mod tests { assert_eq!(mgs.len(), 0); } - fn test(rtol: A::Real) { + fn test_append(rtol: A::Real) { let mut mgs: MGS = MGS::new(N); let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); @@ -100,22 +118,43 @@ mod tests { } #[test] - fn test_f32() { - test::(1e-5); + fn append() { + test_append::(1e-5); + test_append::(1e-5); + test_append::(1e-9); + test_append::(1e-9); } - #[test] - fn test_c32() { - test::(1e-5); - } + fn test_append_if(rtol: A::Real) { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = crate::generate::random((N, 8)); + dbg!(&a); + for col in a.axis_iter(Axis(1)) { + match mgs.append_if(col, rtol) { + Some(res) => { + dbg!(res); + } + None => break, + } + } + let q = mgs.get_q(); + dbg!(&q); + let r = mgs.get_r(); + dbg!(&r); - #[test] - fn test_f64() { - test::(1e-9); + dbg!(q.dot(&r)); + close_l2(&q.dot(&r), &a.slice(s![.., 0..N]), rtol).unwrap(); + + let qt: Array2<_> = conjugate(&q); + dbg!(qt.dot(&q)); + close_l2(&qt.dot(&q), &Array2::eye(N), rtol).unwrap(); } #[test] - fn test_c64() { - test::(1e-9); + fn append_if() { + test_append_if::(1e-5); + test_append_if::(1e-5); + test_append_if::(1e-9); + test_append_if::(1e-9); } } From 7c08005908ed735b7c7e52772e10981d1ab02514 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 18:29:13 +0900 Subject: [PATCH 08/30] Do not manage R matrix in MGS --- src/arnoldi.rs | 72 ++++++++++++++++++++++++-------------------------- 1 file changed, 34 insertions(+), 38 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 3638eab7..c0f69cba 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,51 +1,44 @@ use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; -use num_traits::Zero; +/// Iterative orthogonalizer using modified Gram-Schmit procedure #[derive(Debug, Clone)] pub struct MGS { - dim: usize, + /// Dimension of base space + dimension: usize, + /// Basis of spanned space q: Vec>, - r: Vec>, } -pub enum Dependence { - Orthogonal(A::Real), - Coefficient(Array1), -} +pub type Residual = ArrayBase; +pub type Coefficient = Array1; impl MGS { - pub fn new(dim: usize) -> Self { + pub fn new(dimension: usize) -> Self { Self { - dim, + dimension, q: Vec::new(), - r: Vec::new(), } } pub fn dim(&self) -> usize { - self.dim + self.dimension } pub fn len(&self) -> usize { self.q.len() } - /// Add new vector, return residual norm - pub fn append(&mut self, a: ArrayBase) -> A::Real - where - S: Data, - { - self.append_if(a, A::Real::zero()).unwrap() - } - - /// Add new vector if the residual is larger than relative tolerance - pub fn append_if(&mut self, a: ArrayBase, rtol: A::Real) -> Option + /// Orthogonalize given vector using current basis + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) where - S: Data, + S: DataMut, { assert_eq!(a.len(), self.dim()); - let mut a = a.into_owned(); let mut coef = Array1::zeros(self.len() + 1); for i in 0..self.len() { let q = &self.q[i]; @@ -54,32 +47,35 @@ impl MGS { coef[i] = c; } let nrm = a.norm_l2(); + coef[self.len()] = A::from_real(nrm); + (a, coef) + } + + /// Add new vector if the residual is larger than relative tolerance + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + pub fn append_if_independent(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + where + S: Data, + { + let a = a.into_owned(); + let (mut a, coef) = self.orthogonalize(a); + let nrm = coef[coef.len()].re(); if nrm < rtol { + // Linearly dependent return None; } - coef[self.len()] = A::from_real(nrm); - self.r.push(coef); azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - Some(nrm) + Some(coef) } /// Get orthogonal basis as Q matrix pub fn get_q(&self) -> Array2 { hstack(&self.q).unwrap() } - - /// Get each vector norm and coefficients as R matrix - pub fn get_r(&self) -> Array2 { - let len = self.len(); - let mut r = Array2::zeros((len, len)); - for i in 0..len { - for j in 0..=i { - r[(j, i)] = self.r[i][j]; - } - } - r - } } #[cfg(test)] From c9dce235b575bcc4ee543d34f142dc85bdc67492 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 19:33:20 +0900 Subject: [PATCH 09/30] Rewrite tests for MGS::append --- src/arnoldi.rs | 126 ++++++++++++++++++++----------------------------- 1 file changed, 51 insertions(+), 75 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index c0f69cba..54ae6478 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -10,10 +10,25 @@ pub struct MGS { q: Vec>, } +/// Residual vector of orthogonalization pub type Residual = ArrayBase; +/// Residual vector of orthogonalization pub type Coefficient = Array1; +/// Q-matrix (unitary matrix) +pub type Q = Array2; +/// R-matrix (upper triangle) +pub type R = Array2; -impl MGS { +impl MGS { + /// Create empty linear space + /// + /// ```rust + /// # use ndarray_linalg::*; + /// const N: usize = 5; + /// let mgs = arnoldi::MGS::::new(N); + /// assert_eq!(mgs.dim(), N); + /// assert_eq!(mgs.len(), 0); + /// ``` pub fn new(dimension: usize) -> Self { Self { dimension, @@ -36,6 +51,7 @@ impl MGS { /// - if the size of the input array mismaches to the dimension pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) where + A: Lapack, S: DataMut, { assert_eq!(a.len(), self.dim()); @@ -56,13 +72,27 @@ impl MGS { /// Panic /// ------- /// - if the size of the input array mismaches to the dimension - pub fn append_if_independent(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + /// + /// ```rust + /// # use ndarray::*; + /// # use ndarray_linalg::*; + /// let mut mgs = arnoldi::MGS::new(3); + /// let coef = mgs.append(array![1.0, 0.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); + /// + /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + /// + /// assert!(mgs.append(array![1.0, 1.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// ``` + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where + A: Lapack, S: Data, { let a = a.into_owned(); let (mut a, coef) = self.orthogonalize(a); - let nrm = coef[coef.len()].re(); + let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent return None; @@ -73,84 +103,30 @@ impl MGS { } /// Get orthogonal basis as Q matrix - pub fn get_q(&self) -> Array2 { + pub fn get_q(&self) -> Q { hstack(&self.q).unwrap() } } -#[cfg(test)] -mod tests { - use super::*; - use crate::assert::*; - - const N: usize = 5; - - #[test] - fn new() { - let mgs: MGS = MGS::new(N); - assert_eq!(mgs.dim(), N); - assert_eq!(mgs.len(), 0); - } - - fn test_append(rtol: A::Real) { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = crate::generate::random((N, 3)); - dbg!(&a); - for col in a.axis_iter(Axis(1)) { - let res = mgs.append(col); - dbg!(res); +pub fn mgs(iter: impl Iterator>, dim: usize, rtol: A::Real) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + let mut ortho = MGS::new(dim); + let mut coefs = Vec::new(); + for a in iter { + match ortho.append(a, rtol) { + Some(coef) => coefs.push(coef), + None => break, } - let q = mgs.get_q(); - dbg!(&q); - let r = mgs.get_r(); - dbg!(&r); - - dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a, rtol).unwrap(); - - let qt: Array2<_> = conjugate(&q); - dbg!(qt.dot(&q)); - close_l2(&qt.dot(&q), &Array2::eye(3), rtol).unwrap(); - } - - #[test] - fn append() { - test_append::(1e-5); - test_append::(1e-5); - test_append::(1e-9); - test_append::(1e-9); } - - fn test_append_if(rtol: A::Real) { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = crate::generate::random((N, 8)); - dbg!(&a); - for col in a.axis_iter(Axis(1)) { - match mgs.append_if(col, rtol) { - Some(res) => { - dbg!(res); - } - None => break, - } + let m = coefs.len(); + let mut r = Array2::zeros((m, m)); + for i in 0..m { + for j in 0..=i { + r[(j, i)] = coefs[i][j]; } - let q = mgs.get_q(); - dbg!(&q); - let r = mgs.get_r(); - dbg!(&r); - - dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a.slice(s![.., 0..N]), rtol).unwrap(); - - let qt: Array2<_> = conjugate(&q); - dbg!(qt.dot(&q)); - close_l2(&qt.dot(&q), &Array2::eye(N), rtol).unwrap(); - } - - #[test] - fn append_if() { - test_append_if::(1e-5); - test_append_if::(1e-5); - test_append_if::(1e-9); - test_append_if::(1e-9); } + (ortho.get_q(), r) } From d61793f7c38ceb4396ed86ede899dc513e0df63d Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 19:42:04 +0900 Subject: [PATCH 10/30] Fix test case --- src/arnoldi.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 54ae6478..3a7baeea 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -83,7 +83,7 @@ impl MGS { /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); /// - /// assert!(mgs.append(array![1.0, 1.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector /// ``` pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where From 93118f598173cab6403f888488aab9b38efdbf84 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Fri, 3 May 2019 18:26:11 +0900 Subject: [PATCH 11/30] orthogonalize takes &mut --- src/arnoldi.rs | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 3a7baeea..9f880021 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -10,11 +10,7 @@ pub struct MGS { q: Vec>, } -/// Residual vector of orthogonalization -pub type Residual = ArrayBase; -/// Residual vector of orthogonalization -pub type Coefficient = Array1; -/// Q-matrix (unitary matrix) +/// Q-matrix (unitary) pub type Q = Array2; /// R-matrix (upper triangle) pub type R = Array2; @@ -49,7 +45,7 @@ impl MGS { /// Panic /// ------- /// - if the size of the input array mismaches to the dimension - pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) + pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, S: DataMut, @@ -59,12 +55,12 @@ impl MGS { for i in 0..self.len() { let q = &self.q[i]; let c = q.inner(&a); - azip!(mut a, q (q) in { *a = *a - c * q } ); + azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } ); coef[i] = c; } let nrm = a.norm_l2(); coef[self.len()] = A::from_real(nrm); - (a, coef) + coef } /// Add new vector if the residual is larger than relative tolerance @@ -85,13 +81,13 @@ impl MGS { /// /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector /// ``` - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where A: Lapack, S: Data, { - let a = a.into_owned(); - let (mut a, coef) = self.orthogonalize(a); + let mut a = a.into_owned(); + let coef = self.orthogonalize(&mut a); let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent From 0bc6565efa73c17504b4ebe61a643f546d40f536 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Fri, 3 May 2019 18:47:40 +0900 Subject: [PATCH 12/30] Strategy for linearly dependent vectors --- src/arnoldi.rs | 51 ++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 8 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 9f880021..d9519b17 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -73,15 +73,19 @@ impl MGS { /// # use ndarray::*; /// # use ndarray_linalg::*; /// let mut mgs = arnoldi::MGS::new(3); - /// let coef = mgs.append(array![1.0, 0.0, 0.0], 1e-9).unwrap(); + /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); /// /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); /// - /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend + /// + /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { + /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector + /// } /// ``` - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where A: Lapack, S: Data, @@ -91,11 +95,11 @@ impl MGS { let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent - return None; + return Err(coef); } azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - Some(coef) + Ok(coef) } /// Get orthogonal basis as Q matrix @@ -104,7 +108,34 @@ impl MGS { } } -pub fn mgs(iter: impl Iterator>, dim: usize, rtol: A::Real) -> (Q, R) +/// Strategy for linearly dependent vectors appearing in iterative QR decomposition +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum Strategy { + /// Terminate iteration if dependent vector comes + Terminate, + + /// Skip dependent vector + Skip, + + /// Orghotonalize dependent vector without adding to Q, + /// thus R must be non-regular like following: + /// + /// ```ignore + /// x x x x x + /// 0 x x x x + /// 0 0 0 x x + /// 0 0 0 0 x + /// 0 0 0 0 0 // 0-filled to be square matrix + /// ``` + Full, +} + +pub fn mgs( + iter: impl Iterator>, + dim: usize, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) where A: Scalar + Lapack, S: Data, @@ -113,8 +144,12 @@ where let mut coefs = Vec::new(); for a in iter { match ortho.append(a, rtol) { - Some(coef) => coefs.push(coef), - None => break, + Ok(coef) => coefs.push(coef), + Err(coef) => match strategy { + Strategy::Terminate => break, + Strategy::Skip => continue, + Strategy::Full => coefs.push(coef), + }, } } let m = coefs.len(); From ee59945cb523a2b681cc4140461c5318ec3e9aa6 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 4 May 2019 00:13:17 +0900 Subject: [PATCH 13/30] Add tests of msg --- src/arnoldi.rs | 12 ++++--- tests/arnoldi.rs | 83 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 4 deletions(-) create mode 100644 tests/arnoldi.rs diff --git a/src/arnoldi.rs b/src/arnoldi.rs index d9519b17..1067c10d 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -130,6 +130,7 @@ pub enum Strategy { Full, } +/// Online QR decomposition of vectors using modified Gram-Schmit algorithm pub fn mgs( iter: impl Iterator>, dim: usize, @@ -152,11 +153,14 @@ where }, } } + let n = ortho.len(); let m = coefs.len(); - let mut r = Array2::zeros((m, m)); - for i in 0..m { - for j in 0..=i { - r[(j, i)] = coefs[i][j]; + let mut r = Array2::zeros((n, m).f()); + for j in 0..m { + for i in 0..n { + if i < coefs[j].len() { + r[(i, j)] = coefs[j][i]; + } } } (ortho.get_q(), r) diff --git a/tests/arnoldi.rs b/tests/arnoldi.rs new file mode 100644 index 00000000..26040094 --- /dev/null +++ b/tests/arnoldi.rs @@ -0,0 +1,83 @@ +use ndarray::*; +use ndarray_linalg::{arnoldi::*, *}; + +fn qr_full() { + const N: usize = 5; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn qr_full_real() { + qr_full::(); +} + +#[test] +fn qr_full_complex() { + qr_full::(); +} + +fn qr() { + const N: usize = 4; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N / 2)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); +} + +#[test] +fn qr_real() { + qr::(); +} + +#[test] +fn qr_complex() { + qr::(); +} + +fn qr_over() { + const N: usize = 4; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Skip + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Full + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn qr_over_real() { + qr_over::(); +} + +#[test] +fn qr_over_complex() { + qr_over::(); +} From b4b8df12a1cf7d69d57162de4e9402873f5093b0 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 01:52:54 +0900 Subject: [PATCH 14/30] Householder reflection --- src/arnoldi.rs | 9 ++--- src/householder.rs | 84 ++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + 3 files changed, 90 insertions(+), 4 deletions(-) create mode 100644 src/householder.rs diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 1067c10d..6f55c830 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -65,10 +65,6 @@ impl MGS { /// Add new vector if the residual is larger than relative tolerance /// - /// Panic - /// ------- - /// - if the size of the input array mismaches to the dimension - /// /// ```rust /// # use ndarray::*; /// # use ndarray_linalg::*; @@ -85,6 +81,11 @@ impl MGS { /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector /// } /// ``` + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + /// pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where A: Lapack, diff --git a/src/householder.rs b/src/householder.rs new file mode 100644 index 00000000..6ef0bd2e --- /dev/null +++ b/src/householder.rs @@ -0,0 +1,84 @@ +use crate::{inner::*, norm::*, types::*}; +use ndarray::*; + +/// Iterative orthogonalizer using Householder reflection +#[derive(Debug, Clone)] +pub struct Householder { + dim: usize, + v: Vec>, +} + +impl Householder { + pub fn new(dim: usize) -> Self { + Householder { dim, v: Vec::new() } + } + + pub fn len(&self) -> usize { + self.v.len() + } + + /// Take a Reflection `P = I - 2ww^T` + fn reflect>(&self, k: usize, a: &mut ArrayBase) { + assert!(k < self.v.len()); + assert_eq!(a.len(), self.dim); + let w = self.v[k].slice(s![k..]); + let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..])); + for l in k..self.dim { + a[l] -= c * w[l]; + } + } + + /// Orghotonalize input vector by reflectors + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + pub fn orthogonalize(&self, a: &mut ArrayBase) -> A::Real + where + A: Lapack, + S: DataMut, + { + for k in 0..self.len() { + self.reflect(k, a); + } + // residual norm + a.slice(s![self.len()..]).norm_l2() + } + + /// Orghotonalize input vector by reflectors + /// + /// ```rust + /// # use ndarray::*; + /// # use ndarray_linalg::*; + /// let mut mgs = arnoldi::MGS::new(3); + /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); + /// + /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + /// + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend + /// + /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { + /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector + /// } + /// ``` + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + /// + pub fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> + where + A: Lapack, + S: DataMut, + { + let residual = self.orthogonalize(&mut a); + let coef = a.slice(s![..self.len()]).into_owned(); + if residual < rtol { + return Err(coef); + } + self.v.push(a.into_owned()); + Ok(coef) + } +} diff --git a/src/lib.rs b/src/lib.rs index 68367706..30fcb053 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -44,6 +44,7 @@ pub mod diagonal; pub mod eigh; pub mod error; pub mod generate; +pub mod householder; pub mod inner; pub mod lapack; pub mod layout; From 86d25d8062830e21d5cd087d6fb49fde1fdf2597 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 02:49:52 +0900 Subject: [PATCH 15/30] Split into krylov submodule --- src/krylov/householder.rs | 70 ++++++++++++++++++++++ src/krylov/mgs.rs | 78 ++++++++++++++++++++++++ src/{arnoldi.rs => krylov/mod.rs} | 99 +++++++++---------------------- src/lib.rs | 3 +- tests/{arnoldi.rs => krylov.rs} | 18 +++--- 5 files changed, 186 insertions(+), 82 deletions(-) create mode 100644 src/krylov/householder.rs create mode 100644 src/krylov/mgs.rs rename src/{arnoldi.rs => krylov/mod.rs} (55%) rename tests/{arnoldi.rs => krylov.rs} (76%) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs new file mode 100644 index 00000000..8acd10f4 --- /dev/null +++ b/src/krylov/householder.rs @@ -0,0 +1,70 @@ +use super::*; +use crate::{inner::*, norm::*}; + +/// Iterative orthogonalizer using Householder reflection +#[derive(Debug, Clone)] +pub struct Householder { + dim: usize, + v: Vec>, +} + +impl Householder { + pub fn new(dim: usize) -> Self { + Householder { dim, v: Vec::new() } + } + + /// Take a Reflection `P = I - 2ww^T` + fn reflect>(&self, k: usize, a: &mut ArrayBase) { + assert!(k < self.v.len()); + assert_eq!(a.len(), self.dim); + let w = self.v[k].slice(s![k..]); + let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..])); + for l in k..self.dim { + a[l] -= c * w[l]; + } + } +} + +impl Orthogonalizer for Householder { + type Elem = A; + + fn new(dim: usize) -> Self { + Self { dim, v: Vec::new() } + } + + fn dim(&self) -> usize { + self.dim + } + + fn len(&self) -> usize { + self.v.len() + } + + fn orthogonalize(&self, a: &mut ArrayBase) -> A::Real + where + S: DataMut, + { + for k in 0..self.len() { + self.reflect(k, a); + } + // residual norm + a.slice(s![self.len()..]).norm_l2() + } + + fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> + where + S: DataMut, + { + let residual = self.orthogonalize(&mut a); + let coef = a.slice(s![..self.len()]).into_owned(); + if residual < rtol { + return Err(coef); + } + self.v.push(a.into_owned()); + Ok(coef) + } + + fn get_q(&self) -> Q { + unimplemented!() + } +} diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs new file mode 100644 index 00000000..a6d74c13 --- /dev/null +++ b/src/krylov/mgs.rs @@ -0,0 +1,78 @@ +use super::*; +use crate::{generate::*, inner::*, norm::Norm}; + +/// Iterative orthogonalizer using modified Gram-Schmit procedure +#[derive(Debug, Clone)] +pub struct MGS { + /// Dimension of base space + dim: usize, + /// Basis of spanned space + q: Vec>, +} + +impl MGS { + fn ortho(&self, a: &mut ArrayBase) -> Array1 + where + A: Lapack, + S: DataMut, + { + assert_eq!(a.len(), self.dim); + let mut coef = Array1::zeros(self.q.len() + 1); + for i in 0..self.q.len() { + let q = &self.q[i]; + let c = q.inner(&a); + azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } ); + coef[i] = c; + } + let nrm = a.norm_l2(); + coef[self.q.len()] = A::from(nrm).unwrap(); + coef + } +} + +impl Orthogonalizer for MGS { + type Elem = A; + + fn new(dim: usize) -> Self { + Self { dim, q: Vec::new() } + } + + fn dim(&self) -> usize { + self.dim + } + + fn len(&self) -> usize { + self.q.len() + } + + fn orthogonalize(&self, a: &mut ArrayBase) -> A::Real + where + S: DataMut, + { + let coef = self.ortho(a); + // Write coefficients into `a` + azip!(mut a (a.slice_mut(s![0..self.len()])), coef in { *a = coef }); + // 0-fill for remaining + azip!(mut a (a.slice_mut(s![self.len()..])) in { *a = A::zero() }); + coef[self.len()].re() + } + + fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> + where + S: DataMut, + { + let coef = self.ortho(&mut a); + let nrm = coef[self.len()].re(); + if nrm < rtol { + // Linearly dependent + return Err(coef); + } + azip!(mut a in { *a = *a / A::from_real(nrm) }); + self.q.push(a.into_owned()); + Ok(coef) + } + + fn get_q(&self) -> Q { + hstack(&self.q).unwrap() + } +} diff --git a/src/arnoldi.rs b/src/krylov/mod.rs similarity index 55% rename from src/arnoldi.rs rename to src/krylov/mod.rs index 6f55c830..926db023 100644 --- a/src/arnoldi.rs +++ b/src/krylov/mod.rs @@ -1,74 +1,42 @@ -use crate::{generate::*, inner::*, norm::Norm, types::*}; +use crate::types::*; use ndarray::*; -/// Iterative orthogonalizer using modified Gram-Schmit procedure -#[derive(Debug, Clone)] -pub struct MGS { - /// Dimension of base space - dimension: usize, - /// Basis of spanned space - q: Vec>, -} +mod householder; +mod mgs; + +pub use householder::Householder; +pub use mgs::MGS; /// Q-matrix (unitary) pub type Q = Array2; /// R-matrix (upper triangle) pub type R = Array2; -impl MGS { - /// Create empty linear space - /// - /// ```rust - /// # use ndarray_linalg::*; - /// const N: usize = 5; - /// let mgs = arnoldi::MGS::::new(N); - /// assert_eq!(mgs.dim(), N); - /// assert_eq!(mgs.len(), 0); - /// ``` - pub fn new(dimension: usize) -> Self { - Self { - dimension, - q: Vec::new(), - } - } +pub trait Orthogonalizer { + type Elem: Scalar; - pub fn dim(&self) -> usize { - self.dimension - } + /// Create empty linear space + fn new(dim: usize) -> Self; - pub fn len(&self) -> usize { - self.q.len() - } + fn dim(&self) -> usize; + fn len(&self) -> usize; - /// Orthogonalize given vector using current basis + /// Orthogonalize given vector /// /// Panic /// ------- /// - if the size of the input array mismaches to the dimension - pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + /// + fn orthogonalize(&self, a: &mut ArrayBase) -> ::Real where - A: Lapack, - S: DataMut, - { - assert_eq!(a.len(), self.dim()); - let mut coef = Array1::zeros(self.len() + 1); - for i in 0..self.len() { - let q = &self.q[i]; - let c = q.inner(&a); - azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } ); - coef[i] = c; - } - let nrm = a.norm_l2(); - coef[self.len()] = A::from_real(nrm); - coef - } + S: DataMut; /// Add new vector if the residual is larger than relative tolerance /// /// ```rust /// # use ndarray::*; - /// # use ndarray_linalg::*; - /// let mut mgs = arnoldi::MGS::new(3); + /// # use ndarray_linalg::{*, krylov::*}; + /// let mut mgs = krylov::MGS::new(3); /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); /// @@ -86,27 +54,16 @@ impl MGS { /// ------- /// - if the size of the input array mismaches to the dimension /// - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> + fn append( + &mut self, + a: ArrayBase, + rtol: ::Real, + ) -> Result, Array1> where - A: Lapack, - S: Data, - { - let mut a = a.into_owned(); - let coef = self.orthogonalize(&mut a); - let nrm = coef[coef.len() - 1].re(); - if nrm < rtol { - // Linearly dependent - return Err(coef); - } - azip!(mut a in { *a = *a / A::from_real(nrm) }); - self.q.push(a); - Ok(coef) - } + S: DataMut; - /// Get orthogonal basis as Q matrix - pub fn get_q(&self) -> Q { - hstack(&self.q).unwrap() - } + /// Get Q-matrix of generated basis + fn get_q(&self) -> Q; } /// Strategy for linearly dependent vectors appearing in iterative QR decomposition @@ -132,7 +89,7 @@ pub enum Strategy { } /// Online QR decomposition of vectors using modified Gram-Schmit algorithm -pub fn mgs( +pub fn qr( iter: impl Iterator>, dim: usize, rtol: A::Real, @@ -145,7 +102,7 @@ where let mut ortho = MGS::new(dim); let mut coefs = Vec::new(); for a in iter { - match ortho.append(a, rtol) { + match ortho.append(a.into_owned(), rtol) { Ok(coef) => coefs.push(coef), Err(coef) => match strategy { Strategy::Terminate => break, diff --git a/src/lib.rs b/src/lib.rs index 30fcb053..ed4d910d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,7 +36,6 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) -pub mod arnoldi; pub mod assert; pub mod cholesky; pub mod convert; @@ -44,8 +43,8 @@ pub mod diagonal; pub mod eigh; pub mod error; pub mod generate; -pub mod householder; pub mod inner; +pub mod krylov; pub mod lapack; pub mod layout; pub mod norm; diff --git a/tests/arnoldi.rs b/tests/krylov.rs similarity index 76% rename from tests/arnoldi.rs rename to tests/krylov.rs index 26040094..0bf3f26c 100644 --- a/tests/arnoldi.rs +++ b/tests/krylov.rs @@ -1,12 +1,12 @@ use ndarray::*; -use ndarray_linalg::{arnoldi::*, *}; +use ndarray_linalg::{krylov::*, *}; fn qr_full() { const N: usize = 5; let rtol: A::Real = A::real(1e-9); let a: Array2 = random((N, N)); - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); @@ -23,12 +23,12 @@ fn qr_full_complex() { qr_full::(); } -fn qr() { +fn qr_() { const N: usize = 4; let rtol: A::Real = A::real(1e-9); let a: Array2 = random((N, N / 2)); - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); @@ -37,12 +37,12 @@ fn qr() { #[test] fn qr_real() { - qr::(); + qr_::(); } #[test] fn qr_complex() { - qr::(); + qr_::(); } fn qr_over() { @@ -52,21 +52,21 @@ fn qr_over() { let a: Array2 = random((N, N * 2)); // Terminate - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); let a_sub = a.slice(s![.., 0..N]); assert_close_l2!(&q.dot(&r), &a_sub, rtol); let qc: Array2 = conjugate(&q); assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); // Skip - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); let a_sub = a.slice(s![.., 0..N]); assert_close_l2!(&q.dot(&r), &a_sub, rtol); let qc: Array2 = conjugate(&q); assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); // Full - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); From ac8fa2b82e288578c6c5f18d96301ed666720376 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 17:14:39 +0900 Subject: [PATCH 16/30] Add test for householder (and bug found) --- src/krylov/mgs.rs | 23 ++++++ src/krylov/mod.rs | 37 +++++++++- tests/krylov.rs | 179 ++++++++++++++++++++++++++++++---------------- 3 files changed, 173 insertions(+), 66 deletions(-) diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index a6d74c13..abb8ba0a 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -76,3 +76,26 @@ impl Orthogonalizer for MGS { hstack(&self.q).unwrap() } } + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + + #[test] + fn mgs_append() { + let mut mgs = MGS::new(3); + let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0], 1e-9).unwrap(); + + let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + + assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); + + if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { + close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); + } + } + +} diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 926db023..2006e8d9 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -88,10 +88,10 @@ pub enum Strategy { Full, } -/// Online QR decomposition of vectors using modified Gram-Schmit algorithm +/// Online QR decomposition using arbitary orthogonalizer pub fn qr( iter: impl Iterator>, - dim: usize, + mut ortho: impl Orthogonalizer, rtol: A::Real, strategy: Strategy, ) -> (Q, R) @@ -99,7 +99,8 @@ where A: Scalar + Lapack, S: Data, { - let mut ortho = MGS::new(dim); + assert_eq!(ortho.len(), 0); + let mut coefs = Vec::new(); for a in iter { match ortho.append(a.into_owned(), rtol) { @@ -123,3 +124,33 @@ where } (ortho.get_q(), r) } + +/// Online QR decomposition using modified Gram-Schmit +pub fn mgs( + iter: impl Iterator>, + dim: usize, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + let mgs = MGS::new(dim); + qr(iter, mgs, rtol, strategy) +} + +/// Online QR decomposition using modified Gram-Schmit +pub fn householder( + iter: impl Iterator>, + dim: usize, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + let h = Householder::new(dim); + qr(iter, h, rtol, strategy) +} diff --git a/tests/krylov.rs b/tests/krylov.rs index 0bf3f26c..2f6f5f82 100644 --- a/tests/krylov.rs +++ b/tests/krylov.rs @@ -1,83 +1,136 @@ use ndarray::*; use ndarray_linalg::{krylov::*, *}; -fn qr_full() { - const N: usize = 5; - let rtol: A::Real = A::real(1e-9); - - let a: Array2 = random((N, N)); - let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); - - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); -} - #[test] -fn qr_full_real() { - qr_full::(); +fn mgs_full() { + fn test(rtol: A::Real) { + const N: usize = 5; + let a: Array2 = random((N, N)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } #[test] -fn qr_full_complex() { - qr_full::(); -} - -fn qr_() { - const N: usize = 4; - let rtol: A::Real = A::real(1e-9); - - let a: Array2 = random((N, N / 2)); - let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); - - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); +fn mgs_half() { + fn test(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N / 2)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } #[test] -fn qr_real() { - qr_::(); +fn mgs_over() { + fn test(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Skip + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Full + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } #[test] -fn qr_complex() { - qr_::(); -} - -fn qr_over() { - const N: usize = 4; - let rtol: A::Real = A::real(1e-9); - - let a: Array2 = random((N, N * 2)); - - // Terminate - let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - - // Skip - let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); - let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - - // Full - let (q, r) = qr(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); - assert_close_l2!(&q.dot(&r), &a, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +fn householder_full() { + fn test(rtol: A::Real) { + const N: usize = 5; + let a: Array2 = random((N, N)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } #[test] -fn qr_over_real() { - qr_over::(); +fn householder_half() { + fn test(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N / 2)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } #[test] -fn qr_over_complex() { - qr_over::(); +fn householder_over() { + fn test(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Skip + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Full + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + } + + test::(1e-5); + test::(1e-9); + test::(1e-5); + test::(1e-9); } From ff690dfd2fd6f3233f330436e0c3b8d889e3539f Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 17:50:06 +0900 Subject: [PATCH 17/30] Fix full-rank case --- src/krylov/householder.rs | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 8acd10f4..d05a260b 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -1,5 +1,6 @@ use super::*; use crate::{inner::*, norm::*}; +use num_traits::Zero; /// Iterative orthogonalizer using Householder reflection #[derive(Debug, Clone)] @@ -47,6 +48,10 @@ impl Orthogonalizer for Householder { for k in 0..self.len() { self.reflect(k, a); } + if a.len() >= self.len() { + // full rank + return Zero::zero(); + } // residual norm a.slice(s![self.len()..]).norm_l2() } From bba765d34005c9712223a9edd2ea6cba2ae3276d Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 18:28:45 +0900 Subject: [PATCH 18/30] Bug fix --- src/krylov/householder.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index d05a260b..97fe7820 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -21,7 +21,7 @@ impl Householder { let w = self.v[k].slice(s![k..]); let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..])); for l in k..self.dim { - a[l] -= c * w[l]; + a[l] -= c * w[l - k]; } } } From 4de8c96b5ec165602666090a1211a311a9a5c6df Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 18:31:06 +0900 Subject: [PATCH 19/30] impl Householder::get_q --- src/krylov/householder.rs | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 97fe7820..4330d2e3 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -70,6 +70,13 @@ impl Orthogonalizer for Householder { } fn get_q(&self) -> Q { - unimplemented!() + assert!(self.len() > 0); + let mut a = Array::eye(self.len()); + for mut col in a.axis_iter_mut(Axis(0)) { + for l in 0..self.len() { + self.reflect(l, &mut col); + } + } + a } } From fadcda7c8600b2065509bc71980208197fb35ab6 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 5 May 2019 18:32:41 +0900 Subject: [PATCH 20/30] Add test for Householder::append --- src/krylov/householder.rs | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 4330d2e3..97bc949f 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -80,3 +80,26 @@ impl Orthogonalizer for Householder { a } } + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + + #[test] + fn householder_append() { + let mut householder = Householder::new(3); + let coef = householder.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0], 1e-9).unwrap(); + + let coef = householder.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + + assert!(householder.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); + + if let Err(coef) = householder.append(array![1.0, 2.0, 0.0], 1e-9) { + close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); + } + } + +} From 99839b1eca6f72280b10d6c50a3f465c9713ad0b Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 7 May 2019 03:18:33 +0900 Subject: [PATCH 21/30] Householder debug... --- src/krylov/householder.rs | 42 +++++++++++++++++++++++++++------------ src/krylov/mgs.rs | 8 ++++---- src/krylov/mod.rs | 12 ++++++++--- 3 files changed, 42 insertions(+), 20 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 97bc949f..8d72c0ba 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -4,7 +4,7 @@ use num_traits::Zero; /// Iterative orthogonalizer using Householder reflection #[derive(Debug, Clone)] -pub struct Householder { +pub struct Householder { dim: usize, v: Vec>, } @@ -18,10 +18,12 @@ impl Householder { fn reflect>(&self, k: usize, a: &mut ArrayBase) { assert!(k < self.v.len()); assert_eq!(a.len(), self.dim); + let w = self.v[k].slice(s![k..]); - let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..])); - for l in k..self.dim { - a[l] -= c * w[l - k]; + let mut a_slice = a.slice_mut(s![k..]); + let c = A::from(2.0).unwrap() * w.inner(&a_slice); + for l in 0..self.dim - k { + a_slice[l] -= c * w[l]; } } } @@ -29,10 +31,6 @@ impl Householder { impl Orthogonalizer for Householder { type Elem = A; - fn new(dim: usize) -> Self { - Self { dim, v: Vec::new() } - } - fn dim(&self) -> usize { self.dim } @@ -48,8 +46,7 @@ impl Orthogonalizer for Householder { for k in 0..self.len() { self.reflect(k, a); } - if a.len() >= self.len() { - // full rank + if self.is_full() { return Zero::zero(); } // residual norm @@ -60,12 +57,30 @@ impl Orthogonalizer for Householder { where S: DataMut, { - let residual = self.orthogonalize(&mut a); - let coef = a.slice(s![..self.len()]).into_owned(); - if residual < rtol { + assert_eq!(a.len(), self.dim); + let alpha = self.orthogonalize(&mut a); + + // Generate coefficient vector + let mut coef = Array::zeros(self.len() + 1); + for i in 0..self.len() { + coef[i] = a[i]; + } + coef[self.len()] = A::from_real(alpha); + + if alpha < rtol { return Err(coef); } + + // Add reflector + let k = self.len(); + let xi = a[k].re(); + a[k] = A::from(if xi > Zero::zero() { xi + alpha } else { xi - alpha }).unwrap(); + let norm = a.slice(s![k..]).norm_l2(); + dbg!(alpha); + dbg!(norm); + azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm)} ); self.v.push(a.into_owned()); + dbg!(&self.v); Ok(coef) } @@ -98,6 +113,7 @@ mod tests { assert!(householder.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); if let Err(coef) = householder.append(array![1.0, 2.0, 0.0], 1e-9) { + dbg!(&coef); close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); } } diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index abb8ba0a..506dfae5 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -11,6 +11,10 @@ pub struct MGS { } impl MGS { + pub fn new(dim: usize) -> Self { + Self { dim, q: Vec::new() } + } + fn ortho(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, @@ -33,10 +37,6 @@ impl MGS { impl Orthogonalizer for MGS { type Elem = A; - fn new(dim: usize) -> Self { - Self { dim, q: Vec::new() } - } - fn dim(&self) -> usize { self.dim } diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 2006e8d9..bd4dd602 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -15,12 +15,18 @@ pub type R = Array2; pub trait Orthogonalizer { type Elem: Scalar; - /// Create empty linear space - fn new(dim: usize) -> Self; - fn dim(&self) -> usize; fn len(&self) -> usize; + /// check if the basis spans entire space + fn is_full(&self) -> bool { + self.len() == self.dim() + } + + fn is_empty(&self) -> bool { + self.len() == 0 + } + /// Orthogonalize given vector /// /// Panic From 94725b2e1980d0191627c9586df63d8d75fdc8dc Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Thu, 9 May 2019 00:23:47 +0900 Subject: [PATCH 22/30] Fix test of householder_append --- src/krylov/householder.rs | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 8d72c0ba..30dcd150 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -74,13 +74,11 @@ impl Orthogonalizer for Householder { // Add reflector let k = self.len(); let xi = a[k].re(); - a[k] = A::from(if xi > Zero::zero() { xi + alpha } else { xi - alpha }).unwrap(); + a[k] = A::from(if xi >= Zero::zero() { xi + alpha } else { xi - alpha }).unwrap(); let norm = a.slice(s![k..]).norm_l2(); - dbg!(alpha); - dbg!(norm); - azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm)} ); + azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); + azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); self.v.push(a.into_owned()); - dbg!(&self.v); Ok(coef) } @@ -105,16 +103,17 @@ mod tests { fn householder_append() { let mut householder = Householder::new(3); let coef = householder.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); + dbg!(&coef); close_l2(&coef, &array![1.0], 1e-9).unwrap(); let coef = householder.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + dbg!(&coef); + close_l2(&coef, &array![-1.0, 1.0], 1e-9).unwrap(); assert!(householder.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); - if let Err(coef) = householder.append(array![1.0, 2.0, 0.0], 1e-9) { dbg!(&coef); - close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![-2.0, 1.0, 0.0], 1e-9).unwrap(); } } From 5dce6b6cf331ce66f01b9872d0ca648109b14c0a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 11 May 2019 20:20:28 +0900 Subject: [PATCH 23/30] Add document --- src/krylov/householder.rs | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 30dcd150..23226614 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -5,11 +5,17 @@ use num_traits::Zero; /// Iterative orthogonalizer using Householder reflection #[derive(Debug, Clone)] pub struct Householder { + /// Dimension of orthogonalizer dim: usize, + + /// Store Householder reflector. + /// + /// The coefficient is copied into another array, and this does not contain v: Vec>, } impl Householder { + /// Create a new orthogonalizer pub fn new(dim: usize) -> Self { Householder { dim, v: Vec::new() } } From 8b092f069bceb17e7c1feed1eaaf786fbcf19bd0 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 11 May 2019 20:33:47 +0900 Subject: [PATCH 24/30] Bug fix of Householder::get_q --- src/krylov/householder.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 23226614..d36cac51 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -23,7 +23,7 @@ impl Householder { /// Take a Reflection `P = I - 2ww^T` fn reflect>(&self, k: usize, a: &mut ArrayBase) { assert!(k < self.v.len()); - assert_eq!(a.len(), self.dim); + assert_eq!(a.len(), self.dim, "Input array size mismaches to the dimension"); let w = self.v[k].slice(s![k..]); let mut a_slice = a.slice_mut(s![k..]); @@ -90,8 +90,9 @@ impl Orthogonalizer for Householder { fn get_q(&self) -> Q { assert!(self.len() > 0); - let mut a = Array::eye(self.len()); - for mut col in a.axis_iter_mut(Axis(0)) { + let mut a = Array::zeros((self.dim(), self.len())); + for (i, mut col) in a.axis_iter_mut(Axis(1)).enumerate() { + col[i] = A::one(); for l in 0..self.len() { self.reflect(l, &mut col); } From 83d7af4542e06472a73b8af2b4d815e03632353a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 11 May 2019 20:45:25 +0900 Subject: [PATCH 25/30] Fix assertions --- src/krylov/householder.rs | 11 ++++------- src/krylov/mgs.rs | 6 +++--- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index d36cac51..0f2c2486 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -104,23 +104,20 @@ impl Orthogonalizer for Householder { #[cfg(test)] mod tests { use super::*; - use crate::assert::*; + use crate::assert::close_l2; #[test] fn householder_append() { let mut householder = Householder::new(3); let coef = householder.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - dbg!(&coef); - close_l2(&coef, &array![1.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0], 1e-9); let coef = householder.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - dbg!(&coef); - close_l2(&coef, &array![-1.0, 1.0], 1e-9).unwrap(); + close_l2(&coef, &array![-1.0, 1.0], 1e-9); assert!(householder.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); if let Err(coef) = householder.append(array![1.0, 2.0, 0.0], 1e-9) { - dbg!(&coef); - close_l2(&coef, &array![-2.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![-2.0, 1.0, 0.0], 1e-9); } } diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index 506dfae5..ae9aa1a6 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -86,15 +86,15 @@ mod tests { fn mgs_append() { let mut mgs = MGS::new(3); let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![1.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0], 1e-9); let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + close_l2(&coef, &array![1.0, 1.0], 1e-9); assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { - close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); + close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9); } } From 96348e529086c06ece157a458e773f5bbc072281 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 11 May 2019 20:57:30 +0900 Subject: [PATCH 26/30] Fix reflection order for Q-matrix --- src/krylov/householder.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 0f2c2486..b77437e1 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -93,7 +93,7 @@ impl Orthogonalizer for Householder { let mut a = Array::zeros((self.dim(), self.len())); for (i, mut col) in a.axis_iter_mut(Axis(1)).enumerate() { col[i] = A::one(); - for l in 0..self.len() { + for l in (0..self.len()).rev() { self.reflect(l, &mut col); } } From 8f2300f50d16a9f17ce3d7e28510cb6b4342870b Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 11 May 2019 21:21:09 +0900 Subject: [PATCH 27/30] Fix some bug --- src/krylov/householder.rs | 11 ++++++----- tests/krylov.rs | 16 ++++++++-------- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index b77437e1..28ade020 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -65,6 +65,9 @@ impl Orthogonalizer for Householder { { assert_eq!(a.len(), self.dim); let alpha = self.orthogonalize(&mut a); + let k = self.len(); + let xi = a[k].re(); + let alpha = if xi >= Zero::zero() { -alpha } else { alpha }; // Generate coefficient vector let mut coef = Array::zeros(self.len() + 1); @@ -73,14 +76,12 @@ impl Orthogonalizer for Householder { } coef[self.len()] = A::from_real(alpha); - if alpha < rtol { + if alpha.abs() < rtol { return Err(coef); } // Add reflector - let k = self.len(); - let xi = a[k].re(); - a[k] = A::from(if xi >= Zero::zero() { xi + alpha } else { xi - alpha }).unwrap(); + a[k] = A::from_real(xi - alpha); let norm = a.slice(s![k..]).norm_l2(); azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); @@ -110,7 +111,7 @@ mod tests { fn householder_append() { let mut householder = Householder::new(3); let coef = householder.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![1.0], 1e-9); + close_l2(&coef, &array![-1.0], 1e-9); let coef = householder.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); close_l2(&coef, &array![-1.0, 1.0], 1e-9); diff --git a/tests/krylov.rs b/tests/krylov.rs index 2f6f5f82..6a1703c5 100644 --- a/tests/krylov.rs +++ b/tests/krylov.rs @@ -74,9 +74,9 @@ fn householder_full() { const N: usize = 5; let a: Array2 = random((N, N)); let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); } test::(1e-5); @@ -91,9 +91,9 @@ fn householder_half() { const N: usize = 4; let a: Array2 = random((N, N / 2)); let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); } test::(1e-5); @@ -111,22 +111,22 @@ fn householder_over() { // Terminate let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR"); // Skip let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); let qc: Array2 = conjugate(&q); assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); // Full let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); - assert_close_l2!(&q.dot(&r), &a, rtol); let qc: Array2 = conjugate(&q); assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a, rtol); } test::(1e-5); From 3ec6fa4b07efc500c21e53e12660069026643a43 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 12 May 2019 00:36:09 +0900 Subject: [PATCH 28/30] Check coef size and array size --- src/krylov/householder.rs | 23 ++++++++++++----------- src/krylov/mod.rs | 24 +++++++++--------------- 2 files changed, 21 insertions(+), 26 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 28ade020..e1a06a1b 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -64,26 +64,27 @@ impl Orthogonalizer for Householder { S: DataMut, { assert_eq!(a.len(), self.dim); - let alpha = self.orthogonalize(&mut a); let k = self.len(); - let xi = a[k].re(); - let alpha = if xi >= Zero::zero() { -alpha } else { alpha }; - - // Generate coefficient vector - let mut coef = Array::zeros(self.len() + 1); - for i in 0..self.len() { + let alpha = self.orthogonalize(&mut a); + let mut coef = Array::zeros(k + 1); + for i in 0..k { coef[i] = a[i]; } - coef[self.len()] = A::from_real(alpha); - - if alpha.abs() < rtol { + if alpha < rtol { + // linearly dependent + coef[k] = A::from_real(alpha); return Err(coef); } // Add reflector + assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len() + let xi = a[k].re(); + let alpha = if xi >= Zero::zero() { -alpha } else { alpha }; + coef[k] = A::from_real(alpha); + a[k] = A::from_real(xi - alpha); let norm = a.slice(s![k..]).norm_l2(); - azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); + azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); // this can be omitted azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); self.v.push(a.into_owned()); Ok(coef) diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index bd4dd602..fb5b15a4 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -29,6 +29,10 @@ pub trait Orthogonalizer { /// Orthogonalize given vector /// + /// Returns + /// -------- + /// Residual norm + /// /// Panic /// ------- /// - if the size of the input array mismaches to the dimension @@ -39,22 +43,12 @@ pub trait Orthogonalizer { /// Add new vector if the residual is larger than relative tolerance /// - /// ```rust - /// # use ndarray::*; - /// # use ndarray_linalg::{*, krylov::*}; - /// let mut mgs = krylov::MGS::new(3); - /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); - /// - /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + /// Returns + /// -------- + /// Coefficients to the `i`-th Q-vector /// - /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend - /// - /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { - /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector - /// } - /// ``` + /// - The size of array must be `self.len() + 1` + /// - The last element is the residual norm of input vector /// /// Panic /// ------- From e3e8b8c116d18f305ceba66bc2ff22175b73fb30 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 12 May 2019 04:36:29 +0900 Subject: [PATCH 29/30] Split tests --- tests/householder.rs | 96 ++++++++++++++++++++++++++++++ tests/krylov.rs | 136 ------------------------------------------- tests/mgs.rs | 96 ++++++++++++++++++++++++++++++ 3 files changed, 192 insertions(+), 136 deletions(-) create mode 100644 tests/householder.rs delete mode 100644 tests/krylov.rs create mode 100644 tests/mgs.rs diff --git a/tests/householder.rs b/tests/householder.rs new file mode 100644 index 00000000..adc4d9b2 --- /dev/null +++ b/tests/householder.rs @@ -0,0 +1,96 @@ +use ndarray::*; +use ndarray_linalg::{krylov::*, *}; + +fn over(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR"); + + // Skip + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + + // Full + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a, rtol); +} + +#[test] +fn over_f32() { + over::(1e-5); +} +#[test] +fn over_f64() { + over::(1e-9); +} +#[test] +fn over_c32() { + over::(1e-5); +} +#[test] +fn over_c64() { + over::(1e-9); +} + +fn full(rtol: A::Real) { + const N: usize = 5; + let a: Array2 = random((N, N)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); +} + +#[test] +fn full_f32() { + full::(1e-5); +} +#[test] +fn full_f64() { + full::(1e-9); +} +#[test] +fn full_c32() { + full::(1e-5); +} +#[test] +fn full_c64() { + full::(1e-9); +} + +fn half(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N / 2)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); +} + +#[test] +fn half_f32() { + half::(1e-5); +} +#[test] +fn half_f64() { + half::(1e-9); +} +#[test] +fn half_c32() { + half::(1e-5); +} +#[test] +fn half_c64() { + half::(1e-9); +} diff --git a/tests/krylov.rs b/tests/krylov.rs deleted file mode 100644 index 6a1703c5..00000000 --- a/tests/krylov.rs +++ /dev/null @@ -1,136 +0,0 @@ -use ndarray::*; -use ndarray_linalg::{krylov::*, *}; - -#[test] -fn mgs_full() { - fn test(rtol: A::Real) { - const N: usize = 5; - let a: Array2 = random((N, N)); - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} - -#[test] -fn mgs_half() { - fn test(rtol: A::Real) { - const N: usize = 4; - let a: Array2 = random((N, N / 2)); - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - assert_close_l2!(&q.dot(&r), &a, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} - -#[test] -fn mgs_over() { - fn test(rtol: A::Real) { - const N: usize = 4; - let a: Array2 = random((N, N * 2)); - - // Terminate - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - - // Skip - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); - let a_sub = a.slice(s![.., 0..N]); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - - // Full - let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); - assert_close_l2!(&q.dot(&r), &a, rtol); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} - -#[test] -fn householder_full() { - fn test(rtol: A::Real) { - const N: usize = 5; - let a: Array2 = random((N, N)); - let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); - assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} - -#[test] -fn householder_half() { - fn test(rtol: A::Real) { - const N: usize = 4; - let a: Array2 = random((N, N / 2)); - let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I"); - assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} - -#[test] -fn householder_over() { - fn test(rtol: A::Real) { - const N: usize = 4; - let a: Array2 = random((N, N * 2)); - - // Terminate - let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); - let a_sub = a.slice(s![.., 0..N]); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); - assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR"); - - // Skip - let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); - let a_sub = a.slice(s![.., 0..N]); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - assert_close_l2!(&q.dot(&r), &a_sub, rtol); - - // Full - let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); - let qc: Array2 = conjugate(&q); - assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); - assert_close_l2!(&q.dot(&r), &a, rtol); - } - - test::(1e-5); - test::(1e-9); - test::(1e-5); - test::(1e-9); -} diff --git a/tests/mgs.rs b/tests/mgs.rs new file mode 100644 index 00000000..f80bbf5e --- /dev/null +++ b/tests/mgs.rs @@ -0,0 +1,96 @@ +use ndarray::*; +use ndarray_linalg::{krylov::*, *}; + +fn full(rtol: A::Real) { + const N: usize = 5; + let a: Array2 = random((N, N)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn full_f32() { + full::(1e-5); +} +#[test] +fn full_f64() { + full::(1e-9); +} +#[test] +fn full_c32() { + full::(1e-5); +} +#[test] +fn full_c64() { + full::(1e-9); +} + +fn half(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N / 2)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); +} + +#[test] +fn half_f32() { + half::(1e-5); +} +#[test] +fn half_f64() { + half::(1e-9); +} +#[test] +fn half_c32() { + half::(1e-5); +} +#[test] +fn half_c64() { + half::(1e-9); +} + +fn over(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Skip + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Full + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn over_f32() { + over::(1e-5); +} +#[test] +fn over_f64() { + over::(1e-9); +} +#[test] +fn over_c32() { + over::(1e-5); +} +#[test] +fn over_c64() { + over::(1e-9); +} From 65be8658cfc1a04b667b8eafcfbc2b3beb4638aa Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 12 May 2019 05:19:17 +0900 Subject: [PATCH 30/30] Fix reflector settings --- src/krylov/householder.rs | 35 ++++++++--------------------------- 1 file changed, 8 insertions(+), 27 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index e1a06a1b..9e7cef9d 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -78,11 +78,14 @@ impl Orthogonalizer for Householder { // Add reflector assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len() - let xi = a[k].re(); - let alpha = if xi >= Zero::zero() { -alpha } else { alpha }; - coef[k] = A::from_real(alpha); - - a[k] = A::from_real(xi - alpha); + let alpha = if a[k].abs() > Zero::zero() { + a[k].mul_real(alpha / a[k].abs()) + } else { + A::from_real(alpha) + }; + coef[k] = alpha; + + a[k] -= alpha; let norm = a.slice(s![k..]).norm_l2(); azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); // this can be omitted azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); @@ -102,25 +105,3 @@ impl Orthogonalizer for Householder { a } } - -#[cfg(test)] -mod tests { - use super::*; - use crate::assert::close_l2; - - #[test] - fn householder_append() { - let mut householder = Householder::new(3); - let coef = householder.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![-1.0], 1e-9); - - let coef = householder.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - close_l2(&coef, &array![-1.0, 1.0], 1e-9); - - assert!(householder.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); - if let Err(coef) = householder.append(array![1.0, 2.0, 0.0], 1e-9) { - close_l2(&coef, &array![-2.0, 1.0, 0.0], 1e-9); - } - } - -}