Use std::arch SIMD and runtime target feature detection

.travis.yml

@@ -1,21 +1,49 @@
 language: rust
 sudo: false
 
-# run builds for all the trains (and more)
-rust:
-  - 1.12.0
-  - stable
-  - beta
-  - nightly
+matrix:
+  include:
+    - rust: 1.28.0
+      env:
+        TARGET=x86_64-unknown-linux-gnu
+    - rust: stable
+      env:
+        TARGET=x86_64-unknown-linux-gnu
+    - rust: stable
+      env:
+        TARGET=i686-unknown-linux-gnu
+    - rust: beta
+      env:
+        TARGET=x86_64-unknown-linux-gnu
+    - rust: nightly
+      env:
+        TARGET=x86_64-unknown-linux-gnu
+    - rust: nightly
+      env:
+        TARGET=aarch64-unknown-linux-gnu
+        BUILD_ONLY=1
+env:
+  global:
+    - HOST=x86_64-unknown-linux-gnu
+
+addons:
+  apt:
+    packages:
+      # needed for i686-unknown-linux-gnu target
+      - gcc-multilib
+install:
+  # "rustup error: cannot re-add" without this conditional check
+- if [[ $HOST != $TARGET ]]; then rustup target add $TARGET; fi
 
 # the main build
 script:
   - |
-      cargo build &&
-      cargo test &&
-      cargo test --release &&
-      cargo doc &&
-      cargo bench
+      cargo build --target=$TARGET &&
+      ([ -n "$BUILD_ONLY" ] || (
+      cargo test --target=$TARGET &&
+      cargo test --release --target=$TARGET &&
+      cargo doc --target=$TARGET &&
+      cargo bench --target=$TARGET ))
 
 branches:
   only:

Cargo.toml

@@ -12,8 +12,6 @@ description = "General matrix multiplication of f32 and f64 matrices in Rust. Su
 
 keywords = ["matrix", "sgemm", "dgemm"]
 
-build = "build.rs"
-
 [lib]
 bench = false
 

LICENSE-MIT

@@ -1,4 +1,4 @@
-Copyright (c) 2015
+Copyright (c) 2016 - 2018 Ulrik Sverdrup "bluss"
 
 Permission is hereby granted, free of charge, to any
 person obtaining a copy of this software and associated

benches/benchmarks.rs

@@ -40,43 +40,44 @@ macro_rules! mat_mul {
     };
 }
 
-benchmark_main!(mat_mul_f32, mat_mul_f64, ref_mat_mul_f32);
+benchmark_main!(mat_mul_f32, mat_mul_f64);
 
 mat_mul!{mat_mul_f32, sgemm,
     (m004, 4, 4, 4)
-    (m005, 5, 5, 5)
     (m006, 6, 6, 6)
-    (m007, 7, 7, 7)
     (m008, 8, 8, 8)
-    (m009, 9, 9, 9)
     (m012, 12, 12, 12)
     (m016, 16, 16, 16)
     (m032, 32, 32, 32)
     (m064, 64, 64, 64)
     (m127, 127, 127, 127)
+    /*
     (m256, 256, 256, 256)
     (m512, 512, 512, 512)
     (mix16x4, 32, 4, 32)
     (mix32x2, 32, 2, 32)
     (mix97, 97, 97, 125)
     (mix128x10000x128, 128, 10000, 128)
+    */
 }
 
 mat_mul!{mat_mul_f64, dgemm,
     (m004, 4, 4, 4)
-    (m007, 7, 7, 7)
+    (m006, 6, 6, 6)
     (m008, 8, 8, 8)
     (m012, 12, 12, 12)
     (m016, 16, 16, 16)
     (m032, 32, 32, 32)
     (m064, 64, 64, 64)
     (m127, 127, 127, 127)
+    /*
     (m256, 256, 256, 256)
     (m512, 512, 512, 512)
     (mix16x4, 32, 4, 32)
     (mix32x2, 32, 2, 32)
     (mix97, 97, 97, 125)
     (mix128x10000x128, 128, 10000, 128)
+    */
 }
 
 use std::ops::{Add, Mul};

blas-bench/Cargo.toml

