append_if

Author: termoshtt

src/arnoldi.rs

@@ -1,5 +1,6 @@
 use crate::{generate::*, inner::*, norm::Norm, types::*};
 use ndarray::*;
+use num_traits::Zero;
 
 #[derive(Debug, Clone)]
 pub struct MGS<A> {
@@ -8,6 +9,11 @@ pub struct MGS<A> {
     r: Vec<Array1<A>>,
 }
 
+pub enum Dependence<A: Scalar> {
+    Orthogonal(A::Real),
+    Coefficient(Array1<A>),
+}
+
 impl<A: Scalar + Lapack> MGS<A> {
     pub fn new(dim: usize) -> Self {
         Self {
@@ -25,7 +31,16 @@ impl<A: Scalar + Lapack> MGS<A> {
         self.q.len()
     }
 
+    /// Add new vector, return residual norm
     pub fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> A::Real
+    where
+        S: Data<Elem = A>,
+    {
+        self.append_if(a, A::Real::zero()).unwrap()
+    }
+
+    /// Add new vector if the residual is larger than relative tolerance
+    pub fn append_if<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Option<A::Real>
     where
         S: Data<Elem = A>,
     {
@@ -39,11 +54,14 @@ impl<A: Scalar + Lapack> MGS<A> {
             coef[i] = c;
         }
         let nrm = a.norm_l2();
+        if nrm < rtol {
+            return None;
+        }
         coef[self.len()] = A::from_real(nrm);
         self.r.push(coef);
         azip!(mut a in { *a = *a / A::from_real(nrm) });
         self.q.push(a);
-        nrm
+        Some(nrm)
     }
 
     /// Get orthogonal basis as Q matrix
@@ -78,7 +96,7 @@ mod tests {
         assert_eq!(mgs.len(), 0);
     }
 
-    fn test<A: Scalar + Lapack>(rtol: A::Real) {
+    fn test_append<A: Scalar + Lapack>(rtol: A::Real) {
         let mut mgs: MGS<A> = MGS::new(N);
         let a: Array2<A> = crate::generate::random((N, 3));
         dbg!(&a);
@@ -100,22 +118,43 @@ mod tests {
     }
 
     #[test]
-    fn test_f32() {
-        test::<f32>(1e-5);
+    fn append() {
+        test_append::<f32>(1e-5);
+        test_append::<c32>(1e-5);
+        test_append::<f64>(1e-9);
+        test_append::<c64>(1e-9);
     }
 
-    #[test]
-    fn test_c32() {
-        test::<c32>(1e-5);
-    }
+    fn test_append_if<A: Scalar + Lapack>(rtol: A::Real) {
+        let mut mgs: MGS<A> = MGS::new(N);
+        let a: Array2<A> = crate::generate::random((N, 8));
+        dbg!(&a);
+        for col in a.axis_iter(Axis(1)) {
+            match mgs.append_if(col, rtol) {
+                Some(res) => {
+                    dbg!(res);
+                }
+                None => break,
+            }
+        }
+        let q = mgs.get_q();
+        dbg!(&q);
+        let r = mgs.get_r();
+        dbg!(&r);
 
-    #[test]
-    fn test_f64() {
-        test::<f64>(1e-9);
+        dbg!(q.dot(&r));
+        close_l2(&q.dot(&r), &a.slice(s![.., 0..N]), rtol).unwrap();
+
+        let qt: Array2<_> = conjugate(&q);
+        dbg!(qt.dot(&q));
+        close_l2(&qt.dot(&q), &Array2::eye(N), rtol).unwrap();
     }
 
     #[test]
-    fn test_c64() {
-        test::<c64>(1e-9);
+    fn append_if() {
+        test_append_if::<f32>(1e-5);
+        test_append_if::<c32>(1e-5);
+        test_append_if::<f64>(1e-9);
+        test_append_if::<c64>(1e-9);
     }
 }