Householder debug...

Author: termoshtt

src/krylov/householder.rs

@@ -4,7 +4,7 @@ use num_traits::Zero;
 
 /// Iterative orthogonalizer using Householder reflection
 #[derive(Debug, Clone)]
-pub struct Householder<A> {
+pub struct Householder<A: Scalar> {
     dim: usize,
     v: Vec<Array1<A>>,
 }
@@ -18,21 +18,19 @@ impl<A: Scalar> Householder<A> {
     fn reflect<S: DataMut<Elem = A>>(&self, k: usize, a: &mut ArrayBase<S, Ix1>) {
         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];
         }
     }
 }
 
 impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
     type Elem = A;
 
-    fn new(dim: usize) -> Self {
-        Self { dim, v: Vec::new() }
-    }
-
     fn dim(&self) -> usize {
         self.dim
     }
@@ -48,8 +46,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         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<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
     where
         S: DataMut<Elem = A>,
     {
-        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();
         }
     }

src/krylov/mgs.rs

@@ -11,6 +11,10 @@ pub struct MGS<A> {
 }
 
 impl<A: Scalar> MGS<A> {
+    pub fn new(dim: usize) -> Self {
+        Self { dim, q: Vec::new() }
+    }
+
     fn ortho<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
     where
         A: Lapack,
@@ -33,10 +37,6 @@ impl<A: Scalar> MGS<A> {
 impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
     type Elem = A;
 
-    fn new(dim: usize) -> Self {
-        Self { dim, q: Vec::new() }
-    }
-
     fn dim(&self) -> usize {
         self.dim
     }

src/krylov/mod.rs

@@ -15,12 +15,18 @@ pub type R<A> = Array2<A>;
 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