@@ -0,0 +1,31 @@
+[package]
+name = "blas-bench"
+version = "0.1.0"
+authors = ["bluss"]
+publish = false
+
+license = "MIT/Apache-2.0"
+
+repository = "https://github.com/bluss/matrixmultiply/"
+documentation = ""
+
+description = "Blas benchmarks for comparison with matrixmultiply"
+
+keywords = ["matrix", "sgemm", "dgemm"]
+
+[lib]
+bench = false
+
+[[bench]]
+name = "benchmarks"
+harness = false
+
+[dependencies]
+rawpointer = "0.1"
+matrixmultiply = { path = ".." }
+blas = { version = "0.20", default-features = false }
+blas-src = { version = "0.2.0", default-features = false }
+
+
+[dev-dependencies]
+bencher = "0.1.2"

blas-bench/README.md

@@ -0,0 +1,7 @@
+
+Run BLAS benchmarks to compare with matrixmultiply.
+
+These tests are set up to run vs a system-installed openblas (see the build.rs file),
+because building all of openblas just to benchmark versus it is tedious.
+So make sure openblas is installed, or other library that supports the cblas interface,
+and tweak the build.rs file to suit.

blas-bench/benches/benchmarks.rs

@@ -0,0 +1,81 @@
+extern crate blas_bench;
+extern crate matrixmultiply;
+pub use matrixmultiply::sgemm;
+pub use matrixmultiply::dgemm;
+
+#[macro_use]
+extern crate bencher;
+extern crate blas;
+
+use std::os::raw::c_int;
+
+
+#[allow(non_camel_case_types)]
+type blas_index = c_int; // blas index type
+
+
+// Compute GFlop/s
+// by flop / s = 2 M N K / time
+
+
+benchmark_main!(blas_mat_mul_f32, blas_mat_mul_f64);
+
+macro_rules! blas_mat_mul {
+    ($modname:ident, $gemm:ident, $(($name:ident, $m:expr, $n:expr, $k:expr))+) => {
+        mod $modname {
+            use bencher::{Bencher};
+            use super::blas_index;
+            $(
+            pub fn $name(bench: &mut Bencher)
+            {
+                let a = vec![0.; $m * $n]; 
+                let b = vec![0.; $n * $k];
+                let mut c = vec![0.; $m * $k];
+                bench.iter(|| {
+                    unsafe {
+
+                            blas::$gemm(
+                            b'N',
+                            b'N',
+                            $m as blas_index, // m, rows of Op(a)
+                            $n as blas_index, // n, cols of Op(b)
+                            $k as blas_index, // k, cols of Op(a)
+                            1.,
+                            &a,
+                            $n, // lda
+                            &b,
+                            $k, // ldb
+                            0.,         // beta
+                            &mut c,
+                            $k, // ldc
+                            );
+                    }
+                });
+            }
+            )+
+        }
+        benchmark_group!{ $modname, $($modname::$name),+ }
+    };
+}
+
+blas_mat_mul!{blas_mat_mul_f32, sgemm,
+    (m004, 4, 4, 4)
+    (m006, 6, 6, 6)
+    (m008, 8, 8, 8)
+    (m012, 12, 12, 12)
+    (m016, 16, 16, 16)
+    (m032, 32, 32, 32)
+    (m064, 64, 64, 64)
+   (m127, 127, 127, 127)
+}
+
+blas_mat_mul!{blas_mat_mul_f64, dgemm,
+    (m004, 4, 4, 4)
+    (m006, 6, 6, 6)
+    (m008, 8, 8, 8)
+    (m012, 12, 12, 12)
+    (m016, 16, 16, 16)
+    (m032, 32, 32, 32)
+    (m064, 64, 64, 64)
+   (m127, 127, 127, 127)
+}

blas-bench/build.rs

@@ -0,0 +1,12 @@
+
+///
+/// This build script emits the openblas linking directive if requested
+///
+
+fn main() {
+    // Always linking openblas
+    // Compiling blas just for testing is tedious -- install it on your system
+    // and run this.
+    println!("cargo:rerun-if-changed=build.rs");
+    println!("cargo:rustc-link-lib={}=openblas", "dylib");
+}

blas-bench/src/lib.rs

@@ -0,0 +1,7 @@
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn it_works() {
+        assert_eq!(2 + 2, 4);
+    }
+}

