From 9d198561d90b0f11d826a6b2cb45c1a79286b8df Mon Sep 17 00:00:00 2001 From: "Daniele G. Spampinato" Date: Thu, 22 Sep 2022 16:47:32 +0200 Subject: [PATCH] Fixes wrong symmetry constraints in mxm_band_generic --- include/alp/reference/blas3.hpp | 6 +- tests/unit/dense_mxm.cpp | 153 ++++++++++++++++++++++++++------ 2 files changed, 131 insertions(+), 28 deletions(-) diff --git a/include/alp/reference/blas3.hpp b/include/alp/reference/blas3.hpp index 8d57e2a20..3954c218b 100644 --- a/include/alp/reference/blas3.hpp +++ b/include/alp/reference/blas3.hpp @@ -236,7 +236,7 @@ namespace alp { auto &c_val = internal::access( C, internal::getStorageIndex( C, i, j ) ); // Size + Symmetry constraints - // sym_up_a * i <= l < K * (!sym_up_b) + (j+1) * (sym_up_b) + // sym_up_a * i <= l < K * (!sym_up_b) + ( j + 1 ) * (sym_up_b) // Band constraints // /\ i + l_a <= l < i + u_a // /\ j - u_b + 1 <= l < j - l_b + 1 @@ -268,11 +268,11 @@ namespace alp { auto &c_val = internal::access( C, internal::getStorageIndex( C, i, j ) ); // Size + Symmetry constraints - // max(sym_up_a * i, j) <= l < K + // max(sym_up_a * i, j + 1 ) <= l < K // Band constraints // /\ i + l_a <= l < i + u_a // /\ j - u_b + 1 <= l < j - l_b + 1 - for( std::ptrdiff_t l = std::max( { sym_up_a * i, j, i + l_a, j - u_b + 1 } ); + for( std::ptrdiff_t l = std::max( { sym_up_a * i, j + 1, i + l_a, j - u_b + 1 } ); l < std::min( { K, i + u_a, j - l_b + 1 } ); ++l ) { const auto ta { internal::access( A, internal::getStorageIndex( A, i, l ) ) }; diff --git a/tests/unit/dense_mxm.cpp b/tests/unit/dense_mxm.cpp index 2474807da..ff524c34b 100644 --- a/tests/unit/dense_mxm.cpp +++ b/tests/unit/dense_mxm.cpp @@ -72,20 +72,77 @@ void mxm_stdvec_as_matrix( std::vector< T > & vC, const size_t ldc, template< typename Structure, typename T > void stdvec_build_matrix( std::vector< T > & vA, const size_t m, const size_t n, const size_t lda, const T zero, const T one ) { - if( std::is_same< Structure, structures::General >::value ) { - std::fill( vA.begin(), vA.end(), one ); - } else if( std::is_same< Structure, structures::Symmetric >::value ) { - std::fill( vA.begin(), vA.end(), one ); - } else if( std::is_same< Structure, structures::UpperTriangular >::value ) { - for( size_t row = 0; row < m; ++row ) { - for( size_t col = 0; col < row; ++col ) { - vA[ row * lda + col ] = zero; + if( std::is_same< Structure, structures::General >::value ) { + std::fill( vA.begin(), vA.end(), one ); + } else if( std::is_same< Structure, structures::Symmetric >::value ) { + std::fill( vA.begin(), vA.end(), one ); + } else if( std::is_same< Structure, structures::UpperTriangular >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = 0; col < row; ++col ) { + vA[ row * lda + col ] = zero; + } + for( size_t col = row; col < n; ++col ) { + vA[ row * lda + col ] = one; + } } - for( size_t col = row; col < n; ++col ) { - vA[ row * lda + col ] = one; + } + +} + +template< typename Structure, typename T > +void stdvec_build_matrix( std::vector< T > & vA, const size_t m, const size_t n, const size_t lda, const T zero, const T one, const T inc ) { + + T val = one; + if( std::is_same< Structure, structures::General >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = 0; col < n; ++col ) { + vA[ row * lda + col ] = val; + val += inc; + } + } + } else if( std::is_same< Structure, structures::Symmetric >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = row; col < n; ++col ) { + vA[ row * lda + col ] = vA[ col * lda + row ] = val; + val += inc; + } + } + } else if( std::is_same< Structure, structures::UpperTriangular >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = 0; col < row; ++col ) { + vA[ row * lda + col ] = zero; + } + for( size_t col = row; col < n; ++col ) { + vA[ row * lda + col ] = val; + val += inc; + } + } + } + +} + +template< typename Structure, typename T > +void stdvec_build_matrix_packed( std::vector< T > & vA, const T one ) { + + std::fill( vA.begin(), vA.end(), one ); + +} + +template< typename Structure, typename T > +void stdvec_build_matrix_packed( std::vector< T > & vA, const T one, const T inc ) { + + T val = one; + if( std::is_same< Structure, structures::Symmetric >::value ) { // Assumes Packed Row - Upper + for( auto & elem: vA ) { + elem = val; + val += inc; + } + } else if( std::is_same< Structure, structures::UpperTriangular >::value ) { // Assumes Packed Row - Upper + for( auto & elem: vA ) { + elem = val; + val += inc; } } - } } @@ -141,11 +198,23 @@ void alp_program( const size_t & n, alp::RC & rc ) { T one = ring.getOne< T >(); T zero = ring.getZero< T >(); - std::vector< T > A_data( n * n, one ); - std::vector< T > B_data( n * n, one ); + std::vector< T > A_data( n * n ); + std::vector< T > B_data( n * n ); std::vector< T > C_data( n * n, zero ); + std::vector< T > A_packed( n * ( n + 1 ) / 2 ); + std::vector< T > B_packed( n * ( n + 1 ) / 2 ); + std::vector< T > C_packed( n * ( n + 1 ) / 2, zero ); + + std::vector< T > A_vec( n * n ); + std::vector< T > B_vec( n * n ); + std::vector< T > C_vec( n * n, zero ); + std::cout << "\tTesting dense General mxm " << n << std::endl; + + stdvec_build_matrix< structures::General >( A_data, n, n, n, zero, one, one ); + stdvec_build_matrix< structures::General >( B_data, n, n, n, zero, one, one ); + // initialize test alp::Matrix< T, structures::General > A( n, n ); alp::Matrix< T, structures::General > B( n, n ); @@ -164,9 +233,8 @@ void alp_program( const size_t & n, alp::RC & rc ) { print_matrix("C - POST", C); - std::vector< T > A_vec( n * n, one ); - std::vector< T > B_vec( n * n, one ); - std::vector< T > C_vec( n * n, zero ); + stdvec_build_matrix< structures::General >( A_vec, n, n, n, zero, one, one ); + stdvec_build_matrix< structures::General >( B_vec, n, n, n, zero, one, one ); mxm_stdvec_as_matrix( C_vec, n, A_vec, n, B_vec, n, n, n, n, ring.getMultiplicativeOperator(), ring.getAdditiveMonoid() ); @@ -178,17 +246,21 @@ void alp_program( const size_t & n, alp::RC & rc ) { alp::Matrix< T, structures::UpperTriangular > UB( n ); alp::Matrix< T, structures::UpperTriangular > UC( n ); - rc = alp::buildMatrix( UA, A_data.begin(), A_data.end() ); - rc = alp::buildMatrix( UB, B_data.begin(), B_data.end() ); - stdvec_build_matrix< structures::General >( C_data, n, n, n, zero, zero ); - rc = alp::buildMatrix( UC, C_data.begin(), C_data.end() ); + stdvec_build_matrix_packed< structures::UpperTriangular >( A_packed, one, one ); + stdvec_build_matrix_packed< structures::UpperTriangular >( B_packed, one, one ); + + rc = alp::buildMatrix( UA, A_packed.begin(), A_packed.end() ); + rc = alp::buildMatrix( UB, B_packed.begin(), B_packed.end() ); + rc = alp::buildMatrix( UC, C_packed.begin(), C_packed.end() ); + print_matrix("UA", UA); + print_matrix("UB", UB); print_matrix("UC - PRE", UC); rc = alp::mxm( UC, UA, UB, ring ); print_matrix("UC - POST", UC); - stdvec_build_matrix< structures::UpperTriangular >( A_vec, n, n, n, zero, one ); - stdvec_build_matrix< structures::UpperTriangular >( B_vec, n, n, n, zero, one ); + stdvec_build_matrix< structures::UpperTriangular >( A_vec, n, n, n, zero, one, one ); + stdvec_build_matrix< structures::UpperTriangular >( B_vec, n, n, n, zero, one, one ); stdvec_build_matrix< structures::General >( C_vec, n, n, n, zero, zero ); mxm_stdvec_as_matrix( C_vec, n, A_vec, n, B_vec, n, n, n, n, ring.getMultiplicativeOperator(), ring.getAdditiveMonoid() ); @@ -199,20 +271,51 @@ void alp_program( const size_t & n, alp::RC & rc ) { alp::Matrix< T, structures::Symmetric > SC( n ); - stdvec_build_matrix< structures::General >( C_data, n, n, n, zero, zero ); - rc = alp::buildMatrix( SC, C_data.begin(), C_data.end() ); + stdvec_build_matrix< structures::Symmetric >( A_data, n, n, n, zero, one, one ); + rc = alp::buildMatrix( A, A_data.begin(), A_data.end() ); + rc = alp::buildMatrix( SC, C_packed.begin(), C_packed.end() ); + + print_matrix("A", A ); + print_matrix("A^T", alp::get_view< alp::view::transpose >( A ) ); print_matrix("SC - PRE", SC); rc = alp::mxm( SC, A, alp::get_view< alp::view::transpose >( A ), ring ); print_matrix("SC - POST", SC); - stdvec_build_matrix< structures::General >( A_vec, n, n, n, zero, one ); + stdvec_build_matrix< structures::Symmetric >( A_vec, n, n, n, zero, one, one ); stdvec_build_matrix< structures::Symmetric >( C_vec, n, n, n, zero, zero ); mxm_stdvec_as_matrix( C_vec, n, A_vec, n, A_vec, n, n, n, n, ring.getMultiplicativeOperator(), ring.getAdditiveMonoid() ); diff_stdvec_matrix( C_vec, n, n, n, SC ); + std::cout << "\n\n=========== Testing Symmetric x Symmetric Output ============\n\n"; + + alp::Matrix< T, structures::Symmetric > SA( n ); + alp::Matrix< T, structures::Symmetric > SB( n ); + + stdvec_build_matrix_packed< structures::Symmetric >( A_packed, one, one ); + stdvec_build_matrix_packed< structures::Symmetric >( B_packed, one, one + one ); + + rc = alp::buildMatrix( SA, A_packed.begin(), A_packed.end() ); + rc = alp::buildMatrix( SB, B_packed.begin(), B_packed.end() ); + rc = alp::buildMatrix( C, C_data.begin(), C_data.end() ); + + print_matrix("SA", SA); + print_matrix("SB", SB); + print_matrix("C - PRE", C); + rc = alp::mxm( C, SA, SB, ring ); + print_matrix("C - POST", C); + + stdvec_build_matrix< structures::Symmetric >( A_vec, n, n, n, zero, one, one ); + stdvec_build_matrix< structures::Symmetric >( B_vec, n, n, n, zero, one, one + one ); + stdvec_build_matrix< structures::General >( C_vec, n, n, n, zero, zero ); + + mxm_stdvec_as_matrix( C_vec, n, A_vec, n, B_vec, n, n, n, n, ring.getMultiplicativeOperator(), ring.getAdditiveMonoid() ); + + diff_stdvec_matrix( C_vec, n, n, n, C ); + + } int main( int argc, char ** argv ) {