diff --git a/include/alp/reference/matrix.hpp b/include/alp/reference/matrix.hpp index 9937297b5..3d823e5a0 100644 --- a/include/alp/reference/matrix.hpp +++ b/include/alp/reference/matrix.hpp @@ -528,6 +528,8 @@ namespace alp { /** Expose static properties */ typedef T value_type; + typedef ImfR imf_r_type; + typedef ImfC imf_c_type; /** Type returned by access function */ typedef T &access_type; /** Type of the index used to access the physical storage */ @@ -799,7 +801,9 @@ namespace alp { public std::conditional< internal::is_view_over_functor< View >::value, internal::FunctorBasedMatrix< T, ImfR, ImfC, typename View::applied_to >, - internal::StorageBasedMatrix< T, ImfR, ImfC, storage::polynomials::Full_type, internal::requires_allocation< View >::value > + internal::StorageBasedMatrix< T, ImfR, ImfC, + typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type, + internal::requires_allocation< View >::value > >::type { protected: @@ -823,7 +827,7 @@ namespace alp { public: /** Exposes the types and the static properties. */ typedef structures::General structure; - typedef storage::polynomials::Full_type mapping_polynomial_type; + typedef typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type mapping_polynomial_type; /** * Indicates if a matrix needs to allocate data-related memory * (for the internal container or functor object). @@ -927,6 +931,26 @@ namespace alp { imf::Id( nrows ( target_matrix ) ), imf::Id( ncols ( target_matrix ) ) ) {} + /** + * Constructor for a view over another storage-based matrix. + * + * @tparam ViewType The dummy View type of the constructed matrix. + * Uses SFINAE to enable this constructor only for + * a view over a storage-based matrix. + */ + template< + typename ViewType = View, + std::enable_if_t< + internal::is_view_over_storage< ViewType >::value && + !internal::requires_allocation< ViewType >::value + > * = nullptr + > + Matrix( typename ViewType::applied_to &target_matrix, storage::AMF< ImfR, ImfC, mapping_polynomial_type > amf ) : + internal::StorageBasedMatrix< T, ImfR, ImfC, mapping_polynomial_type, requires_allocation >( + getContainer( target_matrix ), + amf + ) {} + /** * Constructor for a functor-based matrix that allocates memory. * @@ -993,7 +1017,9 @@ namespace alp { public std::conditional< internal::is_view_over_functor< View >::value, internal::FunctorBasedMatrix< T, ImfR, ImfC, typename View::applied_to >, - internal::StorageBasedMatrix< T, ImfR, ImfC, storage::polynomials::Full_type, internal::requires_allocation< View >::value > + internal::StorageBasedMatrix< T, ImfR, ImfC, + typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type, + internal::requires_allocation< View >::value > >::type { protected: @@ -1017,7 +1043,7 @@ namespace alp { public: /** Exposes the types and the static properties. */ typedef structures::Square structure; - typedef storage::polynomials::Full_type mapping_polynomial_type; + typedef typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type mapping_polynomial_type; /** * Indicates if a matrix needs to allocate data-related memory * (for the internal container or functor object). @@ -1109,6 +1135,26 @@ namespace alp { imf::Id( nrows ( target_matrix ) ), imf::Id( ncols ( target_matrix ) ) ) {} + /** + * Constructor for a view over another storage-based matrix. + * + * @tparam ViewType The dummy View type of the constructed matrix. + * Uses SFINAE to enable this constructor only for + * a view over a storage-based matrix. + */ + template< + typename ViewType = View, + std::enable_if_t< + internal::is_view_over_storage< ViewType >::value && + !internal::requires_allocation< ViewType >::value + > * = nullptr + > + Matrix( typename ViewType::applied_to &target_matrix, storage::AMF< ImfR, ImfC, mapping_polynomial_type > amf ) : + internal::StorageBasedMatrix< T, ImfR, ImfC, mapping_polynomial_type, requires_allocation >( + getContainer( target_matrix ), + amf + ) {} + /** * Constructor for a functor-based matrix that allocates memory. */ @@ -1163,7 +1209,9 @@ namespace alp { public std::conditional< internal::is_view_over_functor< View >::value, internal::FunctorBasedMatrix< T, ImfR, ImfC, typename View::applied_to >, - internal::StorageBasedMatrix< T, ImfR, ImfC, storage::polynomials::Full_type, internal::requires_allocation< View >::value > + internal::StorageBasedMatrix< T, ImfR, ImfC, + typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type, + internal::requires_allocation< View >::value > >::type { protected: @@ -1187,7 +1235,7 @@ namespace alp { public: /** Exposes the types and the static properties. */ typedef structures::Symmetric structure; - typedef storage::polynomials::Full_type mapping_polynomial_type; + typedef typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type mapping_polynomial_type; /** * Indicates if a matrix needs to allocate data-related memory * (for the internal container or functor object). @@ -1333,7 +1381,9 @@ namespace alp { public std::conditional< internal::is_view_over_functor< View >::value, internal::FunctorBasedMatrix< T, ImfR, ImfC, typename View::applied_to >, - internal::StorageBasedMatrix< T, ImfR, ImfC, storage::polynomials::Full_type, internal::requires_allocation< View >::value > + internal::StorageBasedMatrix< T, ImfR, ImfC, + typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type, + internal::requires_allocation< View >::value > >::type { protected: @@ -1357,7 +1407,7 @@ namespace alp { public: /** Exposes the types and the static properties. */ typedef structures::UpperTriangular structure; - typedef storage::polynomials::Full_type mapping_polynomial_type; + typedef typename storage::polynomials::apply_view< View::type_id, storage::polynomials::Full_type >::type mapping_polynomial_type; /** * Indicates if a matrix needs to allocate data-related memory * (for the internal container or functor object). @@ -1592,12 +1642,13 @@ namespace alp { * */ template< - typename T, typename Structure, enum Density density, typename View, typename ImfR, typename ImfC, enum Backend backend > - typename Matrix< T, Structure, density, View, ImfR, ImfC, backend >::template view_type< view::original >::type - get_view( Matrix< T, Structure, density, View, ImfR, ImfC, backend > & source ) { + typename SourceMatrixType, + std::enable_if_t< is_matrix< SourceMatrixType >::value, void > * = nullptr + > + typename SourceMatrixType::template view_type< view::original >::type + get_view( SourceMatrixType &source ) { - using source_strmat_t = Matrix< T, Structure, density, View, ImfR, ImfC, backend >; - using target_strmat_t = typename source_strmat_t::template view_type< view::original >::type; + using target_strmat_t = typename SourceMatrixType::template view_type< view::original >::type; target_strmat_t target( source ); @@ -1635,14 +1686,18 @@ namespace alp { */ template< enum view::Views target_view, - typename T, typename Structure, enum Density density, typename View, typename ImfR, typename ImfC, enum Backend backend > - typename Matrix< T, Structure, density, View, ImfR, ImfC, backend >::template view_type< target_view >::type - get_view( Matrix< T, Structure, density, View, ImfR, ImfC, backend > &source ) { + typename SourceMatrixType, + std::enable_if_t< is_matrix< SourceMatrixType >::value, void > * = nullptr + > + typename SourceMatrixType::template view_type< target_view >::type + get_view( SourceMatrixType &source ) { - using source_strmat_t = Matrix< T, Structure, density, View, ImfR, ImfC, backend >; - using target_strmat_t = typename source_strmat_t::template view_type< target_view >::type; + using target_strmat_t = typename SourceMatrixType::template view_type< target_view >::type; - target_strmat_t target( source ); + target_strmat_t target( + source, + storage::AMFFactory::Transform< target_view, decltype( source.amf ) >::Create( source.amf ) + ); return target; } diff --git a/include/alp/reference/vector.hpp b/include/alp/reference/vector.hpp index c255e3472..4d984d949 100644 --- a/include/alp/reference/vector.hpp +++ b/include/alp/reference/vector.hpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include @@ -465,6 +466,30 @@ namespace alp { Vector( typename ViewType::applied_to &vec_view ) : base_type( vec_view ) {} + /** + * Constructor for a view over another storage-based vector. + * + * @tparam ViewType The dummy View type of the constructed vector. + * Uses SFINAE to enable this constructor only for + * a view over a storage-based vector. + */ + template< + typename ViewType = View, + typename std::enable_if_t< + internal::is_view_over_storage< ViewType >::value && + !internal::requires_allocation< ViewType >::value + > * = nullptr + > + Vector( + typename ViewType::applied_to &vec_view, + storage::AMF< + typename base_type::imf_r_type, + typename base_type::imf_c_type, + typename base_type::mapping_polynomial_type + > amf + ) : + base_type( vec_view, amf ) {} + /** * Constructor for a functor-based vector that allocates memory. * diff --git a/include/alp/storage.hpp b/include/alp/storage.hpp index 4b954eac6..0adf54f4e 100644 --- a/include/alp/storage.hpp +++ b/include/alp/storage.hpp @@ -64,11 +64,23 @@ namespace alp { * @tparam A0 Static coefficient corresponding to constant term * @tparam Denominator Static denominator dividing the whole polynomial */ - template< size_t Ax2, size_t Ay2, size_t Axy, size_t Ax, size_t Ay, size_t A0, size_t Denominator > + template< + size_t coeffAx2, size_t coeffAy2, size_t coeffAxy, + size_t coeffAx, size_t coeffAy, + size_t coeffA0, + size_t Denominator + > struct BivariateQuadratic { static_assert( Denominator != 0, "Denominator cannot be zero (division by zero)."); + static constexpr size_t Ax2 = coeffAx2; + static constexpr size_t Ay2 = coeffAy2; + static constexpr size_t Axy = coeffAxy; + static constexpr size_t Ax = coeffAx; + static constexpr size_t Ay = coeffAy; + static constexpr size_t A0 = coeffA0; + static constexpr size_t D = Denominator; const size_t ax2, ay2, axy, ax, ay, a0; BivariateQuadratic( @@ -85,7 +97,7 @@ namespace alp { Axy * axy * x * y + Ax * ax * x + Ay * ay * y + - A0 * a0) / Denominator; + A0 * a0) / D; } }; // BivariateQuadratic @@ -109,6 +121,29 @@ namespace alp { return Full_type( 0, 0, 0, dim, 1, 0 ); } + template< enum view::Views view, typename Polynomial > + struct apply_view {}; + + template< typename Polynomial > + struct apply_view< view::Views::original, Polynomial > { + typedef Polynomial type; + }; + + template< typename Polynomial > + struct apply_view< view::Views::transpose, Polynomial > { + typedef BivariateQuadratic< Polynomial::Ay2, Polynomial::Ax2, Polynomial::Axy, Polynomial::Ay, Polynomial::Ax, Polynomial::A0, Polynomial::D > type; + }; + + template< typename Polynomial > + struct apply_view< view::Views::diagonal, Polynomial > { + typedef Polynomial type; + }; + + template< typename Polynomial > + struct apply_view< view::Views::_internal, Polynomial > { + typedef None_type type; + }; + }; // namespace polynomials /** @@ -149,6 +184,11 @@ namespace alp { private: + /** Expose static properties */ + typedef ImfR imf_r_type; + typedef ImfC imf_c_type; + typedef MappingPolynomial mapping_polynomial_type; + const ImfR imf_r; const ImfC imf_c; const MappingPolynomial map_poly; @@ -156,7 +196,7 @@ namespace alp { public: - AMF( ImfR &&imf_r, ImfC &&imf_c, MappingPolynomial map_poly, const size_t storage_dimensions ) : + AMF( ImfR imf_r, ImfC imf_c, MappingPolynomial map_poly, const size_t storage_dimensions ) : imf_r( imf_r ), imf_c( imf_c ), map_poly( map_poly ), storage_dimensions( storage_dimensions ) {} /** @@ -225,6 +265,11 @@ namespace alp { private: + /** Expose static properties */ + typedef imf::Strided imf_r_type; + typedef imf::Strided imf_c_type; + typedef MappingPolynomial mapping_polynomial_type; + /** For size checks */ const imf::Strided imf_r; const imf::Strided imf_c; @@ -327,6 +372,125 @@ namespace alp { ); } + /** + * Exposes the type of AMF resulting from applying the provided + * view on the provided AMF type. + */ + template< + enum view::Views view, + typename Amf, + typename ImfR = typename Amf::imf_r_type, + typename ImfC = typename Amf::imf_c_type, + typename MappingPolynomial = typename Amf::mapping_polynomial_type + > + struct apply_view {}; + + template< + typename Amf, + typename ImfR, + typename ImfC, + typename MappingPolynomial + > + struct apply_view< view::Views::original, Amf, ImfR, ImfC, MappingPolynomial > { + typedef AMF< + ImfR, + ImfC, + typename polynomials::apply_view< view::Views::original, MappingPolynomial >::type + > type; + }; + + template< + typename Amf, + typename ImfR, + typename ImfC, + typename MappingPolynomial + > + struct apply_view< view::Views::transpose, Amf, ImfR, ImfC, MappingPolynomial > { + typedef AMF< + ImfC, + ImfR, + typename polynomials::apply_view< view::Views::transpose, MappingPolynomial >::type + > type; + }; + + template< + typename Amf, + typename ImfR, + typename ImfC, + typename MappingPolynomial + > + struct apply_view< view::Views::diagonal, Amf, ImfR, ImfC, MappingPolynomial > { + typedef AMF< + ImfR, + ImfC, + typename polynomials::apply_view< view::Views::diagonal, MappingPolynomial >::type + > type; + }; + + template< enum view::Views view, typename AMFType > + struct Transform { + + static + AMFType + Create( const AMFType &amf ) { + throw std::invalid_argument( "Not implemented for the provided view type." ); + return amf; + } + + }; + + template< typename AMFType > + struct Transform< view::Views::original, AMFType > { + + static + AMFType + Create( const AMFType &amf ) { + return amf; + } + + }; + + template< typename AMFType > + struct Transform< view::Views::transpose, AMFType > { + + static + AMF< + typename AMFType::imf_c_type, + typename AMFType::imf_r_type, + typename polynomials::apply_view< + view::Views::transpose, + typename AMFType::mapping_polynomial_type + >::type + > + Create( const AMFType &amf ) { + typedef typename polynomials::apply_view< view::Views::transpose, typename AMFType::mapping_polynomial_type >::type new_mapping_polynomial_type; + return AMF< + typename AMFType::imf_c_type, + typename AMFType::imf_r_type, + new_mapping_polynomial_type + >( + amf.imf_c, + amf.imf_r, + new_mapping_polynomial_type( + amf.map_poly.ay2, amf.map_poly.ax2, amf.map_poly.axy, + amf.map_poly.ay, amf.map_poly.ax, + amf.map_poly.a0 + ), + amf.storage_dimensions + ); + } + }; + + template< typename AMFType > + struct Transform< view::Views::diagonal, AMFType > { + + static + AMFType + Create( const AMFType &amf ) { + return amf; + } + }; + }; // class AMFFactory }; // namespace storage diff --git a/include/alp/views.hpp b/include/alp/views.hpp index db5dd5f7e..f962ef898 100644 --- a/include/alp/views.hpp +++ b/include/alp/views.hpp @@ -36,16 +36,31 @@ namespace alp { namespace view { + /** + * Lists the view types exposed to the user. + * + * \note View type "_internal" shall not be used by the user and + * its use may result in an undefined behaviour. + * + * \note \internal "_internal" value is added so that all view types have + * a defined type_id field, which is used by internal + * type traits. + * + */ enum Views { original, transpose, - diagonal + diagonal, + _internal }; + template< typename OriginalType > struct Original { using applied_to = OriginalType; + static constexpr Views type_id = Views::original; + static std::pair< size_t, size_t > dims( std::pair< size_t, size_t > dims_pair ) { return std::make_pair( dims_pair.first, dims_pair.second ); } @@ -56,6 +71,8 @@ namespace alp { using applied_to = OriginalType; + static constexpr Views type_id = Views::transpose; + static std::pair< size_t, size_t > dims( std::pair< size_t, size_t > dims_pair ) { return std::make_pair( dims_pair.second, dims_pair.first ); } @@ -66,6 +83,8 @@ namespace alp { using applied_to = OriginalType; + static constexpr Views type_id = Views::diagonal; + static size_t getLength( std::pair< size_t, size_t > dims_pair ) { return std::min( dims_pair.first, dims_pair.second ); } @@ -76,6 +95,9 @@ namespace alp { using applied_to = LambdaFunctionType; + /** Functor views are not exposed to the user */ + static constexpr Views type_id = Views::_internal; + static std::pair< size_t, size_t > getLength( std::pair< size_t, size_t > dims_pair ) { return std::make_pair( dims_pair.first, dims_pair.second ); } diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 4bff2bab4..eceae4497 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -34,6 +34,10 @@ add_grb_executables( add15m add15m.cpp # BACKENDS alp_reference #) +add_grb_executables( alp_views alp_views.cpp + BACKENDS alp_reference +) + add_grb_executables( alp_type_traits alp_type_traits.cpp BACKENDS alp_reference ) diff --git a/tests/unit/alp_views.cpp b/tests/unit/alp_views.cpp new file mode 100644 index 000000000..60161877b --- /dev/null +++ b/tests/unit/alp_views.cpp @@ -0,0 +1,205 @@ + +/* + * 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 "alp.hpp" + +template< typename MatrixType > +void print_matrix( std::string name, const MatrixType &A ) { + + if( ! alp::internal::getInitialized( A ) ) { + std::cout << "Matrix " << name << " uninitialized.\n"; + return; + } + + std::cout << "Matrix " << name << " of size " << alp::dims( A ).first << " x " << alp::dims( A ).second << " contains the following elements:\n"; + + for( size_t row = 0; row < alp::nrows( A ); ++row ) { + std::cout << "[\t"; + for( size_t col = 0; col < alp::ncols( A ); ++col ) { + auto pos = alp::internal::getStorageIndex( A, row, col ); + // std::cout << "(" << pos << "): "; + std::cout << alp::internal::access( A, pos ) << "\t"; + } + std::cout << "]\n"; + } +} + +template< typename VectorType > +void print_vector( std::string name, const VectorType &v ) { + + if( ! alp::internal::getInitialized( v ) ) { + std::cout << "Vector " << name << " uninitialized.\n"; + return; + } + + std::cout << "Vector " << name << " of size " << alp::dims( v ).first << " contains the following elements:\n"; + + std::cout << "["; + for( size_t row = 0; row < alp::nrows( v ); ++row ) { + std::cout << v[ row ] << "\t"; + } + std::cout << "]\n"; +} + +template< typename T > +void init_matrix( std::vector< T > &A, const size_t rows, const size_t cols ) { + + for( size_t row = 0; row < rows; ++row ) { + for( size_t col = 0; col < cols; ++col ) { + A[ row * cols + col ] = row + col; + } + } +} + +template< typename T > +void print_stdvec_as_matrix( std::string name, const std::vector< T > &vA, const size_t m, const size_t n, const size_t lda ) { + + std::cout << "Vec " << name << ":" << std::endl; + for( size_t row = 0; row < m; ++row ) { + std::cout << "[\t"; + for( size_t col = 0; col < n; ++col ) { + std::cout << vA[ row * lda + col ] << "\t"; + } + std::cout << "]" << std::endl; + } +} + +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, alp::structures::General >::value ) { + std::fill( vA.begin(), vA.end(), one ); + } else if( std::is_same< Structure, alp::structures::Symmetric >::value ) { + std::fill( vA.begin(), vA.end(), one ); + } +} + +template< typename MatType, typename T > +void diff_stdvec_matrix( + const std::vector< T > &vA, const size_t m, const size_t n, const size_t lda, + const MatType &mA, + double threshold=1e-7 +) { + + if( std::is_same< typename MatType::structure, alp::structures::General >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = 0; col < n; ++col ) { + double va = ( double )( vA[ row * lda + col ] ); + double vm = ( double )( alp::internal::access( mA, alp::internal::getStorageIndex( mA, row, col ) ) ); + double re = std::abs( ( va - vm ) / va ); + if( re > threshold ) { + std::cout << "Error ( " << row << ", " << col << " ): " << va << " v " << vm << std::endl; + } + } + } + } else if( std::is_same< typename MatType::structure, alp::structures::Symmetric >::value ) { + for( size_t row = 0; row < m; ++row ) { + for( size_t col = row; col < n; ++col ) { + double va = ( double )( vA[ row * lda + col ] ); + double vm = ( double )( alp::internal::access( mA, alp::internal::getStorageIndex( mA, row, col ) ) ); + double re = std::abs( ( va - vm ) / va ); + if( re > threshold ) { + std::cout << "Error ( " << row << ", " << col << " ): " << va << " v " << vm << std::endl; + } + } + } + } +} + + + +// alp program +void alpProgram( const size_t &n, alp::RC &rc ) { + + typedef double T; + + alp::Semiring< alp::operators::add< T >, alp::operators::mul< T >, alp::identities::zero, alp::identities::one > ring; + + T zero = ring.getZero< T >(); + + // allocate + const size_t m = 2 * n; + std::vector< T > M_data( m * n, zero ); + init_matrix( M_data, m, n ); + + alp::Matrix< T, alp::structures::General > M( m, n ); + alp::buildMatrix( M, M_data.begin(), M_data.end() ); + print_matrix( "M", M ); + + auto MT = alp::get_view< alp::view::Views::transpose >( M ); + print_matrix( "M^T", MT ); + + // diable diagonal view test as diagonal views are not operational yet + //auto Mdiag = alp::get_view< alp::view::Views::diagonal >( M ); + //print_vector( "Mdiag", Mdiag ); + + rc = alp::SUCCESS; + +} + +int main( int argc, char ** argv ) { + // defaults + bool printUsage = false; + size_t in = 100; + + // 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 == 0 ) { + std::cerr << "n must be a positive number\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 integer, the " + "test size.\n"; + return 1; + } + std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; + + std::cout << "This is functional test " << argv[ 0 ] << "\n"; + alp::Launcher< alp::AUTOMATIC > launcher; + alp::RC out; + if( launcher.exec( &alpProgram, 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; +} +