diff --git a/blis/src/interface.rs b/blis/src/interface.rs index 92af48cc..ca0c1819 100644 --- a/blis/src/interface.rs +++ b/blis/src/interface.rs @@ -4,16 +4,16 @@ pub mod types; use rlst_common::types::Scalar; /// Compute expected size of a data slice from stride and shape. -pub fn get_expected_data_size(stride: (usize, usize), shape: (usize, usize)) -> usize { - if shape.0 == 0 || shape.1 == 0 { +pub fn get_expected_data_size(stride: [usize; 2], shape: [usize; 2]) -> usize { + if shape[0] == 0 || shape[1] == 0 { return 0; } - 1 + (shape.0 - 1) * stride.0 + (shape.1 - 1) * stride.1 + 1 + (shape[0] - 1) * stride[0] + (shape[1] - 1) * stride[1] } /// Panic if expected data size is not identical to actual data size. -pub fn assert_data_size(data: &[T], stride: (usize, usize), shape: (usize, usize)) { +pub fn assert_data_size(data: &[T], stride: [usize; 2], shape: [usize; 2]) { let expected = get_expected_data_size(stride, shape); assert_eq!( diff --git a/blis/src/interface/gemm.rs b/blis/src/interface/gemm.rs index 2f27c754..b1323e2b 100644 --- a/blis/src/interface/gemm.rs +++ b/blis/src/interface/gemm.rs @@ -48,20 +48,20 @@ macro_rules! impl_gemm { csc: usize, ) { match transa { - TransMode::NoTrans => assert_data_size(a, (rsa, csa), (m, k)), - TransMode::ConjNoTrans => assert_data_size(a, (rsa, csa), (m, k)), - TransMode::Trans => assert_data_size(a, (rsa, csa), (k, m)), - TransMode::ConjTrans => assert_data_size(a, (rsa, csa), (k, m)), + TransMode::NoTrans => assert_data_size(a, [rsa, csa], [m, k]), + TransMode::ConjNoTrans => assert_data_size(a, [rsa, csa], [m, k]), + TransMode::Trans => assert_data_size(a, [rsa, csa], [k, m]), + TransMode::ConjTrans => assert_data_size(a, [rsa, csa], [k, m]), } match transb { - TransMode::NoTrans => assert_data_size(b, (rsb, csb), (k, n)), - TransMode::ConjNoTrans => assert_data_size(b, (rsb, csb), (k, n)), - TransMode::Trans => assert_data_size(b, (rsb, csb), (n, k)), - TransMode::ConjTrans => assert_data_size(b, (rsb, csb), (n, k)), + TransMode::NoTrans => assert_data_size(b, [rsb, csb], [k, n]), + TransMode::ConjNoTrans => assert_data_size(b, [rsb, csb], [k, n]), + TransMode::Trans => assert_data_size(b, [rsb, csb], [n, k]), + TransMode::ConjTrans => assert_data_size(b, [rsb, csb], [n, k]), } - assert_data_size(c, (rsc, csc), (m, n)); + assert_data_size(c, [rsc, csc], [m, n]); unsafe { raw::$blis_gemm( diff --git a/common/src/traits/accessors.rs b/common/src/traits/accessors.rs index 0fb766d1..1ea7f436 100644 --- a/common/src/traits/accessors.rs +++ b/common/src/traits/accessors.rs @@ -76,16 +76,16 @@ pub trait ChunkedAccess { /// Get raw access to the underlying data. pub trait RawAccess { - type T: Scalar; + type Item: Scalar; /// Get a slice of the whole data. - fn data(&self) -> &[Self::T]; + fn data(&self) -> &[Self::Item]; } /// Get mutable raw access to the underlying data. pub trait RawAccessMut: RawAccess { /// Get a mutable slice of the whole data. - fn data_mut(&mut self) -> &mut [Self::T]; + fn data_mut(&mut self) -> &mut [Self::Item]; } /// Check if `multi_index` not out of bounds with respect to `shape`. diff --git a/dense/src/array.rs b/dense/src/array.rs index 38ed57e0..661f0ce5 100644 --- a/dense/src/array.rs +++ b/dense/src/array.rs @@ -3,9 +3,13 @@ use crate::base_array::BaseArray; use crate::data_container::VectorContainer; use rlst_common::traits::ChunkedAccess; +use rlst_common::traits::NumberOfElements; use rlst_common::traits::RandomAccessByRef; use rlst_common::traits::RandomAccessMut; +use rlst_common::traits::RawAccess; +use rlst_common::traits::RawAccessMut; use rlst_common::traits::Shape; +use rlst_common::traits::Stride; use rlst_common::traits::UnsafeRandomAccessByRef; use rlst_common::traits::UnsafeRandomAccessByValue; use rlst_common::traits::UnsafeRandomAccessMut; @@ -13,6 +17,7 @@ use rlst_common::types::DataChunk; use rlst_common::types::Scalar; pub mod iterators; +pub mod multiply; pub mod operations; pub mod operators; pub mod random; @@ -174,3 +179,49 @@ pub(crate) fn empty_chunk( valid_entries, }) } + +impl< + Item: Scalar, + ArrayImpl: UnsafeRandomAccessByValue + Shape + RawAccess, + const NDIM: usize, + > RawAccess for Array +{ + type Item = Item; + + fn data(&self) -> &[Self::Item] { + self.0.data() + } +} + +impl< + Item: Scalar, + ArrayImpl: UnsafeRandomAccessByValue + Shape + RawAccessMut, + const NDIM: usize, + > RawAccessMut for Array +{ + fn data_mut(&mut self) -> &mut [Self::Item] { + self.0.data_mut() + } +} + +impl< + Item: Scalar, + ArrayImpl: UnsafeRandomAccessByValue + Shape, + const NDIM: usize, + > NumberOfElements for Array +{ + fn number_of_elements(&self) -> usize { + self.shape().iter().product() + } +} + +impl< + Item: Scalar, + ArrayImpl: UnsafeRandomAccessByValue + Shape + Stride, + const NDIM: usize, + > Stride for Array +{ + fn stride(&self) -> [usize; NDIM] { + self.0.stride() + } +} diff --git a/dense/src/array/multiply.rs b/dense/src/array/multiply.rs new file mode 100644 index 00000000..ffdc2b54 --- /dev/null +++ b/dense/src/array/multiply.rs @@ -0,0 +1 @@ +//! Multiplication of Arrays diff --git a/dense/src/base_array.rs b/dense/src/base_array.rs index 7b6bfd65..5f499cec 100644 --- a/dense/src/base_array.rs +++ b/dense/src/base_array.rs @@ -3,7 +3,10 @@ use crate::data_container::{DataContainer, DataContainerMut}; use crate::layout::{ check_multi_index_in_bounds, convert_1d_nd_from_shape, convert_nd_raw, stride_from_shape, }; -use rlst_common::traits::{ChunkedAccess, UnsafeRandomAccessByValue, UnsafeRandomAccessMut}; +use rlst_common::traits::{ + ChunkedAccess, RawAccess, RawAccessMut, Stride, UnsafeRandomAccessByValue, + UnsafeRandomAccessMut, +}; use rlst_common::{ traits::{Shape, UnsafeRandomAccessByRef}, types::Scalar, @@ -109,3 +112,29 @@ impl, const N: usize, const ND } } } + +impl, const NDIM: usize> RawAccess + for BaseArray +{ + type Item = Item; + + fn data(&self) -> &[Self::Item] { + self.data.data() + } +} + +impl, const NDIM: usize> RawAccessMut + for BaseArray +{ + fn data_mut(&mut self) -> &mut [Self::Item] { + self.data.data_mut() + } +} + +impl, const NDIM: usize> Stride + for BaseArray +{ + fn stride(&self) -> [usize; NDIM] { + self.stride + } +} diff --git a/dense/src/lib.rs b/dense/src/lib.rs index e98abc80..7fd0502e 100644 --- a/dense/src/lib.rs +++ b/dense/src/lib.rs @@ -46,7 +46,7 @@ pub mod layout; pub mod macros; pub mod traits; // pub mod matrix; -// pub mod matrix_multiply; +pub mod matrix_multiply; // pub mod matrix_ref; // pub mod matrix_view; // pub mod op_containers; diff --git a/dense/src/matrix_multiply.rs b/dense/src/matrix_multiply.rs index 39e5a9c5..1077b091 100644 --- a/dense/src/matrix_multiply.rs +++ b/dense/src/matrix_multiply.rs @@ -12,332 +12,534 @@ //! The [MultiplyAdd] trait is currently only implemented for dynamic matrices. //! Implementations for fixed size matrices will be added in the future. -use crate::data_container::{DataContainer, DataContainerMut, VectorContainer}; -use crate::matrix::{GenericBaseMatrix, Matrix}; -use crate::types::*; -use crate::{traits::*, MatrixD}; -use num; -use rlst_blis; -use rlst_blis::interface::{gemm::Gemm, types::TransMode}; -use rlst_common::traits::FillFrom; - -/// This trait provides a high-level interface for the multiplication of a matrix -/// with another matrix. The result is a new matrix, hence memory allocation takes place. -pub trait Dot { - type Output; - - /// Return the matrix product with a right-hand side. - fn dot(&self, rhs: &Rhs) -> Self::Output; -} - -/// This trait is an interface for the `dgemm` operation `mat_c = alpha * mat_a * mat_b + beta * mat_c`. -pub trait MultiplyAdd< - Item: Scalar, - Data1: DataContainer, - Data2: DataContainer, - Data3: DataContainerMut, - S1: SizeIdentifier, - S2: SizeIdentifier, - S3: SizeIdentifier, -> -{ - /// Perform the operation `mat_c = alpha * mat_a * mat_b + beta * mat_c`. - fn multiply_add( - alpha: Item, - mat_a: &GenericBaseMatrix, - mat_b: &GenericBaseMatrix, - beta: Item, - mat_c: &mut GenericBaseMatrix, +// use crate::data_container::{DataContainer, DataContainerMut, VectorContainer}; +// use crate::matrix::{GenericBaseMatrix, Matrix}; +// use crate::types::*; +// use crate::{traits::*, MatrixD}; +// use num; +// use rlst_blis; +// use rlst_blis::interface::{gemm::Gemm, types::TransMode}; +// use rlst_common::traits::FillFrom; + +// /// This trait provides a high-level interface for the multiplication of a matrix +// /// with another matrix. The result is a new matrix, hence memory allocation takes place. +// pub trait Dot { +// type Output; + +// /// Return the matrix product with a right-hand side. +// fn dot(&self, rhs: &Rhs) -> Self::Output; +// } + +// /// This trait is an interface for the `dgemm` operation `mat_c = alpha * mat_a * mat_b + beta * mat_c`. +// pub trait MultiplyAdd< +// Item: Scalar, +// Data1: DataContainer, +// Data2: DataContainer, +// Data3: DataContainerMut, +// S1: SizeIdentifier, +// S2: SizeIdentifier, +// S3: SizeIdentifier, +// > +// { +// /// Perform the operation `mat_c = alpha * mat_a * mat_b + beta * mat_c`. +// fn multiply_add( +// alpha: Item, +// mat_a: &GenericBaseMatrix, +// mat_b: &GenericBaseMatrix, +// beta: Item, +// mat_c: &mut GenericBaseMatrix, +// ); +// } + +// // Matrix x Matrix = Matrix +// impl< +// T: Scalar, +// MatImpl1: MatrixImplTrait, +// MatImpl2: MatrixImplTrait, +// S1: SizeIdentifier, +// S2: SizeIdentifier, +// > Dot> for Matrix +// where +// T: MultiplyAdd< +// T, +// VectorContainer, +// VectorContainer, +// VectorContainer, +// Dynamic, +// Dynamic, +// Dynamic, +// >, +// { +// type Output = MatrixD; + +// fn dot(&self, rhs: &Matrix) -> Self::Output { +// // We evaluate self and the other matrix and then perform the multiplication. +// let mut left = crate::rlst_dynamic_mat![T, self.shape()]; +// let mut right = crate::rlst_dynamic_mat![T, rhs.shape()]; + +// left.fill_from(self); +// right.fill_from(rhs); + +// let mut res = crate::rlst_dynamic_mat!(T, (self.shape().0, rhs.shape().1)); + +// T::multiply_add( +// num::cast::(1.0).unwrap(), +// &left, +// &right, +// num::cast::(0.0).unwrap(), +// &mut res, +// ); +// res +// } +// } + +// impl< +// T: Scalar, +// Data1: DataContainer, +// Data2: DataContainer, +// Data3: DataContainerMut, +// > MultiplyAdd for T +// where +// T: Gemm, +// { +// fn multiply_add( +// alpha: T, +// mat_a: &GenericBaseMatrix, +// mat_b: &GenericBaseMatrix, +// beta: T, +// mat_c: &mut GenericBaseMatrix, +// ) { +// let dim1 = mat_a.layout().dim(); +// let dim2 = mat_b.layout().dim(); +// let dim3 = mat_c.layout().dim(); + +// assert!( +// (dim1.1 == dim2.0) & (dim3.0 == dim1.0) & (dim3.1 == dim2.1), +// "Matrix multiply incompatible dimensions for C = A * B: A = {:#?}, B = {:#?}, C = {:#?}", +// dim1, +// dim2, +// dim3 +// ); + +// let stride_c = mat_c.stride(); + +// ::gemm( +// TransMode::NoTrans, +// TransMode::NoTrans, +// dim3.0, +// dim3.1, +// dim1.1, +// alpha, +// mat_a.data(), +// mat_a.stride().0, +// mat_a.stride().1, +// mat_b.data(), +// mat_b.stride().0, +// mat_b.stride().1, +// beta, +// mat_c.data_mut(), +// stride_c.0, +// stride_c.1, +// ); +// } +// } + +// #[cfg(test)] +// mod test { + +// use super::*; +// use approx::assert_ulps_eq; +// use rand_distr::StandardNormal; +// use rlst_common::assert_matrix_relative_eq; +// use rlst_common::tools::RandScalar; +// use rlst_common::traits::Eval; + +// use rand::prelude::*; + +// fn matmul_expect< +// Item: Scalar, +// Data1: DataContainer, +// Data2: DataContainer, +// Data3: DataContainerMut, +// S1: SizeIdentifier, +// S2: SizeIdentifier, +// S3: SizeIdentifier, +// >( +// alpha: Item, +// mat_a: &GenericBaseMatrix, +// mat_b: &GenericBaseMatrix, +// beta: Item, +// mat_c: &mut GenericBaseMatrix, +// ) { +// let m = mat_a.layout().dim().0; +// let k = mat_a.layout().dim().1; +// let n = mat_b.layout().dim().1; + +// for m_index in 0..m { +// for n_index in 0..n { +// *mat_c.get_mut(m_index, n_index).unwrap() *= beta; +// for k_index in 0..k { +// *mat_c.get_mut(m_index, n_index).unwrap() += alpha +// * mat_a.get_value(m_index, k_index).unwrap() +// * mat_b.get_value(k_index, n_index).unwrap(); +// } +// } +// } +// } + +// macro_rules! matmul_test { +// ($Scalar:ty, $fname:ident) => { +// #[test] +// fn $fname() { +// let mut mat_a = crate::rlst_dynamic_mat!($Scalar, (4, 6)); +// let mut mat_b = crate::rlst_dynamic_mat!($Scalar, (6, 5)); +// let mut mat_c_actual = crate::rlst_dynamic_mat!($Scalar, (4, 5)); +// let mut mat_c_expect = crate::rlst_dynamic_mat!($Scalar, (4, 5)); + +// let dist = StandardNormal; + +// let mut rng = rand::rngs::StdRng::seed_from_u64(0); + +// mat_a.fill_from_equally_distributed(&mut rng); +// mat_b.fill_from_equally_distributed(&mut rng); +// mat_c_actual.fill_from_equally_distributed(&mut rng); + +// for index in 0..mat_c_actual.layout().number_of_elements() { +// *mat_c_expect.get1d_mut(index).unwrap() = +// mat_c_actual.get1d_value(index).unwrap(); +// } + +// let alpha = <$Scalar>::random_scalar(&mut rng, &dist); +// let beta = <$Scalar>::random_scalar(&mut rng, &dist); + +// matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); +// <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); + +// for index in 0..mat_c_expect.layout().number_of_elements() { +// let val1 = mat_c_actual.get1d_value(index).unwrap(); +// let val2 = mat_c_expect.get1d_value(index).unwrap(); +// assert_ulps_eq!(&val1, &val2, max_ulps = 100); +// } +// } +// }; +// } + +// macro_rules! col_matvec_test { +// ($Scalar:ty, $fname:ident) => { +// #[test] +// fn $fname() { +// let mut mat_a = crate::rlst_dynamic_mat![$Scalar, (4, 6)]; +// let mut mat_b = crate::rlst_col_vec![$Scalar, 6]; +// let mut mat_c_actual = crate::rlst_col_vec![$Scalar, 4]; +// let mut mat_c_expect = crate::rlst_col_vec![$Scalar, 4]; + +// let dist = StandardNormal; + +// let mut rng = rand::rngs::StdRng::seed_from_u64(0); + +// mat_a.fill_from_equally_distributed(&mut rng); +// mat_b.fill_from_equally_distributed(&mut rng); +// mat_c_actual.fill_from_equally_distributed(&mut rng); + +// for index in 0..mat_c_actual.layout().number_of_elements() { +// *mat_c_expect.get1d_mut(index).unwrap() = +// mat_c_actual.get1d_value(index).unwrap(); +// } + +// let alpha = <$Scalar>::random_scalar(&mut rng, &dist); +// let beta = <$Scalar>::random_scalar(&mut rng, &dist); + +// matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); +// <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); + +// for index in 0..mat_c_expect.layout().number_of_elements() { +// let val1 = mat_c_actual.get1d_value(index).unwrap(); +// let val2 = mat_c_expect.get1d_value(index).unwrap(); +// assert_ulps_eq!(&val1, &val2, max_ulps = 100); +// } +// } +// }; +// } + +// macro_rules! row_matvec_test { +// ($Scalar:ty, $fname:ident) => { +// #[test] +// fn $fname() { +// let mut mat_a = crate::rlst_row_vec![$Scalar, 4]; +// let mut mat_b = crate::rlst_dynamic_mat![$Scalar, (4, 6)]; +// let mut mat_c_actual = crate::rlst_row_vec![$Scalar, 6]; +// let mut mat_c_expect = crate::rlst_row_vec![$Scalar, 6]; + +// let dist = StandardNormal; + +// let mut rng = rand::rngs::StdRng::seed_from_u64(0); + +// mat_a.fill_from_equally_distributed(&mut rng); +// mat_b.fill_from_equally_distributed(&mut rng); +// mat_c_actual.fill_from_equally_distributed(&mut rng); + +// for index in 0..mat_c_actual.layout().number_of_elements() { +// *mat_c_expect.get1d_mut(index).unwrap() = +// mat_c_actual.get1d_value(index).unwrap(); +// } + +// let alpha = <$Scalar>::random_scalar(&mut rng, &dist); +// let beta = <$Scalar>::random_scalar(&mut rng, &dist); + +// matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); +// <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); + +// for index in 0..mat_c_expect.layout().number_of_elements() { +// let val1 = mat_c_actual.get1d_value(index).unwrap(); +// let val2 = mat_c_expect.get1d_value(index).unwrap(); +// assert_ulps_eq!(&val1, &val2, max_ulps = 10); +// } +// } +// }; +// } + +// // matmul_test!(f32, test_matmul_f32); +// matmul_test!(f64, test_matmul_f64); +// matmul_test!(c32, test_matmul_c32); +// matmul_test!(c64, test_matmul_c64); + +// row_matvec_test!(f32, test_row_matvec_f32); +// row_matvec_test!(f64, test_row_matvec_f64); +// row_matvec_test!(c32, test_row_matvec_c32); +// row_matvec_test!(c64, test_row_matvec_c64); + +// col_matvec_test!(f64, test_col_matvec_f64); +// col_matvec_test!(f32, test_col_matvec_f32); +// col_matvec_test!(c32, test_col_matvec_c32); +// col_matvec_test!(c64, test_col_matvec_c64); + +// #[test] +// fn test_dot_matvec() { +// let mut mat1 = crate::rlst_dynamic_mat![f64, (2, 2)]; +// let mut mat2 = crate::rlst_dynamic_mat![f64, (2, 2)]; + +// mat1.fill_from_seed_equally_distributed(0); +// mat2.fill_from_seed_equally_distributed(1); + +// let mut mat3 = crate::rlst_dynamic_mat![f64, (2, 3)]; + +// mat3.fill_from_seed_equally_distributed(2); + +// let actual = (mat1.view() + mat2.view()).dot(&mat3); +// let mut expect = crate::rlst_dynamic_mat![f64, (2, 3)]; + +// let mat_sum = (mat1.view() + mat2.view()).eval(); + +// for row in 0..2 { +// for col in 0..3 { +// for k in 0..2 { +// expect[[row, col]] += mat_sum[[row, k]] * mat3[[k, col]]; +// } +// } +// } + +// assert_matrix_relative_eq!(expect, actual, 1E-14); +// } +// } + +use rlst_blis::interface::gemm::Gemm; +use rlst_blis::interface::types::TransMode; +use rlst_common::traits::*; +use rlst_common::{traits::RawAccess, types::Scalar}; + +pub fn matrix_multiply< + Item: Scalar + Gemm, + MatA: RawAccess + Shape<2> + Stride<2>, + MatB: RawAccess + Shape<2> + Stride<2>, + MatC: RawAccessMut + Shape<2> + Stride<2>, +>( + transa: TransMode, + transb: TransMode, + alpha: Item, + mat_a: &MatA, + mat_b: &MatB, + beta: Item, + mat_c: &mut MatC, +) { + let m = mat_c.shape()[0]; + let n = mat_c.shape()[1]; + + let a_shape = match transa { + TransMode::NoTrans => mat_a.shape(), + TransMode::ConjNoTrans => mat_a.shape(), + TransMode::Trans => [mat_a.shape()[1], mat_a.shape()[0]], + TransMode::ConjTrans => [mat_a.shape()[1], mat_a.shape()[0]], + }; + + let b_shape = match transb { + TransMode::NoTrans => mat_b.shape(), + TransMode::ConjNoTrans => mat_b.shape(), + TransMode::Trans => [mat_b.shape()[1], mat_b.shape()[0]], + TransMode::ConjTrans => [mat_b.shape()[1], mat_b.shape()[0]], + }; + + assert_eq!(m, a_shape[0], "Wrong dimension. {} != {}", m, a_shape[0]); + assert_eq!(n, b_shape[1], "Wrong dimension. {} != {}", n, b_shape[1]); + assert_eq!( + a_shape[1], b_shape[0], + "Wrong dimension. {} != {}", + a_shape[1], b_shape[0] ); -} - -// Matrix x Matrix = Matrix -impl< - T: Scalar, - MatImpl1: MatrixImplTrait, - MatImpl2: MatrixImplTrait, - S1: SizeIdentifier, - S2: SizeIdentifier, - > Dot> for Matrix -where - T: MultiplyAdd< - T, - VectorContainer, - VectorContainer, - VectorContainer, - Dynamic, - Dynamic, - Dynamic, - >, -{ - type Output = MatrixD; - - fn dot(&self, rhs: &Matrix) -> Self::Output { - // We evaluate self and the other matrix and then perform the multiplication. - let mut left = crate::rlst_dynamic_mat![T, self.shape()]; - let mut right = crate::rlst_dynamic_mat![T, rhs.shape()]; - - left.fill_from(self); - right.fill_from(rhs); - - let mut res = crate::rlst_dynamic_mat!(T, (self.shape().0, rhs.shape().1)); - - T::multiply_add( - num::cast::(1.0).unwrap(), - &left, - &right, - num::cast::(0.0).unwrap(), - &mut res, - ); - res - } -} -impl< - T: Scalar, - Data1: DataContainer, - Data2: DataContainer, - Data3: DataContainerMut, - > MultiplyAdd for T -where - T: Gemm, -{ - fn multiply_add( - alpha: T, - mat_a: &GenericBaseMatrix, - mat_b: &GenericBaseMatrix, - beta: T, - mat_c: &mut GenericBaseMatrix, - ) { - let dim1 = mat_a.layout().dim(); - let dim2 = mat_b.layout().dim(); - let dim3 = mat_c.layout().dim(); - - assert!( - (dim1.1 == dim2.0) & (dim3.0 == dim1.0) & (dim3.1 == dim2.1), - "Matrix multiply incompatible dimensions for C = A * B: A = {:#?}, B = {:#?}, C = {:#?}", - dim1, - dim2, - dim3 - ); - - let stride_c = mat_c.stride(); - - ::gemm( - TransMode::NoTrans, - TransMode::NoTrans, - dim3.0, - dim3.1, - dim1.1, - alpha, - mat_a.data(), - mat_a.stride().0, - mat_a.stride().1, - mat_b.data(), - mat_b.stride().0, - mat_b.stride().1, - beta, - mat_c.data_mut(), - stride_c.0, - stride_c.1, - ); - } + let [rsa, csa] = mat_a.stride(); + let [rsb, csb] = mat_b.stride(); + let [rsc, csc] = mat_c.stride(); + + ::gemm( + transa, + transb, + m, + n, + a_shape[1], + alpha, + mat_a.data(), + rsa, + csa, + mat_b.data(), + rsb, + csb, + beta, + mat_c.data_mut(), + rsc, + csc, + ) } #[cfg(test)] mod test { + use rlst_common::assert_array_relative_eq; + use rlst_common::types::{c32, c64}; + use super::*; - use approx::assert_ulps_eq; - use rand_distr::StandardNormal; - use rlst_common::assert_matrix_relative_eq; - use rlst_common::tools::RandScalar; - use rlst_common::traits::Eval; - - use rand::prelude::*; - - fn matmul_expect< - Item: Scalar, - Data1: DataContainer, - Data2: DataContainer, - Data3: DataContainerMut, - S1: SizeIdentifier, - S2: SizeIdentifier, - S3: SizeIdentifier, - >( - alpha: Item, - mat_a: &GenericBaseMatrix, - mat_b: &GenericBaseMatrix, - beta: Item, - mat_c: &mut GenericBaseMatrix, - ) { - let m = mat_a.layout().dim().0; - let k = mat_a.layout().dim().1; - let n = mat_b.layout().dim().1; - - for m_index in 0..m { - for n_index in 0..n { - *mat_c.get_mut(m_index, n_index).unwrap() *= beta; - for k_index in 0..k { - *mat_c.get_mut(m_index, n_index).unwrap() += alpha - * mat_a.get_value(m_index, k_index).unwrap() - * mat_b.get_value(k_index, n_index).unwrap(); + use crate::array::Array; + use crate::base_array::BaseArray; + use crate::data_container::VectorContainer; + use crate::rlst_dynamic_array2; + use paste::paste; + + macro_rules! mat_mul_test_impl { + ($ScalarType:ty, $eps:expr) => { + paste! { + fn [](transa: TransMode, transb: TransMode, shape_a: [usize; 2], shape_b: [usize; 2], shape_c: [usize; 2]) { + + let mut mat_a = rlst_dynamic_array2!($ScalarType, shape_a); + let mut mat_b = rlst_dynamic_array2!($ScalarType, shape_b); + let mut mat_c = rlst_dynamic_array2!($ScalarType, shape_c); + let mut expected = rlst_dynamic_array2!($ScalarType, shape_c); + + mat_a.fill_from_seed_equally_distributed(0); + mat_b.fill_from_seed_equally_distributed(1); + //mat_c.fill_from_seed_equally_distributed(2); + + expected.fill_from(mat_c.view_mut()); + + matrix_multiply( + transa, + transb, + <$ScalarType as Scalar>::from_real(1.), + &mat_a, + &mat_b, + <$ScalarType as Scalar>::from_real(1.), + &mut mat_c, + ); + matrix_multiply_compare( + transa, + transb, + <$ScalarType as Scalar>::from_real(1.), + &mat_a, + &mat_b, + <$ScalarType as Scalar>::from_real(1.), + &mut expected, + ); + + assert_array_relative_eq!(mat_c, expected, $eps); } - } - } - } - macro_rules! matmul_test { - ($Scalar:ty, $fname:ident) => { - #[test] - fn $fname() { - let mut mat_a = crate::rlst_dynamic_mat!($Scalar, (4, 6)); - let mut mat_b = crate::rlst_dynamic_mat!($Scalar, (6, 5)); - let mut mat_c_actual = crate::rlst_dynamic_mat!($Scalar, (4, 5)); - let mut mat_c_expect = crate::rlst_dynamic_mat!($Scalar, (4, 5)); + #[test] + fn []() { - let dist = StandardNormal; + [](TransMode::NoTrans, TransMode::NoTrans, [3, 5], [5, 7], [3, 7]); + [](TransMode::ConjNoTrans, TransMode::ConjNoTrans, [3, 5], [5, 7], [3, 7]); + [](TransMode::ConjNoTrans, TransMode::NoTrans, [3, 5], [5, 7], [3, 7]); + [](TransMode::NoTrans, TransMode::ConjNoTrans, [3, 5], [5, 7], [3, 7]); - let mut rng = rand::rngs::StdRng::seed_from_u64(0); + [](TransMode::NoTrans, TransMode::Trans, [3, 5], [7, 5], [3, 7]); + [](TransMode::Trans, TransMode::NoTrans, [2, 1], [2, 1], [1, 1]); + [](TransMode::Trans, TransMode::Trans, [5, 3], [7, 5], [3, 7]); - mat_a.fill_from_equally_distributed(&mut rng); - mat_b.fill_from_equally_distributed(&mut rng); - mat_c_actual.fill_from_equally_distributed(&mut rng); + [](TransMode::NoTrans, TransMode::ConjTrans, [3, 5], [7, 5], [3, 7]); + [](TransMode::ConjTrans, TransMode::NoTrans, [5, 3], [5, 7], [3, 7]); + [](TransMode::ConjTrans, TransMode::ConjTrans, [5, 3], [7, 5], [3, 7]); - for index in 0..mat_c_actual.layout().number_of_elements() { - *mat_c_expect.get1d_mut(index).unwrap() = - mat_c_actual.get1d_value(index).unwrap(); - } - let alpha = <$Scalar>::random_scalar(&mut rng, &dist); - let beta = <$Scalar>::random_scalar(&mut rng, &dist); - - matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); - <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); - - for index in 0..mat_c_expect.layout().number_of_elements() { - let val1 = mat_c_actual.get1d_value(index).unwrap(); - let val2 = mat_c_expect.get1d_value(index).unwrap(); - assert_ulps_eq!(&val1, &val2, max_ulps = 100); } - } - }; - } - - macro_rules! col_matvec_test { - ($Scalar:ty, $fname:ident) => { - #[test] - fn $fname() { - let mut mat_a = crate::rlst_dynamic_mat![$Scalar, (4, 6)]; - let mut mat_b = crate::rlst_col_vec![$Scalar, 6]; - let mut mat_c_actual = crate::rlst_col_vec![$Scalar, 4]; - let mut mat_c_expect = crate::rlst_col_vec![$Scalar, 4]; - - let dist = StandardNormal; - let mut rng = rand::rngs::StdRng::seed_from_u64(0); - - mat_a.fill_from_equally_distributed(&mut rng); - mat_b.fill_from_equally_distributed(&mut rng); - mat_c_actual.fill_from_equally_distributed(&mut rng); - - for index in 0..mat_c_actual.layout().number_of_elements() { - *mat_c_expect.get1d_mut(index).unwrap() = - mat_c_actual.get1d_value(index).unwrap(); - } - - let alpha = <$Scalar>::random_scalar(&mut rng, &dist); - let beta = <$Scalar>::random_scalar(&mut rng, &dist); - - matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); - <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); - - for index in 0..mat_c_expect.layout().number_of_elements() { - let val1 = mat_c_actual.get1d_value(index).unwrap(); - let val2 = mat_c_expect.get1d_value(index).unwrap(); - assert_ulps_eq!(&val1, &val2, max_ulps = 100); - } } }; } - macro_rules! row_matvec_test { - ($Scalar:ty, $fname:ident) => { - #[test] - fn $fname() { - let mut mat_a = crate::rlst_row_vec![$Scalar, 4]; - let mut mat_b = crate::rlst_dynamic_mat![$Scalar, (4, 6)]; - let mut mat_c_actual = crate::rlst_row_vec![$Scalar, 6]; - let mut mat_c_expect = crate::rlst_row_vec![$Scalar, 6]; - - let dist = StandardNormal; - - let mut rng = rand::rngs::StdRng::seed_from_u64(0); - - mat_a.fill_from_equally_distributed(&mut rng); - mat_b.fill_from_equally_distributed(&mut rng); - mat_c_actual.fill_from_equally_distributed(&mut rng); - - for index in 0..mat_c_actual.layout().number_of_elements() { - *mat_c_expect.get1d_mut(index).unwrap() = - mat_c_actual.get1d_value(index).unwrap(); - } - - let alpha = <$Scalar>::random_scalar(&mut rng, &dist); - let beta = <$Scalar>::random_scalar(&mut rng, &dist); - - matmul_expect(alpha, &mat_a, &mat_b, beta, &mut mat_c_expect); - <$Scalar>::multiply_add(alpha, &mat_a, &mat_b, beta, &mut mat_c_actual); - - for index in 0..mat_c_expect.layout().number_of_elements() { - let val1 = mat_c_actual.get1d_value(index).unwrap(); - let val2 = mat_c_expect.get1d_value(index).unwrap(); - assert_ulps_eq!(&val1, &val2, max_ulps = 10); - } - } + fn matrix_multiply_compare( + transa: TransMode, + transb: TransMode, + alpha: Item, + mat_a: &Array, 2>, 2>, + mat_b: &Array, 2>, 2>, + beta: Item, + mat_c: &mut Array, 2>, 2>, + ) { + let a_shape = match transa { + TransMode::NoTrans => mat_a.shape(), + TransMode::ConjNoTrans => mat_a.shape(), + TransMode::Trans => [mat_a.shape()[1], mat_a.shape()[0]], + TransMode::ConjTrans => [mat_a.shape()[1], mat_a.shape()[0]], }; - } - - // matmul_test!(f32, test_matmul_f32); - matmul_test!(f64, test_matmul_f64); - matmul_test!(c32, test_matmul_c32); - matmul_test!(c64, test_matmul_c64); - - row_matvec_test!(f32, test_row_matvec_f32); - row_matvec_test!(f64, test_row_matvec_f64); - row_matvec_test!(c32, test_row_matvec_c32); - row_matvec_test!(c64, test_row_matvec_c64); - - col_matvec_test!(f64, test_col_matvec_f64); - col_matvec_test!(f32, test_col_matvec_f32); - col_matvec_test!(c32, test_col_matvec_c32); - col_matvec_test!(c64, test_col_matvec_c64); - #[test] - fn test_dot_matvec() { - let mut mat1 = crate::rlst_dynamic_mat![f64, (2, 2)]; - let mut mat2 = crate::rlst_dynamic_mat![f64, (2, 2)]; - - mat1.fill_from_seed_equally_distributed(0); - mat2.fill_from_seed_equally_distributed(1); + let b_shape = match transb { + TransMode::NoTrans => mat_b.shape(), + TransMode::ConjNoTrans => mat_b.shape(), + TransMode::Trans => [mat_b.shape()[1], mat_b.shape()[0]], + TransMode::ConjTrans => [mat_b.shape()[1], mat_b.shape()[0]], + }; - let mut mat3 = crate::rlst_dynamic_mat![f64, (2, 3)]; + let mut a_actual = rlst_dynamic_array2!(Item, a_shape); + let mut b_actual = rlst_dynamic_array2!(Item, b_shape); - mat3.fill_from_seed_equally_distributed(2); + match transa { + TransMode::NoTrans => a_actual.fill_from(mat_a.view()), + TransMode::ConjNoTrans => a_actual.fill_from(mat_a.view().conj()), + TransMode::Trans => a_actual.fill_from(mat_a.view().transpose()), + TransMode::ConjTrans => a_actual.fill_from(mat_a.view().conj().transpose()), + } - let actual = (mat1.view() + mat2.view()).dot(&mat3); - let mut expect = crate::rlst_dynamic_mat![f64, (2, 3)]; + match transb { + TransMode::NoTrans => b_actual.fill_from(mat_b.view()), + TransMode::ConjNoTrans => b_actual.fill_from(mat_b.view().conj()), + TransMode::Trans => b_actual.fill_from(mat_b.view().transpose()), + TransMode::ConjTrans => b_actual.fill_from(mat_b.view().conj().transpose()), + } - let mat_sum = (mat1.view() + mat2.view()).eval(); + let m = mat_c.shape()[0]; + let n = mat_c.shape()[1]; + let k = a_actual.shape()[1]; - for row in 0..2 { - for col in 0..3 { - for k in 0..2 { - expect[[row, col]] += mat_sum[[row, k]] * mat3[[k, col]]; + for row in 0..m { + for col in 0..n { + mat_c[[row, col]] *= beta; + for index in 0..k { + mat_c[[row, col]] += alpha * a_actual[[row, index]] * b_actual[[index, col]]; } } } - - assert_matrix_relative_eq!(expect, actual, 1E-14); } + + mat_mul_test_impl!(f64, 1E-14); + mat_mul_test_impl!(f32, 1E-5); + mat_mul_test_impl!(c32, 1E-5); + mat_mul_test_impl!(c64, 1E-14); }