Use matrixmultiply cgemm and improve tests

Refactor tests to use generics in a bit smarter way.

We test both .dot() and general_mat_mul. Tests now use relative accuracy
more explicitly (this works better with generics, instead of using
approx).

Author: bluss

Cargo.toml

@@ -40,7 +40,7 @@ approx-0_5 = { package = "approx", version = "0.5", optional = true , default-fe
 cblas-sys = { version = "0.1.4", optional = true, default-features = false }
 libc = { version = "0.2.82", optional = true }
 
-matrixmultiply = { version = "0.3.0", default-features = false}
+matrixmultiply = { version = "0.3.2", default-features = false, features=["cgemm"] }
 
 serde = { version = "1.0", optional = true, default-features = false, features = ["alloc"] }
 rawpointer = { version = "0.2" }

benches/gemv.rs renamed to benches/gemv_gemm.rs

@@ -9,8 +9,13 @@
 extern crate test;
 use test::Bencher;
 
+use num_complex::Complex;
+use num_traits::{Float, One, Zero};
+
 use ndarray::prelude::*;
 
+use ndarray::LinalgScalar;
+use ndarray::linalg::general_mat_mul;
 use ndarray::linalg::general_mat_vec_mul;
 
 #[bench]
@@ -45,3 +50,27 @@ fn gemv_64_32(bench: &mut Bencher) {
         general_mat_vec_mul(1.0, &a, &x, 1.0, &mut y);
     });
 }
+
+#[bench]
+fn cgemm_100(bench: &mut Bencher) {
+    cgemm_bench::<f32>(100, bench);
+}
+
+#[bench]
+fn zgemm_100(bench: &mut Bencher) {
+    cgemm_bench::<f64>(100, bench);
+}
+
+fn cgemm_bench<A>(size: usize, bench: &mut Bencher)
+where
+    A: LinalgScalar + Float,
+{
+    let (m, k, n) = (size, size, size);
+    let a = Array::<Complex<A>, _>::zeros((m, k));
+
+    let x = Array::zeros((k, n));
+    let mut y = Array::zeros((m, n));
+    bench.iter(|| {
+        general_mat_mul(Complex::one(), &a, &x, Complex::zero(), &mut y);
+    });
+}

src/linalg/impl_linalg.rs

@@ -18,6 +18,9 @@ use std::any::TypeId;
 use std::mem::MaybeUninit;
 use alloc::vec::Vec;
 
+use num_complex::Complex;
+use num_complex::{Complex32 as c32, Complex64 as c64};
+
 #[cfg(feature = "blas")]
 use libc::c_int;
 #[cfg(feature = "blas")]
@@ -30,9 +33,6 @@ use cblas_sys as blas_sys;
 #[cfg(feature = "blas")]
 use cblas_sys::{CblasNoTrans, CblasRowMajor, CblasTrans, CBLAS_LAYOUT};
 
-#[cfg(feature = "blas")]
-use num_complex::{Complex32 as c32, Complex64 as c64};
-
 /// len of vector before we use blas
 #[cfg(feature = "blas")]
 const DOT_BLAS_CUTOFF: usize = 32;
@@ -505,7 +505,7 @@ fn mat_mul_general<A>(
     let (rsc, csc) = (c.strides()[0], c.strides()[1]);
     if same_type::<A, f32>() {
         unsafe {
-            ::matrixmultiply::sgemm(
+            matrixmultiply::sgemm(
                 m,
                 k,
                 n,
@@ -524,7 +524,7 @@ fn mat_mul_general<A>(
         }
     } else if same_type::<A, f64>() {
         unsafe {
-            ::matrixmultiply::dgemm(
+            matrixmultiply::dgemm(
                 m,
                 k,
                 n,
@@ -541,6 +541,48 @@ fn mat_mul_general<A>(
                 csc,
             );
         }
+    } else if same_type::<A, c32>() {
+        unsafe {
+            matrixmultiply::cgemm(
+                matrixmultiply::CGemmOption::Standard,
+                matrixmultiply::CGemmOption::Standard,
+                m,
+                k,
+                n,
+                complex_array(cast_as(&alpha)),
+                ap as *const _,
+                lhs.strides()[0],
+                lhs.strides()[1],
+                bp as *const _,
+                rhs.strides()[0],
+                rhs.strides()[1],
+                complex_array(cast_as(&beta)),
+                cp as *mut _,
+                rsc,
+                csc,
+            );
+        }
+    } else if same_type::<A, c64>() {
+        unsafe {
+            matrixmultiply::zgemm(
+                matrixmultiply::CGemmOption::Standard,
+                matrixmultiply::CGemmOption::Standard,
+                m,
+                k,
+                n,
+                complex_array(cast_as(&alpha)),
+                ap as *const _,
+                lhs.strides()[0],
+                lhs.strides()[1],
+                bp as *const _,
+                rhs.strides()[0],
+                rhs.strides()[1],
+                complex_array(cast_as(&beta)),
+                cp as *mut _,
+                rsc,
+                csc,
+            );
+        }
     } else {
         // It's a no-op if `c` has zero length.
         if c.is_empty() {
@@ -768,10 +810,17 @@ fn same_type<A: 'static, B: 'static>() -> bool {
 //
 // **Panics** if `A` and `B` are not the same type
 fn cast_as<A: 'static + Copy, B: 'static + Copy>(a: &A) -> B {
-    assert!(same_type::<A, B>());
+    assert!(same_type::<A, B>(), "expect type {} and {} to match",
+            std::any::type_name::<A>(), std::any::type_name::<B>());
     unsafe { ::std::ptr::read(a as *const _ as *const B) }
 }
 
+/// Return the complex in the form of an array [re, im]
+#[inline]
+fn complex_array<A: 'static + Copy>(z: Complex<A>) -> [A; 2] {
+    [z.re, z.im]
+}
+
 #[cfg(feature = "blas")]
 fn blas_compat_1d<A, S>(a: &ArrayBase<S, Ix1>) -> bool
 where

xtest-numeric/Cargo.toml

@@ -3,6 +3,7 @@ name = "numeric-tests"
 version = "0.1.0"
 authors = ["bluss"]
 publish = false
+edition = "2018"
 
 [dependencies]
 approx = "0.4"
@@ -17,6 +18,10 @@ openblas-src = { optional = true, version = "0.10", default-features = false, fe
 version = "0.8.0"
 features = ["small_rng"]
 
+[dev-dependencies]
+num-traits = { version = "0.2.14", default-features = false }
+num-complex = { version = "0.4", default-features = false }
+
 [lib]
 test = false