BUG: Fix blas + axpy by only using it for certain cases

Only use blas axpy for 1D arrays or contiguous arrays.

Author: bluss

src/linalg/impl_linalg.rs

@@ -281,8 +281,10 @@ impl<A, S, D> ArrayBase<S, D>
           D: Dimension,
 {
     /// Perform the operation `self += alpha * rhs` efficiently, where
-    /// `alpha` is a scalar and `rhs` is another array. This operation is
-    /// also known as `axpy` in BLAS.
+    /// `alpha` is a scalar and `rhs` is another array.
+    ///
+    /// This operation calls BLAS `axpy` if blas is enabled, the scalar type
+    /// is `f32` or `f64` and the dimensions are compatible.
     ///
     /// If their shapes disagree, `rhs` is broadcast to the shape of `self`.
     ///
@@ -322,25 +324,19 @@ impl<A, S, D> ArrayBase<S, D>
               A: LinalgScalar,
               E: Dimension,
     {
-        debug_assert_eq!(self.len(), rhs.len());
-        assert!(self.len() == rhs.len());
-        {
+        if (same_type::<A, f32>() || same_type::<A, f64>()) &&
+            self.len() == rhs.len() {
             macro_rules! axpy {
                 ($ty:ty, $func:ident) => {{
-            if blas_compat::<$ty, _, _>(self) && blas_compat::<$ty, _, _>(rhs) {
-                let order = Dimension::_fastest_varying_stride_order(&self.strides);
-                let incx = self.strides()[order[0]];
-
-                let order = Dimension::_fastest_varying_stride_order(&rhs.strides);
-                let incy = self.strides()[order[0]];
-
+        if let Some(self_stride) = blas_compat_axpy::<$ty, _, _>(self) {
+            if let Some(rhs_stride) = blas_compat_axpy::<$ty, _, _>(rhs) {
                 unsafe {
                     let (lhs_ptr, n, incx) = blas_1d_params(self.ptr,
                                                             self.len(),
-                                                            incx);
+                                                            self_stride);
                     let (rhs_ptr, _, incy) = blas_1d_params(rhs.ptr,
-                                                            rhs.len(),
-                                                            incy);
+                                                            self.len(),
+                                                            rhs_stride);
                     blas_sys::c::$func(
                         n,
                         cast_as(&alpha),
@@ -351,8 +347,10 @@ impl<A, S, D> ArrayBase<S, D>
                     return;
                 }
             }
-                }}
+        }
+            }}
             }
+
             axpy!{f32, cblas_saxpy};
             axpy!{f64, cblas_daxpy};
         }
@@ -595,29 +593,29 @@ fn blas_compat_1d<A, S>(a: &ArrayBase<S, Ix1>) -> bool
 }
 
 #[cfg(feature="blas")]
-fn blas_compat<A, S, D>(a: &ArrayBase<S, D>) -> bool
+// Return the interelement stride if the memory layout is blas compatible
+fn blas_compat_axpy<A, S, D>(a: &ArrayBase<S, D>) -> Option<isize>
     where S: Data,
           A: 'static,
           S::Elem: 'static,
           D: Dimension,
 {
     if !same_type::<A, S::Elem>() {
-        return false;
+        return None
     }
 
-    match D::equispaced_stride(&a.raw_dim(), &a.strides) {
-        Some(stride) => {
-            if a.len() as isize * stride > blas_index::max_value() as isize ||
-                stride < blas_index::min_value() as isize {
-                return false;
-            }
-        },
-        None => {
-            return false;
+    // FIXME: Support more memory layouts
+    if a.ndim() == 1 {
+        let stride = a.strides()[0];
+        if (a.len() * stride.abs() as usize) < blas_index::max_value() as usize {
+            return Some(stride)
+        }
+    } else if D::is_contiguous(&a.dim, &a.strides) {
+        if a.len() < blas_index::max_value() as usize {
+            return Some(1)
         }
     }
-
-    true
+    None
 }
 
 #[cfg(feature="blas")]