spare kernels/x86_sse_sgemm.rs

@@ -0,0 +1,84 @@
+
+// 4x4 sse sgemm
+macro_rules! mm_transpose4 {
+    ($c0:expr, $c1:expr, $c2:expr, $c3:expr) => {{
+        // This is _MM_TRANSPOSE4_PS except we take variables, not references
+        let tmp0 = _mm_unpacklo_ps($c0, $c1);
+        let tmp2 = _mm_unpacklo_ps($c2, $c3);
+        let tmp1 = _mm_unpackhi_ps($c0, $c1);
+        let tmp3 = _mm_unpackhi_ps($c2, $c3);
+
+        $c0 = _mm_movelh_ps(tmp0, tmp2);
+        $c1 = _mm_movehl_ps(tmp2, tmp0);
+        $c2 = _mm_movelh_ps(tmp1, tmp3);
+        $c3 = _mm_movehl_ps(tmp3, tmp1);
+    }}
+}
+
+#[inline(always)]
+#[cfg(any(target_arch="x86", target_arch="x86_64"))]
+unsafe fn kernel_x86_sse(k: usize, alpha: T, a: *const T, b: *const T,
+                         beta: T, c: *mut T, rsc: isize, csc: isize)
+{
+    let mut ab = [_mm_setzero_ps(); MR];
+
+    let mut bv;
+    let (mut a, mut b) = (a, b);
+
+    // Compute A B
+    for _ in 0..k {
+        bv = _mm_load_ps(b as _); // aligned due to GemmKernel::align_to
+
+        loop_m!(i, {
+            // Compute ab_i += [ai b_j+0, ai b_j+1, ai b_j+2, ai b_j+3]
+            let aiv = _mm_set1_ps(at(a, i));
+            ab[i] = _mm_add_ps(ab[i], _mm_mul_ps(aiv, bv));
+        });
+
+        a = a.add(MR);
+        b = b.add(NR);
+    }
+
+    // Compute α (A B)
+    let alphav = _mm_set1_ps(alpha);
+    loop_m!(i, ab[i] = _mm_mul_ps(alphav, ab[i]));
+
+    macro_rules! c {
+        ($i:expr, $j:expr) => (c.offset(rsc * $i as isize + csc * $j as isize));
+    }
+
+    // C ← α A B + β C
+    let mut c = [_mm_setzero_ps(); MR];
+    let betav = _mm_set1_ps(beta);
+    if beta != 0. {
+        // Read C
+        if csc == 1 {
+            loop_m!(i, c[i] = _mm_loadu_ps(c![i, 0]));
+        } else if rsc == 1 {
+            loop_m!(i, c[i] = _mm_loadu_ps(c![0, i]));
+            mm_transpose4!(c[0], c[1], c[2], c[3]);
+        } else {
+            loop_m!(i, c[i] = _mm_set_ps(*c![i, 3], *c![i, 2], *c![i, 1], *c![i, 0]));
+        }
+        // Compute β C
+        loop_m!(i, c[i] = _mm_mul_ps(c[i], betav));
+    }
+
+    // Compute (α A B) + (β C)
+    loop_m!(i, c[i] = _mm_add_ps(c[i], ab[i]));
+
+    // Store C back to memory
+    if csc == 1 {
+        loop_m!(i, _mm_storeu_ps(c![i, 0], c[i]));
+    } else if rsc == 1 {
+        mm_transpose4!(c[0], c[1], c[2], c[3]);
+        loop_m!(i, _mm_storeu_ps(c![0, i], c[i]));
+    } else {
+        // extract the nth value of a vector using _mm_cvtss_f32 (extract lowest)
+        // in combination with shuffle (move nth value to first position)
+        loop_m!(i, *c![i, 0] = _mm_cvtss_f32(c[i]));
+        loop_m!(i, *c![i, 1] = _mm_cvtss_f32(_mm_shuffle_ps(c[i], c[i], 1)));
+        loop_m!(i, *c![i, 2] = _mm_cvtss_f32(_mm_shuffle_ps(c[i], c[i], 2)));
+        loop_m!(i, *c![i, 3] = _mm_cvtss_f32(_mm_shuffle_ps(c[i], c[i], 3)));
+    }
+}