From 2eb12b054866f5b27ae16d015e0f4823c48565cd Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 14:30:49 +0200 Subject: [PATCH 1/7] Add conjugate operation with specializations for square and non-square matrices --- include/alp/reference/blas2.hpp | 113 ++++++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/include/alp/reference/blas2.hpp b/include/alp/reference/blas2.hpp index 085b33c11..b1e522942 100644 --- a/include/alp/reference/blas2.hpp +++ b/include/alp/reference/blas2.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #define NO_CAST_OP_ASSERT( x, y, z ) \ static_assert( x, \ @@ -1054,6 +1055,118 @@ namespace alp { return internal::fold_matrix_generic< left, scalar, descr >( &A, no_matrix, &beta, op ) ; } + /** + * Returns a view over the input matrix returning conjugate of the accessed element. + * This avoids materializing the resulting container. + * The elements are calculated lazily on access. + * + * @tparam descr The descriptor to be used (descriptors::no_operation + * if left unspecified). + * @tparam InputType The value type of the input matrix. + * @tparam InputStructure The Structure type applied to the input matrix. + * @tparam InputView The view type applied to the input matrix. + * + * @param A The input matrix + * + * @return Matrix view over a lambda function defined in this function. + * + * Specialization for non-square matrices. This distinction is necessary due + * to different constructor signature for square and non-square matrices. + * + */ + template< + Descriptor descr = descriptors::no_operation, + typename DataType, typename Structure, typename View, typename ImfR, typename ImfC, + std::enable_if_t< + !structures::is_in< structures::Square, typename Structure::inferred_structures >::value + > * = nullptr + > + Matrix< + DataType, Structure, Density::Dense, + view::Functor< std::function< void( DataType &, const size_t, const size_t ) > >, + imf::Id, imf::Id, + reference + > + conjugate( + const Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A, + const std::enable_if_t< + !alp::is_object< DataType >::value + > * const = nullptr + ) { + + std::function< void( DataType &, const size_t, const size_t ) > data_lambda = + [ &A ]( DataType &result, const size_t i, const size_t j ) { + result = grb::utils::is_complex< DataType >::conjugate( + internal::access( A, internal::getStorageIndex( A, i, j ) ) + ); + }; + std::function< bool() > init_lambda = + [ &A ]() -> bool { + return internal::getInitialized( A ); + }; + + return Matrix< + DataType, + Structure, + Density::Dense, + view::Functor< std::function< void( DataType &, const size_t, const size_t ) > >, + imf::Id, imf::Id, + reference + >( + init_lambda, + nrows( A ), + ncols( A ), + data_lambda + ); + + } + + /** Specialization for square matrices */ + template< + Descriptor descr = descriptors::no_operation, + typename DataType, typename Structure, typename View, typename ImfR, typename ImfC, + std::enable_if_t< + structures::is_in< structures::Square, typename Structure::inferred_structures >::value + > * = nullptr + > + Matrix< + DataType, Structure, Density::Dense, + view::Functor< std::function< void( DataType &, const size_t, const size_t ) > >, + imf::Id, imf::Id, + reference + > + conjugate( + const Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A, + const std::enable_if_t< + !alp::is_object< DataType >::value + > * const = nullptr + ) { + + std::function< void( DataType &, const size_t, const size_t ) > data_lambda = + [ &A ]( DataType &result, const size_t i, const size_t j ) { + result = grb::utils::is_complex< DataType >::conjugate( + internal::access( A, internal::getStorageIndex( A, i, j ) ) + ); + }; + std::function< bool() > init_lambda = + [ &A ]() -> bool { + return internal::getInitialized( A ); + }; + + return Matrix< + DataType, + Structure, + Density::Dense, + view::Functor< std::function< void( DataType &, const size_t, const size_t ) > >, + imf::Id, imf::Id, + reference + >( + init_lambda, + nrows( A ), + data_lambda + ); + + } /** @} */ } // end namespace ``alp'' From 880dc1bffa2d2954fd2396f57108a142e23b97e2 Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 14:31:20 +0200 Subject: [PATCH 2/7] Add unit test for conjugate It checks the correctness of conjugate by comparing the result with the tranpose of the original matrix. Tests complex and real matrices. --- tests/unit/CMakeLists.txt | 4 + tests/unit/dense_conjugate.cpp | 194 +++++++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 tests/unit/dense_conjugate.cpp diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index ad00f4166..69dc5f637 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -224,6 +224,10 @@ add_grb_executables( spy spy.cpp BACKENDS reference reference_omp ) +add_grb_executables( dense_conjugate dense_conjugate.cpp + BACKENDS alp_reference +) + add_grb_executables( dense_matrix_access dense_matrix_access.cpp BACKENDS alp_reference ) diff --git a/tests/unit/dense_conjugate.cpp b/tests/unit/dense_conjugate.cpp new file mode 100644 index 000000000..909671db7 --- /dev/null +++ b/tests/unit/dense_conjugate.cpp @@ -0,0 +1,194 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include "../utils/print_alp_containers.hpp" + +constexpr float tol = 1.e-10; + +template< typename T > +T random_value(); + +template<> +float random_value< float >() { + return static_cast< float >( rand() ) / RAND_MAX; +} + +template<> +std::complex< float > random_value< std::complex< float > >() { + const float re = random_value< float >(); + const float im = random_value< float >(); + return std::complex< float >( re, im ); +} + +template< typename MatrixType > +void init_matrix( MatrixType &M ) { + // Temporary until proper matrix building is implemented + typedef typename MatrixType::value_type value_type; + alp::internal::setInitialized( M, true ); + const size_t height = alp::ncols( M ); + const size_t width = alp::nrows( M ); + for( size_t r = 0; r < height; ++r ) { + for( size_t c = 0; c < width; ++c ) { + const value_type val = random_value< value_type >(); + if( r < c ) { + alp::internal::access( M, alp::internal::getStorageIndex( M, r, c ) ) = val; + if( r != c ) { + alp::internal::access( M, alp::internal::getStorageIndex( M, c, r ) ) = grb::utils::is_complex< value_type >::conjugate( val ); + } + } else if ( r == c ) { + alp::internal::access( M, alp::internal::getStorageIndex( M, r, c ) ) = std::real( val ); + } + } + } +} + +template< + typename MatrixType1, + typename MatrixType2, + typename T = typename MatrixType1::value_type, + typename Ring +> +alp::RC check_if_same( const MatrixType1 &A, const MatrixType2 &B, const Ring &ring ) { + + alp::RC rc = alp::SUCCESS; + + alp::Matrix< T, alp::structures::Square > E( nrows( A ) ); + rc = rc ? rc : set( E, alp::Scalar< T >( ring.template getZero< T >() ) ); + + rc = rc ? rc : alp::foldl( E, A, ring.getAdditiveOperator() ); + rc = rc ? rc : alp::foldl( E, B, alp::operators::subtract< T >() ); + + float fnorm = ring.template getZero< float >(); + rc = rc ? rc : alp::eWiseLambda( + [ &fnorm, &ring ]( const size_t i, const size_t j, T &val ) { + (void) i; + (void) j; + const float valsquare = static_cast< float >( std::real( val * grb::utils::is_complex< T >::conjugate( val ) ) ); + alp::internal::foldl( + fnorm, + valsquare, + alp::operators::add< float >() + ); + }, + E + ); + fnorm = std::sqrt( fnorm ); + + if( fnorm < tol ) { + return alp::SUCCESS; + } else { + return alp::FAILED; + } +} + +template< + typename T, + typename Structure = typename std::conditional< + grb::utils::is_complex< T >::value, + alp::structures::Hermitian, + alp::structures::Square + >::type +> +alp::RC test_conjugate( const size_t n ) { + + const alp::Semiring< alp::operators::add< T >, alp::operators::mul< T >, alp::identities::zero, alp::identities::one > ring; + + alp::RC rc = alp::SUCCESS; + + // create the original matrix + alp::Matrix< T, Structure > H( n, n ); + // set matrix elements using the internal interface + init_matrix( H ); + + // create a conjugated matrix + auto H_conj = alp::conjugate( H ); + + // create a transpose view over original matrix (used for error checking) + auto H_T = alp::get_view< alp::view::transpose >( H ); + + // check if conjugated and transposed matrix are the same + rc = rc ? rc : check_if_same( H_conj, H_T, ring ); + + return rc; +} + +void alp_program( const size_t &n, alp::RC &rc ) { + + rc = alp::SUCCESS; + + rc = rc ? rc : test_conjugate< std::complex< float > >( n ); + rc = rc ? rc : test_conjugate< float >( n ); + + rc = alp::SUCCESS; +} + +int main( int argc, char ** argv ) { + // defaults + bool printUsage = false; + size_t in = 5; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( read % 2 != 0 ) { + std::cerr << "Given value for n is odd\n"; + printUsage = true; + } else { + // all OK + in = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an even integer, the " + "test size.\n"; + return 1; + } + + std::cout << "This is functional test " << argv[ 0 ] << "\n"; + alp::Launcher< alp::AUTOMATIC > launcher; + alp::RC out; + if( launcher.exec( &alp_program, in, out, true ) != alp::SUCCESS ) { + std::cerr << "Launching test FAILED\n"; + return 255; + } + if( out != alp::SUCCESS ) { + std::cerr << "Test FAILED (" << alp::toString( out ) << ")" << std::endl; + } else { + std::cout << "Test OK" << std::endl; + } + return 0; +} From 000829e31792efe7400acc8b6ee1be21a55e2bb0 Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 17:22:28 +0200 Subject: [PATCH 3/7] Add an explanation for a temporary use of square instead of symmetric structure for real matrices --- tests/unit/dense_conjugate.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit/dense_conjugate.cpp b/tests/unit/dense_conjugate.cpp index 909671db7..6931a98b6 100644 --- a/tests/unit/dense_conjugate.cpp +++ b/tests/unit/dense_conjugate.cpp @@ -109,6 +109,8 @@ template< typename Structure = typename std::conditional< grb::utils::is_complex< T >::value, alp::structures::Hermitian, + // Should be Symmetric. + // Temporarily using Square until fold is fixed to support folding symmetric onto more general structures alp::structures::Square >::type > From e7139c644663d68a870dc3fcd7d90b19882d03aa Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 17:30:14 +0200 Subject: [PATCH 4/7] Use a dedicated type trait for the SFINAE check --- include/alp/reference/blas2.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/alp/reference/blas2.hpp b/include/alp/reference/blas2.hpp index b1e522942..4da796d0d 100644 --- a/include/alp/reference/blas2.hpp +++ b/include/alp/reference/blas2.hpp @@ -1078,7 +1078,7 @@ namespace alp { Descriptor descr = descriptors::no_operation, typename DataType, typename Structure, typename View, typename ImfR, typename ImfC, std::enable_if_t< - !structures::is_in< structures::Square, typename Structure::inferred_structures >::value + !structures::is_a< Structure, structures::Square >::value > * = nullptr > Matrix< @@ -1126,7 +1126,7 @@ namespace alp { Descriptor descr = descriptors::no_operation, typename DataType, typename Structure, typename View, typename ImfR, typename ImfC, std::enable_if_t< - structures::is_in< structures::Square, typename Structure::inferred_structures >::value + structures::is_a< Structure, structures::Square >::value > * = nullptr > Matrix< From 54eb2b02f2d25b1c7a2404574c97b52ac7237afd Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 17:33:12 +0200 Subject: [PATCH 5/7] Replace float with a typedef defined at the beginning of the unit test --- tests/unit/dense_conjugate.cpp | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/tests/unit/dense_conjugate.cpp b/tests/unit/dense_conjugate.cpp index 6931a98b6..94e3612bd 100644 --- a/tests/unit/dense_conjugate.cpp +++ b/tests/unit/dense_conjugate.cpp @@ -26,21 +26,22 @@ #include #include "../utils/print_alp_containers.hpp" -constexpr float tol = 1.e-10; +typedef float BaseScalarType; +constexpr BaseScalarType tol = 1.e-10; template< typename T > T random_value(); template<> -float random_value< float >() { - return static_cast< float >( rand() ) / RAND_MAX; +BaseScalarType random_value< BaseScalarType >() { + return static_cast< BaseScalarType >( rand() ) / RAND_MAX; } template<> -std::complex< float > random_value< std::complex< float > >() { - const float re = random_value< float >(); - const float im = random_value< float >(); - return std::complex< float >( re, im ); +std::complex< BaseScalarType > random_value< std::complex< BaseScalarType > >() { + const BaseScalarType re = random_value< BaseScalarType >(); + const BaseScalarType im = random_value< BaseScalarType >(); + return std::complex< BaseScalarType >( re, im ); } template< typename MatrixType > @@ -81,16 +82,16 @@ alp::RC check_if_same( const MatrixType1 &A, const MatrixType2 &B, const Ring &r rc = rc ? rc : alp::foldl( E, A, ring.getAdditiveOperator() ); rc = rc ? rc : alp::foldl( E, B, alp::operators::subtract< T >() ); - float fnorm = ring.template getZero< float >(); + BaseScalarType fnorm = ring.template getZero< BaseScalarType >(); rc = rc ? rc : alp::eWiseLambda( [ &fnorm, &ring ]( const size_t i, const size_t j, T &val ) { (void) i; (void) j; - const float valsquare = static_cast< float >( std::real( val * grb::utils::is_complex< T >::conjugate( val ) ) ); + const BaseScalarType valsquare = static_cast< BaseScalarType >( std::real( val * grb::utils::is_complex< T >::conjugate( val ) ) ); alp::internal::foldl( fnorm, valsquare, - alp::operators::add< float >() + alp::operators::add< BaseScalarType >() ); }, E @@ -141,8 +142,8 @@ void alp_program( const size_t &n, alp::RC &rc ) { rc = alp::SUCCESS; - rc = rc ? rc : test_conjugate< std::complex< float > >( n ); - rc = rc ? rc : test_conjugate< float >( n ); + rc = rc ? rc : test_conjugate< std::complex< BaseScalarType > >( n ); + rc = rc ? rc : test_conjugate< BaseScalarType >( n ); rc = alp::SUCCESS; } From aacfba1f8181ceac2c5d9681d0955f6f5d15e381 Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 17:34:26 +0200 Subject: [PATCH 6/7] Use std::norm instead of manually calculating it --- tests/unit/dense_conjugate.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/dense_conjugate.cpp b/tests/unit/dense_conjugate.cpp index 94e3612bd..c53282d4d 100644 --- a/tests/unit/dense_conjugate.cpp +++ b/tests/unit/dense_conjugate.cpp @@ -87,7 +87,7 @@ alp::RC check_if_same( const MatrixType1 &A, const MatrixType2 &B, const Ring &r [ &fnorm, &ring ]( const size_t i, const size_t j, T &val ) { (void) i; (void) j; - const BaseScalarType valsquare = static_cast< BaseScalarType >( std::real( val * grb::utils::is_complex< T >::conjugate( val ) ) ); + const BaseScalarType valsquare = std::norm( val ); alp::internal::foldl( fnorm, valsquare, From 4d0f6ae8b7fafb7226c785898229970238e21db1 Mon Sep 17 00:00:00 2001 From: Vladimir Dimic Date: Fri, 30 Sep 2022 17:36:06 +0200 Subject: [PATCH 7/7] Remove the override for rc value --- tests/unit/dense_conjugate.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit/dense_conjugate.cpp b/tests/unit/dense_conjugate.cpp index c53282d4d..752b00901 100644 --- a/tests/unit/dense_conjugate.cpp +++ b/tests/unit/dense_conjugate.cpp @@ -145,7 +145,6 @@ void alp_program( const size_t &n, alp::RC &rc ) { rc = rc ? rc : test_conjugate< std::complex< BaseScalarType > >( n ); rc = rc ? rc : test_conjugate< BaseScalarType >( n ); - rc = alp::SUCCESS; } int main( int argc, char ** argv ) {