Fix blas usage with vector of stride <= 0

Stride == 0 is unsuppored for vector increments;
Stride < 0 would be supported but the code needs to be adapted to pass
the right pointer for this case (lowest in memory pointer).

Co-authored-by: bluss <[email protected]>

Author: jturner314

src/linalg/impl_linalg.rs

@@ -8,6 +8,8 @@
 
 use crate::imp_prelude::*;
 use crate::numeric_util;
+#[cfg(feature = "blas")]
+use crate::dimension::offset_from_low_addr_ptr_to_logical_ptr;
 
 use crate::{LinalgScalar, Zip};
 
@@ -649,6 +651,12 @@ unsafe fn general_mat_vec_mul_impl<A, S1, S2>(
                             }
                         };
 
+                        // Low addr in memory pointers required for x, y
+                        let x_offset = offset_from_low_addr_ptr_to_logical_ptr(&x.dim, &x.strides);
+                        let x_ptr = x.ptr.as_ptr().sub(x_offset);
+                        let y_offset = offset_from_low_addr_ptr_to_logical_ptr(&y.dim, &y.strides);
+                        let y_ptr = y.ptr.as_ptr().sub(y_offset);
+
                         let x_stride = x.strides()[0] as blas_index;
                         let y_stride = y.strides()[0] as blas_index;
 
@@ -660,10 +668,10 @@ unsafe fn general_mat_vec_mul_impl<A, S1, S2>(
                             cast_as(&alpha),            // alpha
                             a.ptr.as_ptr() as *const _, // a
                             a_stride,                   // lda
-                            x.ptr.as_ptr() as *const _, // x
+                            x_ptr as *const _,          // x
                             x_stride,
-                            cast_as(&beta),           // beta
-                            y.ptr.as_ptr() as *mut _, // x
+                            cast_as(&beta),             // beta
+                            y_ptr as *mut _,            // y
                             y_stride,
                         );
                         return;
@@ -719,7 +727,10 @@ where
         return false;
     }
     let stride = a.strides()[0];
-    if stride > blas_index::max_value() as isize || stride < blas_index::min_value() as isize {
+    if stride == 0
+        || stride > blas_index::max_value() as isize
+        || stride < blas_index::min_value() as isize
+    {
         return false;
     }
     true

xtest-blas/tests/oper.rs

@@ -21,6 +21,25 @@ fn mat_vec_product_1d() {
     assert_eq!(a.t().dot(&b), ans);
 }
 
+#[test]
+fn mat_vec_product_1d_broadcast() {
+    let a = arr2(&[[1.], [2.], [3.]]);
+    let b = arr1(&[1.]);
+    let b = b.broadcast(3).unwrap();
+    let ans = arr1(&[6.]);
+    assert_eq!(a.t().dot(&b), ans);
+}
+
+#[test]
+fn mat_vec_product_1d_inverted_axis() {
+    let a = arr2(&[[1.], [2.], [3.]]);
+    let mut b = arr1(&[1., 2., 3.]);
+    b.invert_axis(Axis(0));
+
+    let ans = arr1(&[3. + 4. + 3.]);
+    assert_eq!(a.t().dot(&b), ans);
+}
+
 fn range_mat(m: Ix, n: Ix) -> Array2<f32> {
     Array::linspace(0., (m * n) as f32 - 1., m * n)
         .into_shape((m, n))