Introduce AppendResult and move tolerance into orthogonalizer

Author: termoshtt

src/krylov/householder.rs

@@ -50,12 +50,19 @@ pub struct Householder<A: Scalar> {
     ///
     /// The coefficient is copied into another array, and this does not contain
     v: Vec<Array1<A>>,
+
+    /// Tolerance
+    tol: A::Real,
 }
 
 impl<A: Scalar + Lapack> Householder<A> {
     /// Create a new orthogonalizer
-    pub fn new(dim: usize) -> Self {
-        Householder { dim, v: Vec::new() }
+    pub fn new(dim: usize, tol: A::Real) -> Self {
+        Householder {
+            dim,
+            v: Vec::new(),
+            tol,
+        }
     }
 
     /// Take a Reflection `P = I - 2ww^T`
@@ -127,6 +134,10 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.v.len()
     }
 
+    fn tolerance(&self) -> A::Real {
+        self.tol
+    }
+
     fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: DataMut<Elem = A>,
@@ -146,7 +157,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.compose_coefficients(&a)
     }
 
-    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>) -> AppendResult<A>
     where
         S: DataMut<Elem = A>,
     {
@@ -159,18 +170,18 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
             coef[i] = a[i];
         }
         if self.is_full() {
-            return Err(coef); // coef[k] must be zero in this case
+            return AppendResult::Dependent(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 {
+        if alpha.abs() < self.tol {
             // linearly dependent
-            return Err(coef);
+            return AppendResult::Dependent(coef);
         }
         self.v.push(a.into_owned());
-        Ok(coef)
+        AppendResult::Added(coef)
     }
 
     fn get_q(&self) -> Q<A> {
@@ -195,8 +206,8 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    let h = Householder::new(dim);
-    qr(iter, h, rtol, strategy)
+    let h = Householder::new(dim, rtol);
+    qr(iter, h, strategy)
 }
 
 #[cfg(test)]

src/krylov/mgs.rs

@@ -5,19 +5,24 @@ use crate::{generate::*, inner::*, norm::Norm};
 
 /// Iterative orthogonalizer using modified Gram-Schmit procedure
 #[derive(Debug, Clone)]
-pub struct MGS<A> {
+pub struct MGS<A: Scalar> {
     /// Dimension of base space
-    dimension: usize,
+    dim: usize,
+
     /// Basis of spanned space
     q: Vec<Array1<A>>,
+
+    /// Tolerance
+    tol: A::Real,
 }
 
 impl<A: Scalar + Lapack> MGS<A> {
     /// Create an empty orthogonalizer
-    pub fn new(dimension: usize) -> Self {
+    pub fn new(dim: usize, tol: A::Real) -> Self {
         Self {
-            dimension,
+            dim,
             q: Vec::new(),
+            tol,
         }
     }
 }
@@ -26,13 +31,17 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
     type Elem = A;
 
     fn dim(&self) -> usize {
-        self.dimension
+        self.dim
     }
 
     fn len(&self) -> usize {
         self.q.len()
     }
 
+    fn tolerance(&self) -> A::Real {
+        self.tol
+    }
+
     fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: DataMut<Elem = A>,
@@ -59,21 +68,21 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         self.decompose(&mut a)
     }
 
-    fn append<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
     where
         A: Lapack,
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
         let coef = self.decompose(&mut a);
         let nrm = coef[coef.len() - 1].re();
-        if nrm < rtol {
+        if nrm < self.tol {
             // Linearly dependent
-            return Err(coef);
+            return AppendResult::Dependent(coef);
         }
         azip!(mut a in { *a = *a / A::from_real(nrm) });
         self.q.push(a);
-        Ok(coef)
+        AppendResult::Added(coef)
     }
 
     fn get_q(&self) -> Q<A> {
@@ -92,6 +101,6 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    let mgs = MGS::new(dim);
-    qr(iter, mgs, rtol, strategy)
+    let mgs = MGS::new(dim, rtol);
+    qr(iter, mgs, strategy)
 }

src/krylov/mod.rs

@@ -43,18 +43,18 @@ pub type Coefficients<A> = Array1<A>;
 /// ```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();
+/// let mut mgs = MGS::new(3, 1e-9);
+/// let coef = mgs.append(array![0.0, 1.0, 0.0]).into_coeff();
 /// close_l2(&coef, &array![1.0], 1e-9);
 ///
-/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap();
+/// let coef = mgs.append(array![1.0, 1.0, 0.0]).into_coeff();
 /// 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());
+/// assert!(mgs.append(array![1.0, 2.0, 0.0]).is_dependent());
 ///
 /// // You can get coefficients of dependent vector
-/// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) {
+/// if let AppendResult::Dependent(coef) = mgs.append(array![1.0, 2.0, 0.0]) {
 ///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9);
 /// }
 /// ```
@@ -76,6 +76,8 @@ pub trait Orthogonalizer {
         self.len() == 0
     }
 
+    fn tolerance(&self) -> <Self::Elem as Scalar>::Real;
+
     /// Decompose given vector into the span of current basis and
     /// its tangent space
     ///
@@ -96,18 +98,47 @@ pub trait Orthogonalizer {
         S: Data<Elem = Self::Elem>;
 
     /// Add new vector if the residual is larger than relative tolerance
-    fn append<S>(
-        &mut self,
-        a: ArrayBase<S, Ix1>,
-        rtol: <Self::Elem as Scalar>::Real,
-    ) -> Result<Coefficients<Self::Elem>, Coefficients<Self::Elem>>
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
     where
         S: DataMut<Elem = Self::Elem>;
 
     /// Get Q-matrix of generated basis
     fn get_q(&self) -> Q<Self::Elem>;
 }
 
+pub enum AppendResult<A> {
+    Added(Coefficients<A>),
+    Dependent(Coefficients<A>),
+}
+
+impl<A: Scalar> AppendResult<A> {
+    pub fn into_coeff(self) -> Coefficients<A> {
+        match self {
+            AppendResult::Added(c) => c,
+            AppendResult::Dependent(c) => c,
+        }
+    }
+
+    pub fn is_dependent(&self) -> bool {
+        match self {
+            AppendResult::Added(_) => false,
+            AppendResult::Dependent(_) => true,
+        }
+    }
+
+    pub fn coeff(&self) -> &Coefficients<A> {
+        match self {
+            AppendResult::Added(c) => c,
+            AppendResult::Dependent(c) => c,
+        }
+    }
+
+    pub fn residual_norm(&self) -> A::Real {
+        let c = self.coeff();
+        c[c.len() - 1].abs()
+    }
+}
+
 /// Strategy for linearly dependent vectors appearing in iterative QR decomposition
 #[derive(Clone, Copy, Debug, Eq, PartialEq)]
 pub enum Strategy {
@@ -133,7 +164,6 @@ pub enum Strategy {
 pub fn qr<A, S>(
     iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
     mut ortho: impl Orthogonalizer<Elem = A>,
-    rtol: A::Real,
     strategy: Strategy,
 ) -> (Q<A>, R<A>)
 where
@@ -144,9 +174,9 @@ where
 
     let mut coefs = Vec::new();
     for a in iter {
-        match ortho.append(a.into_owned(), rtol) {
-            Ok(coef) => coefs.push(coef),
-            Err(coef) => match strategy {
+        match ortho.append(a.into_owned()) {
+            AppendResult::Added(coef) => coefs.push(coef),
+            AppendResult::Dependent(coef) => match strategy {
                 Strategy::Terminate => break,
                 Strategy::Skip => continue,
                 Strategy::Full => coefs.push(coef),