From 54ff3b1c4964af504a2a011dc6caab97440d4357 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 03:44:46 +0900 Subject: [PATCH 1/9] Import householder reflection from 'arnoldi' branch --- src/krylov/householder.rs | 107 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 107 insertions(+) create mode 100644 src/krylov/householder.rs diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs new file mode 100644 index 00000000..9e7cef9d --- /dev/null +++ b/src/krylov/householder.rs @@ -0,0 +1,107 @@ +use super::*; +use crate::{inner::*, norm::*}; +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() } + } + + /// 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, "Input array size mismaches to the dimension"); + + let w = self.v[k].slice(s![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]; + } + } +} + +impl Orthogonalizer for Householder { + type Elem = A; + + 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); + } + if self.is_full() { + return Zero::zero(); + } + // residual norm + a.slice(s![self.len()..]).norm_l2() + } + + fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> + where + S: DataMut, + { + assert_eq!(a.len(), self.dim); + let k = self.len(); + let alpha = self.orthogonalize(&mut a); + let mut coef = Array::zeros(k + 1); + for i in 0..k { + coef[i] = a[i]; + } + 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 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) }); + self.v.push(a.into_owned()); + Ok(coef) + } + + fn get_q(&self) -> Q { + assert!(self.len() > 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()).rev() { + self.reflect(l, &mut col); + } + } + a + } +} From 2c4ed3d290b8c526a07b9d91f5599f3812c96b37 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 04:59:27 +0900 Subject: [PATCH 2/9] Replace orthogonalize into coeff --- src/krylov/mgs.rs | 38 ++++++++++++++++++++++++++------------ src/krylov/mod.rs | 9 ++++++--- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index 87e57c9b..c586b77c 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -39,19 +39,12 @@ impl MGS { q: Vec::new(), } } -} - -impl Orthogonalizer for MGS { - type Elem = A; - - fn dim(&self) -> usize { - self.dimension - } - - fn len(&self) -> usize { - self.q.len() - } + /// Orthogonalize given vector against to the current basis + /// + /// - Returned array is coefficients and residual norm + /// - `a` will contain the residual vector + /// fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, @@ -69,6 +62,27 @@ impl Orthogonalizer for MGS { coef[self.len()] = A::from_real(nrm); coef } +} + +impl Orthogonalizer for MGS { + type Elem = A; + + fn dim(&self) -> usize { + self.dimension + } + + fn len(&self) -> usize { + self.q.len() + } + + fn coeff(&self, a: ArrayBase) -> Array1 + where + A: Lapack, + S: Data, + { + let mut a = a.into_owned(); + self.orthogonalize(&mut a) + } fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 677b1601..6c56fd25 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -40,15 +40,18 @@ pub trait Orthogonalizer { self.len() == 0 } - /// Orthogonalize given vector using current basis + /// Calculate the coefficient to the given basis and residual norm + /// + /// - The length of the returned array must be `self.len() + 1` + /// - Last component is the residual norm /// /// Panic /// ------- /// - if the size of the input array mismatches to the dimension /// - fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + fn coeff(&self, a: ArrayBase) -> Array1 where - S: DataMut; + S: Data; /// Add new vector if the residual is larger than relative tolerance /// From bcfee7b6e863b5abe20f375f144b7235b27cead4 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 05:40:35 +0900 Subject: [PATCH 3/9] Refactoring Householder --- src/krylov/householder.rs | 67 +++++++++++++++++++++++++++++---------- src/krylov/mod.rs | 2 ++ 2 files changed, 53 insertions(+), 16 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 9e7cef9d..f55a90cb 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -14,14 +14,14 @@ pub struct Householder { v: Vec>, } -impl Householder { +impl Householder { /// Create a new orthogonalizer 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) { + fn fundamental_reflection>(&self, k: usize, a: &mut ArrayBase) { assert!(k < self.v.len()); assert_eq!(a.len(), self.dim, "Input array size mismaches to the dimension"); @@ -32,6 +32,42 @@ impl Householder { a_slice[l] -= c * w[l]; } } + + /// Take forward reflection `P = P_l ... P_1` + pub fn forward_reflection(&self, a: &mut ArrayBase) + where + S: DataMut, + { + assert!(a.len() == self.dim); + let l = self.v.len(); + for k in 0..l { + self.fundamental_reflection(k, a); + } + } + + /// Take backward reflection `P = P_1 ... P_l` + pub fn backward_reflection(&self, a: &mut ArrayBase) + where + S: DataMut, + { + assert!(a.len() == self.dim); + let l = self.v.len(); + for k in (0..l).rev() { + self.fundamental_reflection(k, a); + } + } + + fn eval_residual(&self, a: &ArrayBase) -> A::Real + where + S: Data, + { + let l = self.v.len(); + if l == self.dim { + Zero::zero() + } else { + a.slice(s![l..]).norm_l2() + } + } } impl Orthogonalizer for Householder { @@ -45,18 +81,18 @@ impl Orthogonalizer for Householder { self.v.len() } - fn orthogonalize(&self, a: &mut ArrayBase) -> A::Real + fn coeff(&self, a: ArrayBase) -> Array1 where - S: DataMut, + S: Data, { - for k in 0..self.len() { - self.reflect(k, a); - } - if self.is_full() { - return Zero::zero(); - } - // residual norm - a.slice(s![self.len()..]).norm_l2() + let mut a = a.into_owned(); + self.forward_reflection(&mut a); + let res = self.eval_residual(&a); + let k = self.len(); + let mut c = Array1::zeros(k + 1); + azip!(mut c(c.slice_mut(s![..k])), a(a.slice(s![..k])) in { *c = a }); + c[k] = A::from_real(res); + c } fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> @@ -65,7 +101,8 @@ impl Orthogonalizer for Householder { { assert_eq!(a.len(), self.dim); let k = self.len(); - let alpha = self.orthogonalize(&mut a); + self.forward_reflection(&mut a); + let alpha = self.eval_residual(&a); let mut coef = Array::zeros(k + 1); for i in 0..k { coef[i] = a[i]; @@ -98,9 +135,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()).rev() { - self.reflect(l, &mut col); - } + self.backward_reflection(&mut col); } a } diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 6c56fd25..e13d2ecd 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -3,8 +3,10 @@ use crate::types::*; use ndarray::*; +mod householder; mod mgs; +pub use householder::*; pub use mgs::{mgs, MGS}; /// Q-matrix From e4af53765a5e80b1f3ae701174a00dfeae1c28c2 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 06:30:39 +0900 Subject: [PATCH 4/9] Add householder iterative decomposition --- src/krylov/householder.rs | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index f55a90cb..9c94a42f 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -21,7 +21,10 @@ impl Householder { } /// Take a Reflection `P = I - 2ww^T` - fn fundamental_reflection>(&self, k: usize, a: &mut ArrayBase) { + fn fundamental_reflection(&self, k: usize, a: &mut ArrayBase) + where + S: DataMut, + { assert!(k < self.v.len()); assert_eq!(a.len(), self.dim, "Input array size mismaches to the dimension"); @@ -140,3 +143,18 @@ impl Orthogonalizer for Householder { a } } + +/// Online QR decomposition of vectors +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) +} From b6f64dda625b8eaebba81ffe6f129d49ec07db21 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 06:31:22 +0900 Subject: [PATCH 5/9] Test for householder --- tests/householder.rs | 96 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 tests/householder.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); +} From fb671333a1f46ee299a688e5dee1455d24d20c87 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 09:39:33 +0900 Subject: [PATCH 6/9] Split primal functions --- src/krylov/householder.rs | 62 ++++++++++++++++++++++++++++++++------- 1 file changed, 52 insertions(+), 10 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 9c94a42f..47232056 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -1,6 +1,34 @@ use super::*; use crate::{inner::*, norm::*}; -use num_traits::Zero; +use num_traits::{One, Zero}; + +/// Calc a reflactor `w` from a vector `x` +pub fn calc_reflector(x: &mut ArrayBase) +where + A: Scalar + Lapack, + S: DataMut, +{ + let norm = x.norm_l2(); + let alpha = x[0].mul_real(norm / x[0].abs()); + x[0] -= alpha; + let inv_rev_norm = A::Real::one() / x.norm_l2(); + azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)}); +} + +/// Take a reflection using `w` +pub fn reflect(w: &ArrayBase, a: &mut ArrayBase) +where + A: Scalar + Lapack, + S1: Data, + S2: DataMut, +{ + assert_eq!(w.len(), a.len()); + let n = a.len(); + let c = A::from(2.0).unwrap() * w.inner(&a); + for l in 0..n { + a[l] -= c * w[l]; + } +} /// Iterative orthogonalizer using Householder reflection #[derive(Debug, Clone)] @@ -27,13 +55,7 @@ impl Householder { { assert!(k < self.v.len()); 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..]); - let c = A::from(2.0).unwrap() * w.inner(&a_slice); - for l in 0..self.dim - k { - a_slice[l] -= c * w[l]; - } + reflect(&self.v[k].slice(s![k..]), &mut a.slice_mut(s![k..])); } /// Take forward reflection `P = P_l ... P_1` @@ -110,14 +132,15 @@ impl Orthogonalizer for Householder { for i in 0..k { coef[i] = a[i]; } + coef[k] = A::from_real(alpha); 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() + + // Add reflector let alpha = if a[k].abs() > Zero::zero() { a[k].mul_real(alpha / a[k].abs()) } else { @@ -158,3 +181,22 @@ where let h = Householder::new(dim); qr(iter, h, rtol, strategy) } + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + + #[test] + fn check_reflector() { + let mut a = array![c64::new(1.0, 1.0), c64::new(1.0, 0.0), c64::new(0.0, 1.0)]; + let mut w = a.clone(); + calc_reflector(&mut w); + reflect(&w, &mut a); + close_l2( + &a, + &array![c64::new(2.0.sqrt(), 2.0.sqrt()), c64::zero(), c64::zero()], + 1e-9, + ); + } +} From 4c1275aee9486e4512dcd1416a50be6b5a38654c Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 10:06:25 +0900 Subject: [PATCH 7/9] Fix reflection bug --- src/krylov/householder.rs | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 47232056..2941a8a9 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -9,7 +9,7 @@ where S: DataMut, { let norm = x.norm_l2(); - let alpha = x[0].mul_real(norm / x[0].abs()); + let alpha = -x[0].mul_real(norm / x[0].abs()); x[0] -= alpha; let inv_rev_norm = A::Real::one() / x.norm_l2(); azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)}); @@ -87,11 +87,7 @@ impl Householder { S: Data, { let l = self.v.len(); - if l == self.dim { - Zero::zero() - } else { - a.slice(s![l..]).norm_l2() - } + a.slice(s![l..]).norm_l2() } } @@ -125,15 +121,24 @@ impl Orthogonalizer for Householder { S: DataMut, { assert_eq!(a.len(), self.dim); - let k = self.len(); + self.forward_reflection(&mut a); + + let k = self.len(); let alpha = self.eval_residual(&a); + let alpha = if k < a.len() && a[k].abs() > Zero::zero() { + -a[k].mul_real(alpha / a[k].abs()) + } else { + A::from_real(alpha) + }; + let mut coef = Array::zeros(k + 1); for i in 0..k { coef[i] = a[i]; } - coef[k] = A::from_real(alpha); - if alpha < rtol { + coef[k] = alpha; + + if alpha.abs() < rtol { // linearly dependent return Err(coef); } @@ -141,16 +146,9 @@ impl Orthogonalizer for Householder { assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len() // Add reflector - 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::zero() }); azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); self.v.push(a.into_owned()); Ok(coef) @@ -195,7 +193,7 @@ mod tests { reflect(&w, &mut a); close_l2( &a, - &array![c64::new(2.0.sqrt(), 2.0.sqrt()), c64::zero(), c64::zero()], + &array![-c64::new(2.0.sqrt(), 2.0.sqrt()), c64::zero(), c64::zero()], 1e-9, ); } From cc94837a5e9f73f4850a2fec3f5205f80241be72 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 12 Jun 2019 10:26:12 +0900 Subject: [PATCH 8/9] Refactoring --- src/krylov/householder.rs | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 2941a8a9..2b429ec6 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -1,9 +1,9 @@ use super::*; use crate::{inner::*, norm::*}; -use num_traits::{One, Zero}; +use num_traits::One; /// Calc a reflactor `w` from a vector `x` -pub fn calc_reflector(x: &mut ArrayBase) +pub fn calc_reflector(x: &mut ArrayBase) -> A where A: Scalar + Lapack, S: DataMut, @@ -13,6 +13,7 @@ where x[0] -= alpha; let inv_rev_norm = A::Real::one() / x.norm_l2(); azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)}); + alpha } /// Take a reflection using `w` @@ -121,35 +122,24 @@ impl Orthogonalizer for Householder { S: DataMut, { assert_eq!(a.len(), self.dim); - - self.forward_reflection(&mut a); - let k = self.len(); - let alpha = self.eval_residual(&a); - let alpha = if k < a.len() && a[k].abs() > Zero::zero() { - -a[k].mul_real(alpha / a[k].abs()) - } else { - A::from_real(alpha) - }; + self.forward_reflection(&mut a); let mut coef = Array::zeros(k + 1); for i in 0..k { coef[i] = a[i]; } + if self.is_full() { + return Err(coef); // coef[k] must be zero in this case + } + + let alpha = calc_reflector(&mut a.slice_mut(s![k..])); coef[k] = alpha; if alpha.abs() < rtol { // linearly dependent return Err(coef); } - - assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len() - - // Add reflector - a[k] -= alpha; - let norm = a.slice(s![k..]).norm_l2(); - azip!(mut a (a.slice_mut(s![..k])) in { *a = A::zero() }); - azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) }); self.v.push(a.into_owned()); Ok(coef) } @@ -184,6 +174,7 @@ where mod tests { use super::*; use crate::assert::*; + use num_traits::Zero; #[test] fn check_reflector() { From 840bf262dd6930d5755db049307bcb51c6361b6b Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Thu, 13 Jun 2019 03:07:15 +0900 Subject: [PATCH 9/9] Fix document and pub setting --- src/krylov/householder.rs | 13 +++++++++++-- src/krylov/mgs.rs | 26 +++----------------------- src/krylov/mod.rs | 28 +++++++++++++++++++++++++--- 3 files changed, 39 insertions(+), 28 deletions(-) diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs index 2b429ec6..391b18c4 100644 --- a/src/krylov/householder.rs +++ b/src/krylov/householder.rs @@ -1,3 +1,8 @@ +//! Householder reflection +//! +//! - [Householder transformation - Wikipedia](https://en.wikipedia.org/wiki/Householder_transformation) +//! + use super::*; use crate::{inner::*, norm::*}; use num_traits::One; @@ -16,7 +21,11 @@ where alpha } -/// Take a reflection using `w` +/// Take a reflection `P = I - 2ww^T` +/// +/// Panic +/// ------ +/// - if the size of `w` and `a` mismaches pub fn reflect(w: &ArrayBase, a: &mut ArrayBase) where A: Scalar + Lapack, @@ -155,7 +164,7 @@ impl Orthogonalizer for Householder { } } -/// Online QR decomposition of vectors +/// Online QR decomposition using Householder reflection pub fn householder( iter: impl Iterator>, dim: usize, diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index c586b77c..9caff111 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -4,25 +4,6 @@ use super::*; use crate::{generate::*, inner::*, norm::Norm}; /// Iterative orthogonalizer using modified Gram-Schmit procedure -/// -/// ```rust -/// # use ndarray::*; -/// # use ndarray_linalg::{krylov::*, *}; -/// 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); -/// -/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); -/// close_l2(&coef, &array![1.0, 1.0], 1e-9); -/// -/// // Fail if the vector is linearly dependent -/// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); -/// -/// // You can get coefficients of dependent vector -/// 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); -/// } -/// ``` #[derive(Debug, Clone)] pub struct MGS { /// Dimension of base space @@ -31,7 +12,7 @@ pub struct MGS { q: Vec>, } -impl MGS { +impl MGS { /// Create an empty orthogonalizer pub fn new(dimension: usize) -> Self { Self { @@ -45,9 +26,8 @@ impl MGS { /// - Returned array is coefficients and residual norm /// - `a` will contain the residual vector /// - fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where - A: Lapack, S: DataMut, { assert_eq!(a.len(), self.dim()); @@ -106,7 +86,7 @@ impl Orthogonalizer for MGS { } } -/// Online QR decomposition of vectors using modified Gram-Schmit algorithm +/// Online QR decomposition using modified Gram-Schmit algorithm pub fn mgs( iter: impl Iterator>, dim: usize, diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index e13d2ecd..b53df0c7 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -3,10 +3,10 @@ use crate::types::*; use ndarray::*; -mod householder; -mod mgs; +pub mod householder; +pub mod mgs; -pub use householder::*; +pub use householder::{householder, Householder}; pub use mgs::{mgs, MGS}; /// Q-matrix @@ -24,6 +24,28 @@ pub type Q = Array2; pub type R = Array2; /// Trait for creating orthogonal basis from iterator of arrays +/// +/// Example +/// ------- +/// +/// ```rust +/// # use ndarray::*; +/// # use ndarray_linalg::{krylov::*, *}; +/// 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); +/// +/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); +/// close_l2(&coef, &array![1.0, 1.0], 1e-9); +/// +/// // Fail if the vector is linearly dependent +/// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); +/// +/// // You can get coefficients of dependent vector +/// 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); +/// } +/// ``` pub trait Orthogonalizer { type Elem: Scalar;