Add tests of msg

Author: termoshtt

src/arnoldi.rs

@@ -130,6 +130,7 @@ pub enum Strategy {
     Full,
 }
 
+/// Online QR decomposition of vectors using modified Gram-Schmit algorithm
 pub fn mgs<A, S>(
     iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
     dim: usize,
@@ -152,11 +153,14 @@ where
             },
         }
     }
+    let n = ortho.len();
     let m = coefs.len();
-    let mut r = Array2::zeros((m, m));
-    for i in 0..m {
-        for j in 0..=i {
-            r[(j, i)] = coefs[i][j];
+    let mut r = Array2::zeros((n, m).f());
+    for j in 0..m {
+        for i in 0..n {
+            if i < coefs[j].len() {
+                r[(i, j)] = coefs[j][i];
+            }
         }
     }
     (ortho.get_q(), r)

tests/arnoldi.rs

@@ -0,0 +1,83 @@
+use ndarray::*;
+use ndarray_linalg::{arnoldi::*, *};
+
+fn qr_full<A: Scalar + Lapack>() {
+    const N: usize = 5;
+    let rtol: A::Real = A::real(1e-9);
+
+    let a: Array2<A> = random((N, N));
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+}
+
+#[test]
+fn qr_full_real() {
+    qr_full::<f64>();
+}
+
+#[test]
+fn qr_full_complex() {
+    qr_full::<c64>();
+}
+
+fn qr<A: Scalar + Lapack>() {
+    const N: usize = 4;
+    let rtol: A::Real = A::real(1e-9);
+
+    let a: Array2<A> = random((N, N / 2));
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol);
+}
+
+#[test]
+fn qr_real() {
+    qr::<f64>();
+}
+
+#[test]
+fn qr_complex() {
+    qr::<c64>();
+}
+
+fn qr_over<A: Scalar + Lapack>() {
+    const N: usize = 4;
+    let rtol: A::Real = A::real(1e-9);
+
+    let a: Array2<A> = random((N, N * 2));
+
+    // Terminate
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    let a_sub = a.slice(s![.., 0..N]);
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+
+    // Skip
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip);
+    let a_sub = a.slice(s![.., 0..N]);
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+
+    // Full
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+}
+
+#[test]
+fn qr_over_real() {
+    qr_over::<f64>();
+}
+
+#[test]
+fn qr_over_complex() {
+    qr_over::<c64>();
+}