[lab] Move tests into respective files

russell_lab/src/matrix/mat_mat_mul.rs

@@ -101,7 +101,35 @@ pub fn mat_mat_mul(c: &mut Matrix, alpha: f64, a: &Matrix, b: &Matrix, beta: f64
 #[cfg(test)]
 mod tests {
     use super::{mat_mat_mul, Matrix};
-    use crate::mat_approx_eq;
+    use crate::{mat_approx_eq, mat_norm, Norm};
+
+    fn naive_mat_mat_mul(c: &mut Matrix, alpha: f64, a: &Matrix, b: &Matrix) {
+        let (m, n) = c.dims();
+        let k = a.ncol();
+        if a.nrow() != m || b.nrow() != k || b.ncol() != n {
+            panic!("matrices are incompatible");
+        }
+        if m == 0 || n == 0 || k == 0 {
+            return;
+        }
+        for i in 0..m {
+            for j in 0..n {
+                c.set(i, j, 0.0);
+                for p in 0..k {
+                    c.add(i, j, alpha * a.get(i, p) * b.get(p, j));
+                }
+            }
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "matrices are incompatible")]
+    fn naive_mat_mat_mul_capture_errors() {
+        let a = Matrix::new(1, 0);
+        let b = Matrix::new(0, 0);
+        let mut c = Matrix::new(0, 0);
+        naive_mat_mat_mul(&mut c, 1.0, &a, &b);
+    }
 
     #[test]
     fn mat_mat_mul_fails_on_wrong_dims() {
@@ -188,4 +216,28 @@ mod tests {
         ];
         mat_approx_eq(&c, correct, 1e-15);
     }
+
+    #[test]
+    fn mat_mat_mul_works_range() {
+        //   c  :=  a  ⋅  b
+        // (m,n)  (m,k) (k,n)
+        for m in [0, 5, 7_usize] {
+            for n in [0, 6, 12_usize] {
+                let mut c = Matrix::new(m, n);
+                let mut c_local = Matrix::new(m, n);
+                for k in [0, 5, 10, 15_usize] {
+                    let a = Matrix::filled(m, k, 1.0);
+                    let b = Matrix::filled(k, n, 1.0);
+                    mat_mat_mul(&mut c, 1.0, &a, &b, 0.0).unwrap();
+                    naive_mat_mat_mul(&mut c_local, 1.0, &a, &b);
+                    if m == 0 || n == 0 {
+                        assert_eq!(mat_norm(&c, Norm::Max), 0.0);
+                    } else {
+                        assert_eq!(mat_norm(&c, Norm::Max), k as f64);
+                    }
+                    mat_approx_eq(&c, &c_local, 1e-15);
+                }
+            }
+        }
+    }
 }

russell_lab/src/matvec/mat_vec_mul.rs

@@ -102,7 +102,7 @@ pub fn mat_vec_mul(v: &mut Vector, alpha: f64, a: &Matrix, u: &Vector) -> Result
 #[cfg(test)]
 mod tests {
     use super::{mat_vec_mul, Matrix, Vector};
-    use crate::vec_approx_eq;
+    use crate::{vec_approx_eq, vec_norm, Norm};
 
     #[test]
     fn mat_vec_mul_fails_on_wrong_dims() {
@@ -151,4 +151,23 @@ mod tests {
         mat_vec_mul(&mut v1, 1.0, &a_1x0, &u0).unwrap();
         assert_eq!(v1.as_data(), &[0.0]);
     }
+
+    #[test]
+    fn mat_vec_mul_works_range() {
+        // v  :=  a  ⋅ u
+        // (m)  (m,n) (n)
+        for m in [0, 7, 15_usize] {
+            for n in [0, 4, 8_usize] {
+                let a = Matrix::filled(m, n, 1.0);
+                let u = Vector::filled(n, 1.0);
+                let mut v = Vector::new(m);
+                mat_vec_mul(&mut v, 1.0, &a, &u).unwrap();
+                if m == 0 {
+                    assert_eq!(vec_norm(&v, Norm::Max), 0.0);
+                } else {
+                    assert_eq!(vec_norm(&v, Norm::Max), n as f64);
+                }
+            }
+        }
+    }
 }

russell_lab/src/matvec/solve_lin_sys.rs

@@ -113,7 +113,7 @@ pub fn solve_lin_sys(b: &mut Vector, a: &mut Matrix) -> Result<(), StrError> {
 #[cfg(test)]
 mod tests {
     use super::{solve_lin_sys, Matrix, Vector};
-    use crate::vec_approx_eq;
+    use crate::{mat_vec_mul, vec_add, vec_approx_eq, vec_norm, Norm};
 
     #[test]
     fn solve_lin_sys_fails_on_non_square() {
@@ -195,6 +195,45 @@ mod tests {
         vec_approx_eq(b.as_data(), x_correct, 1e-14);
     }
 
+    #[test]
+    fn solve_lin_sys_2_works() {
+        const TARGET: f64 = 1234.0;
+        for m in [0, 1, 5, 7, 12_usize] {
+            // prepare matrix and rhs
+            let mut a = Matrix::filled(m, m, 1.0);
+            let mut b = Vector::filled(m, TARGET);
+            for i in 0..m {
+                for j in (i + 1)..m {
+                    a.mul(i, j, -1.0);
+                }
+            }
+
+            // take copies
+            let a_copy = a.clone();
+            let b_copy = b.clone();
+
+            // solve linear system: b := a⁻¹⋅b == x
+            solve_lin_sys(&mut b, &mut a).unwrap();
+
+            // compare solution
+            if m == 0 {
+                let empty: [f64; 0] = [];
+                assert_eq!(b.as_data(), &empty);
+            } else {
+                let mut x_correct = Vector::new(m);
+                x_correct[0] = TARGET;
+                assert_eq!(b.as_data(), x_correct.as_data());
+            }
+
+            // check solution a ⋅ x == rhs (with x == b)
+            let mut rhs = Vector::new(m);
+            let mut diff = Vector::new(m);
+            mat_vec_mul(&mut rhs, 1.0, &a_copy, &b).unwrap();
+            vec_add(&mut diff, 1.0, &rhs, -1.0, &b_copy).unwrap();
+            assert_eq!(vec_norm(&diff, Norm::Max), 0.0);
+        }
+    }
+
     #[test]
     fn solve_lin_sys_singular_handles_error() {
         let mut a = Matrix::from(&[