Add some unit tests

Cargo.toml

@@ -21,3 +21,6 @@ statrs = "0.10"
 [dev-dependencies]
 approx = "0.3.1"
 openblas-src = "0.7"
+pyo3 = "0.8"
+numpy = "0.7.0"
+num-traits = "0.2"

src/lib.rs

@@ -240,7 +240,7 @@ impl PadeOrder for PadeOrder_13 {
               S2: DataMut<Elem=f64>,
               S3: DataMut<Elem=f64>,
     {
-        assert_eq!(a_powers.len(), (13 - 1)/2 + 1);
+        //assert_eq!(a_powers.len(), (13 - 1)/2 + 1);
 
         let (n_rows, n_cols) = a.dim();
         assert_eq!(n_rows, n_cols, "Pade sum only defined for square matrices.");
@@ -726,8 +726,11 @@ fn as_slice_with_layout_mut<S, T, D>(a: &mut ArrayBase<S, D>) -> Option<(&mut [T
 #[cfg(test)]
 mod tests {
     extern crate openblas_src;
+    extern crate ndarray;
     use ndarray::prelude::*;
-    use approx::assert_ulps_eq;
+    use approx::{assert_ulps_eq, relative_eq};
+    use std::f64::consts::E;
+    use ndarray::array;
 
     use crate::PadeOrder;
 
@@ -744,6 +747,46 @@ mod tests {
         factorial(2*m - i) * factorial(m) / factorial(2*m) / factorial(m-i) / factorial(i)
     }
 
+    use pyo3::prelude::{ObjectProtocol, Python};
+    use pyo3::types::PyDict;
+    use numpy::{PyArray, PyArray2};
+
+    macro_rules! python_testing {
+        ( $( $testname:ident, $scalar:ty, $arr:expr ),+ ) => {
+            $(
+
+            #[test]
+            fn $testname() {
+                let values = $arr;
+                let sqrt = (values.len() as f64).sqrt().round() as usize;
+                assert_eq!(sqrt * sqrt, values.len());
+                let n = sqrt;
+                let a = Array2::<$scalar>::from_shape_vec((n, n), values).unwrap();
+                let gil = Python::acquire_gil();
+                let py = gil.python();
+                let sla = py.import("scipy.linalg").expect("failed to import scipy");
+                let dict = PyDict::new(py);
+                let matrix = PyArray::from_array(py, &a);
+                dict.set_item("sla", sla).expect("failed to load sla");
+                dict.set_item("matrix", matrix).expect("failed to set matrix");
+                let pyarray: &PyArray2<$scalar> = py
+                    .eval("sla.expm(matrix)", Some(&dict), None)
+                    .expect("failed to execute exponentiation")
+                    .extract()
+                    .expect("faild to extract data form python");
+                let python_result: ndarray::Array2<$scalar> = pyarray.to_owned_array();
+
+                let mut b = Array2::zeros((n, n));
+                crate::expm(&a, &mut b);
+
+                for (b1, b2) in b.iter().zip(python_result.iter()){
+                    relative_eq!(b1, b2, epsilon=0.00000000001);
+                }
+            }
+            )+
+        }
+    }
+
     macro_rules! verify_pade_coefficients {
         ( $( $testname:ident, $padeorder:ty ),+ ) => {
             $(
@@ -779,16 +822,37 @@ mod tests {
         verify_pade_13, crate::PadeOrder_13
     );
 
+    python_testing!(
+        simple, f64, vec![1.0, 0.0, 1.0, 0.0],
+        double, f64, vec![2.0, 0.0, 2.0, 0.0],
+        random, f64, vec![1.02, -3.2, 4.2, 100.0]
+    );
+
     #[test]
     fn exp_of_unit() {
-        let n = 5;
-        let a = Array2::eye(n);
-        let mut b = Array2::<f64>::zeros((n, n));
+        for n in 2..10 {
+            let a = Array2::eye(n);
+            let mut b = Array2::<f64>::zeros((n, n));
+
+            crate::expm(&a, &mut b);
+
+            for &elem in &b.diag() {
+                assert_ulps_eq!(elem, E, max_ulps = 1);
+            }
+        }
+    }
 
-        crate::expm(&a, &mut b);
+    #[test]
+    fn exp_of_doubled_unit() {
+        for n in 2..10 {
+            let a = 2.0 * Array2::eye(n);
+            let mut b = Array2::<f64>::zeros((n, n));
 
-        for &elem in &b.diag() {
-            assert_ulps_eq!(elem, 1f64.exp(), max_ulps=1);
+            crate::expm(&a, &mut b);
+
+            for &elem in &b.diag() {
+                assert_ulps_eq!(elem, E * E, max_ulps = 30);
+            }
         }
     }
 }