Orthogonalizer::div_append

Author: termoshtt

src/krylov/householder.rs

@@ -8,17 +8,17 @@ use crate::{inner::*, norm::*};
 use num_traits::One;
 
 /// Calc a reflactor `w` from a vector `x`
-pub fn calc_reflector<A, S>(x: &mut ArrayBase<S, Ix1>) -> A
+pub fn calc_reflector<A, S>(x: &mut ArrayBase<S, Ix1>)
 where
     A: Scalar + Lapack,
     S: DataMut<Elem = A>,
 {
+    assert!(x.len() > 0);
     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)});
-    alpha
 }
 
 /// Take a reflection `P = I - 2ww^T`
@@ -108,7 +108,12 @@ impl<A: Scalar + Lapack> Householder<A> {
         let res = a.slice(s![k..]).norm_l2();
         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);
+        if k < a.len() {
+            let ak = a[k];
+            c[k] = -ak.mul_real(res / ak.abs());
+        } else {
+            c[k] = A::from_real(res);
+        }
         c
     }
 
@@ -157,30 +162,37 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.compose_coefficients(&a)
     }
 
-    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>) -> AppendResult<A>
+    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
     where
         S: DataMut<Elem = A>,
     {
         assert_eq!(a.len(), self.dim);
         let k = self.len();
-
-        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 AppendResult::Dependent(coef); // coef[k] must be zero in this case
+        self.forward_reflection(a);
+        let coef = self.compose_coefficients(a);
+        if coef[k].abs() < self.tol {
+            return AppendResult::Dependent(coef);
         }
+        calc_reflector(&mut a.slice_mut(s![k..]));
+        self.v.push(a.to_owned());
+        self.construct_residual(a);
+        AppendResult::Added(coef)
+    }
 
-        let alpha = calc_reflector(&mut a.slice_mut(s![k..]));
-        coef[k] = alpha;
-
-        if alpha.abs() < self.tol {
-            // linearly dependent
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
+    where
+        S: Data<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let mut a = a.into_owned();
+        let k = self.len();
+        self.forward_reflection(&mut a);
+        let coef = self.compose_coefficients(&a);
+        if coef[k].abs() < self.tol {
             return AppendResult::Dependent(coef);
         }
-        self.v.push(a.into_owned());
+        calc_reflector(&mut a.slice_mut(s![k..]));
+        self.v.push(a.to_owned());
         AppendResult::Added(coef)
     }
 

src/krylov/mgs.rs

@@ -74,14 +74,22 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
+        self.div_append(&mut a)
+    }
+
+    fn div_append<S>(&mut self, mut a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
         let coef = self.decompose(&mut a);
         let nrm = coef[coef.len() - 1].re();
         if nrm < self.tol {
             // Linearly dependent
             return AppendResult::Dependent(coef);
         }
-        azip!(mut a in { *a = *a / A::from_real(nrm) });
-        self.q.push(a);
+        azip!(mut a(&mut *a) in { *a = *a / A::from_real(nrm) });
+        self.q.push(a.to_owned());
         AppendResult::Added(coef)
     }
 

src/krylov/mod.rs

@@ -99,6 +99,12 @@ pub trait Orthogonalizer {
 
     /// Add new vector if the residual is larger than relative tolerance
     fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
+    where
+        S: Data<Elem = Self::Elem>;
+
+    /// Add new vector if the residual is larger than relative tolerance,
+    /// and return the residual vector
+    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
     where
         S: DataMut<Elem = Self::Elem>;