Introduce decompose method for Orthogonalizer

Author: termoshtt

src/krylov/householder.rs

@@ -92,12 +92,27 @@ impl<A: Scalar + Lapack> Householder<A> {
         }
     }
 
-    fn eval_residual<S>(&self, a: &ArrayBase<S, Ix1>) -> A::Real
+    /// Compose coefficients array using reflected vector
+    fn compose_coefficients<S>(&self, a: &ArrayBase<S, Ix1>) -> Coefficients<A>
     where
         S: Data<Elem = A>,
     {
-        let l = self.v.len();
-        a.slice(s![l..]).norm_l2()
+        let k = self.len();
+        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);
+        c
+    }
+
+    /// Construct the residual vector from reflected vector
+    fn construct_residual<S>(&self, a: &mut ArrayBase<S, Ix1>)
+    where
+        S: DataMut<Elem = A>,
+    {
+        let k = self.len();
+        azip!(mut a( a.slice_mut(s![..k])) in { *a = A::zero() });
+        self.backward_reflection(a);
     }
 }
 
@@ -112,18 +127,23 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.v.len()
     }
 
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    where
+        S: DataMut<Elem = A>,
+    {
+        self.forward_reflection(a);
+        let coef = self.compose_coefficients(a);
+        self.construct_residual(a);
+        coef
+    }
+
     fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: Data<Elem = A>,
     {
         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
+        self.compose_coefficients(&a)
     }
 
     fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>

src/krylov/mgs.rs

@@ -20,13 +20,20 @@ impl<A: Scalar + Lapack> MGS<A> {
             q: Vec::new(),
         }
     }
+}
+
+impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
+    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
-    ///
-    pub fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: DataMut<Elem = A>,
     {
@@ -42,26 +49,14 @@ impl<A: Scalar + Lapack> MGS<A> {
         coef[self.len()] = A::from_real(nrm);
         coef
     }
-}
-
-impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
-    type Elem = A;
-
-    fn dim(&self) -> usize {
-        self.dimension
-    }
-
-    fn len(&self) -> usize {
-        self.q.len()
-    }
 
     fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
     where
         A: Lapack,
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
-        self.orthogonalize(&mut a)
+        self.decompose(&mut a)
     }
 
     fn append<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
@@ -70,7 +65,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
-        let coef = self.orthogonalize(&mut a);
+        let coef = self.decompose(&mut a);
         let nrm = coef[coef.len() - 1].re();
         if nrm < rtol {
             // Linearly dependent

src/krylov/mod.rs

@@ -23,8 +23,19 @@ pub type Q<A> = Array2<A>;
 ///
 pub type R<A> = Array2<A>;
 
+/// Array type for coefficients to the current basis
+///
+/// - The length must be `self.len() + 1`
+/// - Last component is the residual norm
+///
+pub type Coefficients<A> = Array1<A>;
+
 /// Trait for creating orthogonal basis from iterator of arrays
 ///
+/// Panic
+/// -------
+/// - if the size of the input array mismatches to the dimension
+///
 /// Example
 /// -------
 ///
@@ -64,37 +75,31 @@ pub trait Orthogonalizer {
         self.len() == 0
     }
 
-    /// Calculate the coefficient to the given basis and residual norm
+    /// Decompose given vector into the span of current basis and
+    /// its tangent space
     ///
-    /// - The length of the returned array must be `self.len() + 1`
-    /// - Last component is the residual norm
+    /// - `a` becomes the tangent vector
+    /// - The Coefficients to the current basis is returned.
     ///
-    /// Panic
-    /// -------
-    /// - if the size of the input array mismatches to the dimension
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
+    where
+        S: DataMut<Elem = Self::Elem>;
+
+    /// Calculate the coefficient to the current basis basis
+    ///
+    /// - This will be faster than `decompose` because the construction of the residual vector may
+    ///   requires more Calculation
     ///
-    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<Self::Elem>
+    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
     where
         S: Data<Elem = Self::Elem>;
 
     /// Add new vector if the residual is larger than relative tolerance
-    ///
-    /// Returns
-    /// --------
-    /// Coefficients to the `i`-th Q-vector
-    ///
-    /// - The size of array must be `self.len() + 1`
-    /// - The last element is the residual norm of input vector
-    ///
-    /// Panic
-    /// -------
-    /// - if the size of the input array mismatches to the dimension
-    ///
     fn append<S>(
         &mut self,
         a: ArrayBase<S, Ix1>,
         rtol: <Self::Elem as Scalar>::Real,
-    ) -> Result<Array1<Self::Elem>, Array1<Self::Elem>>
+    ) -> Result<Coefficients<Self::Elem>, Coefficients<Self::Elem>>
     where
         S: DataMut<Elem = Self::Elem>;