Merge branch 'master' into ci_windows

Author: termoshtt

src/assert.rs

@@ -1,87 +1,87 @@
 //! Assertions for array
 
 use ndarray::*;
+use std::fmt::Debug;
 
 use super::norm::*;
 use super::types::*;
 
-/// check two values are close in terms of the relative torrence
-pub fn rclose<A: Scalar>(test: A, truth: A, rtol: A::Real) -> Result<A::Real, A::Real> {
+/// check two values are close in terms of the relative tolerance
+pub fn rclose<A: Scalar>(test: A, truth: A, rtol: A::Real) {
     let dev = (test - truth).abs() / truth.abs();
-    if dev < rtol {
-        Ok(dev)
-    } else {
-        Err(dev)
+    if dev > rtol {
+        eprintln!("==== Assetion Failed ====");
+        eprintln!("Expected = {}", truth);
+        eprintln!("Actual   = {}", test);
+        panic!("Too large deviation in relative tolerance: {}", dev);
     }
 }
 
-/// check two values are close in terms of the absolute torrence
-pub fn aclose<A: Scalar>(test: A, truth: A, atol: A::Real) -> Result<A::Real, A::Real> {
+/// check two values are close in terms of the absolute tolerance
+pub fn aclose<A: Scalar>(test: A, truth: A, atol: A::Real) {
     let dev = (test - truth).abs();
-    if dev < atol {
-        Ok(dev)
-    } else {
-        Err(dev)
+    if dev > atol {
+        eprintln!("==== Assetion Failed ====");
+        eprintln!("Expected = {}", truth);
+        eprintln!("Actual   = {}", test);
+        panic!("Too large deviation in absolute tolerance: {}", dev);
     }
 }
 
 /// check two arrays are close in maximum norm
-pub fn close_max<A, S1, S2, D>(
-    test: &ArrayBase<S1, D>,
-    truth: &ArrayBase<S2, D>,
-    atol: A::Real,
-) -> Result<A::Real, A::Real>
+pub fn close_max<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, atol: A::Real)
 where
     A: Scalar + Lapack,
     S1: Data<Elem = A>,
     S2: Data<Elem = A>,
     D: Dimension,
+    D::Pattern: PartialEq + Debug,
 {
+    assert_eq!(test.dim(), truth.dim());
     let tol = (test - truth).norm_max();
-    if tol < atol {
-        Ok(tol)
-    } else {
-        Err(tol)
+    if tol > atol {
+        eprintln!("==== Assetion Failed ====");
+        eprintln!("Expected:\n{}", truth);
+        eprintln!("Actual:\n{}", test);
+        panic!("Too large deviation in maximum norm: {} > {}", tol, atol);
     }
 }
 
 /// check two arrays are close in L1 norm
-pub fn close_l1<A, S1, S2, D>(
-    test: &ArrayBase<S1, D>,
-    truth: &ArrayBase<S2, D>,
-    rtol: A::Real,
-) -> Result<A::Real, A::Real>
+pub fn close_l1<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: A::Real)
 where
     A: Scalar + Lapack,
     S1: Data<Elem = A>,
     S2: Data<Elem = A>,
     D: Dimension,
+    D::Pattern: PartialEq + Debug,
 {
+    assert_eq!(test.dim(), truth.dim());
     let tol = (test - truth).norm_l1() / truth.norm_l1();
-    if tol < rtol {
-        Ok(tol)
-    } else {
-        Err(tol)
+    if tol > rtol {
+        eprintln!("==== Assetion Failed ====");
+        eprintln!("Expected:\n{}", truth);
+        eprintln!("Actual:\n{}", test);
+        panic!("Too large deviation in L1-norm: {} > {}", tol, rtol);
     }
 }
 
 /// check two arrays are close in L2 norm
-pub fn close_l2<A, S1, S2, D>(
-    test: &ArrayBase<S1, D>,
-    truth: &ArrayBase<S2, D>,
-    rtol: A::Real,
-) -> Result<A::Real, A::Real>
+pub fn close_l2<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: A::Real)
 where
     A: Scalar + Lapack,
     S1: Data<Elem = A>,
     S2: Data<Elem = A>,
     D: Dimension,
+    D::Pattern: PartialEq + Debug,
 {
+    assert_eq!(test.dim(), truth.dim());
     let tol = (test - truth).norm_l2() / truth.norm_l2();
-    if tol < rtol {
-        Ok(tol)
-    } else {
-        Err(tol)
+    if tol > rtol {
+        eprintln!("==== Assetion Failed ====");
+        eprintln!("Expected:\n{}", truth);
+        eprintln!("Actual:\n{}", test);
+        panic!("Too large deviation in L2-norm: {} > {} ", tol, rtol);
     }
 }
 
@@ -90,10 +90,11 @@ macro_rules! generate_assert {
         #[macro_export]
         macro_rules! $assert {
             ($test: expr,$truth: expr,$tol: expr) => {
-                $crate::$close($test, $truth, $tol).unwrap();
+                $crate::$close($test, $truth, $tol);
             };
             ($test: expr,$truth: expr,$tol: expr; $comment: expr) => {
-                $crate::$close($test, $truth, $tol).expect($comment);
+                eprintln!($comment);
+                $crate::$close($test, $truth, $tol);
             };
         }
     };

src/inner.rs

@@ -0,0 +1,30 @@
+use crate::types::*;
+use ndarray::*;
+
+/// Inner Product
+///
+/// Differenct from `Dot` trait, this take complex conjugate of `self` elements
+///
+pub trait InnerProduct {
+    type Elem: Scalar;
+
+    /// Inner product `(self.conjugate, rhs)
+    fn inner<S>(&self, rhs: &ArrayBase<S, Ix1>) -> Self::Elem
+    where
+        S: Data<Elem = Self::Elem>;
+}
+
+impl<A, S> InnerProduct for ArrayBase<S, Ix1>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Elem = A;
+    fn inner<St: Data<Elem = A>>(&self, rhs: &ArrayBase<St, Ix1>) -> A {
+        assert_eq!(self.len(), rhs.len());
+        Zip::from(self)
+            .and(rhs)
+            .fold_while(A::zero(), |acc, s, r| FoldWhile::Continue(acc + s.conj() * *r))
+            .into_inner()
+    }
+}

src/krylov/mgs.rs

@@ -0,0 +1,108 @@
+//! Modified Gram-Schmit orthogonalizer
+
+use super::*;
+use crate::{generate::*, inner::*, norm::Norm};
+
+/// Iterative orthogonalizer using modified Gram-Schmit procedure
+///
+/// ```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();
+/// close_l2(&coef, &array![1.0], 1e-9);
+///
+/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap();
+/// 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());
+///
+/// // You can get coefficients of dependent vector
+/// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) {
+///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9);
+/// }
+/// ```
+#[derive(Debug, Clone)]
+pub struct MGS<A> {
+    /// Dimension of base space
+    dimension: usize,
+    /// Basis of spanned space
+    q: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> MGS<A> {
+    /// Create an empty orthogonalizer
+    pub fn new(dimension: usize) -> Self {
+        Self {
+            dimension,
+            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()
+    }
+
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim());
+        let mut coef = Array1::zeros(self.len() + 1);
+        for i in 0..self.len() {
+            let q = &self.q[i];
+            let c = q.inner(&a);
+            azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } );
+            coef[i] = c;
+        }
+        let nrm = a.norm_l2();
+        coef[self.len()] = A::from_real(nrm);
+        coef
+    }
+
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        A: Lapack,
+        S: Data<Elem = A>,
+    {
+        let mut a = a.into_owned();
+        let coef = self.orthogonalize(&mut a);
+        let nrm = coef[coef.len() - 1].re();
+        if nrm < rtol {
+            // Linearly dependent
+            return Err(coef);
+        }
+        azip!(mut a in { *a = *a / A::from_real(nrm) });
+        self.q.push(a);
+        Ok(coef)
+    }
+
+    fn get_q(&self) -> Q<A> {
+        hstack(&self.q).unwrap()
+    }
+}
+
+/// Online QR decomposition of vectors using modified Gram-Schmit algorithm
+pub fn mgs<A, S>(
+    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
+    dim: usize,
+    rtol: A::Real,
+    strategy: Strategy,
+) -> (Q<A>, R<A>)
+where
+    A: Scalar + Lapack,
+    S: Data<Elem = A>,
+{
+    let mgs = MGS::new(dim);
+    qr(iter, mgs, rtol, strategy)
+}