diff --git a/src/base/blas_uninit.rs b/src/base/blas_uninit.rs index 8f05dd193..5da5c2bc1 100644 --- a/src/base/blas_uninit.rs +++ b/src/base/blas_uninit.rs @@ -108,14 +108,26 @@ pub unsafe fn axcpy_uninit( let y = y.data.as_mut_slice_unchecked(); let x = x.data.as_slice_unchecked(); + let xlen = number_of_elements(x.len(), rstride2); + if !b.is_zero() { - array_axcpy(status, y, a, x, c, b, rstride1, rstride2, x.len()); + array_axcpy(status, y, a, x, c, b, rstride1, rstride2, xlen); } else { - array_axc(status, y, a, x, c, rstride1, rstride2, x.len()); + array_axc(status, y, a, x, c, rstride1, rstride2, xlen); } } } +fn number_of_elements(len: usize, stride: usize) -> usize { + if len == 0 { + 0 + } else if len % stride == 0 { + len / stride + } else { + len / stride + 1 + } +} + /// Computes `y = alpha * a * x + beta * y`, where `a` is a matrix, `x` a vector, and /// `alpha, beta` two scalars. /// @@ -331,3 +343,36 @@ pub unsafe fn gemm_uninit< } } } + +#[cfg(test)] +mod test { + use super::*; + #[test] + fn number_of_elements_test() { + assert_eq!(number_of_elements(0, 0), 0); + } +} + +#[cfg(test)] +#[cfg(feature = "proptest-support")] +mod blas_proptest { + use super::*; + use proptest::{prop_assert, proptest}; + + proptest! { + #[test] + fn number_of_elements_proptest(len in 0..1024usize, stride in 1..32usize) { + let predicted_num = number_of_elements(len, stride); + println!("{}", predicted_num); + // last element must be inside the length + if predicted_num == 0 { + prop_assert!(len == 0); + } else { + prop_assert!((predicted_num-1) * stride < len); + } + + // element after that must be outside the length + prop_assert!(predicted_num * stride >= len); + } + } +} diff --git a/tests/core/blas.rs b/tests/core/blas.rs index 2c3ca6df3..42d0ded0b 100644 --- a/tests/core/blas.rs +++ b/tests/core/blas.rs @@ -1,4 +1,4 @@ -use na::{Matrix2, Vector3, geometry::Quaternion}; +use na::{Matrix, Matrix2, U2, Vector3, geometry::Quaternion}; use num_traits::{One, Zero}; #[test] @@ -21,6 +21,14 @@ fn gemm_noncommutative() { assert_eq!(res, Matrix2::zero()); } +#[test] +fn array_axcpy_nonstandard_stride() { + let matrix_a = Matrix::<_, U2, U2, _>::from_slice_with_strides(&[0, 0, 0, 0], 2, 1); + let matrix_b = Matrix::<_, U2, U2, _>::from_slice_with_strides(&[0, 0, 0, 0], 2, 1); + let res = matrix_b * matrix_a; + assert_eq!(res, Matrix2::zero()); +} + #[cfg(feature = "proptest-support")] mod blas_proptest { use crate::proptest::{PROPTEST_F64, PROPTEST_MATRIX_DIM};