diff --git a/include/alp/reference/blas1.hpp b/include/alp/reference/blas1.hpp index b6a61de4f..4985c908a 100644 --- a/include/alp/reference/blas1.hpp +++ b/include/alp/reference/blas1.hpp @@ -23,11 +23,17 @@ #ifndef _H_ALP_REFERENCE_BLAS1 #define _H_ALP_REFERENCE_BLAS1 +#include + #include #include #include #include #include +#include +#include +#include + #ifndef NO_CAST_ASSERT #define NO_CAST_ASSERT( x, y, z ) \ @@ -817,9 +823,8 @@ namespace alp { Scalar< IOType, IOStructure, reference > & beta, const Monoid & monoid = Monoid(), const typename std::enable_if< ! alp::is_object< InputType >::value && ! alp::is_object< IOType >::value && alp::is_monoid< Monoid >::value, void >::type * const = NULL ) { - - throw std::runtime_error( "Needs an implementation." ); - return SUCCESS; + + return foldl( beta, x, monoid ); } /** C++ scalar variant */ @@ -927,8 +932,7 @@ namespace alp { "called on a vector x of a type that does not match the third domain " "of the given operator" ); - throw std::runtime_error( "Needs an implementation." ); - return SUCCESS; + return foldl( y, alpha, monoid ); } /** @@ -1332,19 +1336,37 @@ namespace alp { RC foldl( Vector< IOType, IOStructure, Density::Dense, IOView, IOImfR, IOImfC, reference > & x, const Scalar< InputType, InputStructure, reference > beta, const Op & op = Op(), - const typename std::enable_if< ! alp::is_object< IOType >::value && ! alp::is_object< InputType >::value && alp::is_operator< Op >::value, void >::type * = NULL ) { + //const typename std::enable_if< ! alp::is_object< IOType >::value && ! alp::is_object< InputType >::value && alp::is_operator< Op >::value, void >::type * = NULL ) { + const typename std::enable_if< + ! alp::is_object< IOType >::value && ! alp::is_object< InputType >::value && alp::is_operator< Op >::value, void + >::type * = NULL ) { // static sanity checks - NO_CAST_OP_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< typename Op::D1, IOType >::value ), "alp::foldl", + NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) + || std::is_same< typename Op::D1, IOType >::value ), "alp::foldl", "called with a vector x of a type that does not match the first domain " "of the given operator" ); - NO_CAST_OP_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< typename Op::D2, InputType >::value ), "alp::foldl", + NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) + || std::is_same< typename Op::D2, InputType >::value ), "alp::foldl", "called on a vector y of a type that does not match the second domain " "of the given operator" ); - NO_CAST_OP_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< typename Op::D3, IOType >::value ), "alp::foldl", + NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) + || std::is_same< typename Op::D3, IOType >::value ), "alp::foldl", "called on a vector x of a type that does not match the third domain " "of the given operator" ); - throw std::runtime_error( "Needs an implementation." ); + internal::setInitialized( + x, + internal::getInitialized( x ) && internal::getInitialized( beta ) + ); + + if( !internal::getInitialized( x ) ) { + return SUCCESS; + } + + const size_t n = getLength( x ); + for ( size_t i = 0; i < n ; i++ ) { + (void) internal::foldl( x[ i ], *beta, op ); + } return SUCCESS; } @@ -4185,19 +4207,13 @@ namespace alp { template< Descriptor descr = descriptors::no_operation, typename IOType, typename IOStructure, typename InputType, typename InputStructure, typename InputView, typename InputImfR, typename InputImfC, - typename MaskType, typename MaskStructure, typename MaskView, typename MaskImfR, typename MaskImfC, class Monoid > - RC foldl( Scalar< IOType, IOStructure, reference > &x, + RC foldl( Scalar< IOType, IOStructure, reference > &alpha, const Vector< InputType, InputStructure, Density::Dense, InputView, InputImfR, InputImfC, reference > & y, - const Vector< MaskType, MaskStructure, Density::Dense, MaskView, MaskImfR, MaskImfC, reference > & mask, const Monoid & monoid = Monoid(), - const typename std::enable_if< ! alp::is_object< IOType >::value && ! alp::is_object< InputType >::value && ! alp::is_object< MaskType >::value && alp::is_monoid< Monoid >::value, - void >::type * const = NULL ) { - #ifdef _DEBUG - std::cout << "foldl: IOType <- [InputType] with a monoid called. Array has size " << size( y ) << " with " << nnz( y ) << " nonzeroes. It has a mask of size " << size( mask ) << " with " - << nnz( mask ) << " nonzeroes.\n"; - #endif + const typename std::enable_if< ! alp::is_object< IOType >::value && ! alp::is_object< InputType >::value && alp::is_monoid< Monoid >::value, + void >::type * const = NULL ) { // static sanity checks NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< IOType, InputType >::value ), "alp::reduce", "called with a scalar IO type that does not match the input vector type" ); @@ -4210,9 +4226,21 @@ namespace alp { NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< InputType, typename Monoid::D3 >::value ), "alp::reduce", "called with an input vector type that does not match the third domain of " "the given monoid" ); - NO_CAST_ASSERT( ( ! ( descr & descriptors::no_casting ) || std::is_same< bool, MaskType >::value ), "alp::reduce", "called with a vector mask type that is not boolean" ); + NO_CAST_ASSERT( ! ( descr & descriptors::no_casting ), "alp::reduce", "called with a vector mask type that is not boolean" ); - throw std::runtime_error( "Needs an implementation." ); + internal::setInitialized( + alpha, + internal::getInitialized( alpha ) && internal::getInitialized( y ) + ); + + if( !internal::getInitialized( alpha ) ) { + return SUCCESS; + } + + const size_t n = getLength( y ); + for ( size_t i = 0; i < n ; i++ ) { + (void) internal::foldl( *alpha, y[ i ], monoid.getOperator() ); + } return SUCCESS; } diff --git a/include/alp/reference/io.hpp b/include/alp/reference/io.hpp index f8d4c0e14..7e6418bd1 100644 --- a/include/alp/reference/io.hpp +++ b/include/alp/reference/io.hpp @@ -189,7 +189,6 @@ namespace alp { const fwd_iterator &start, const fwd_iterator &end ) noexcept { - // Temporarily assuming 1-1 mapping with user container internal::setInitialized(v, true); @@ -202,7 +201,7 @@ namespace alp { *p = *it; } - return PANIC; + return SUCCESS; } } // end namespace ``alp'' diff --git a/include/alp/reference/scalar.hpp b/include/alp/reference/scalar.hpp index 6b67bdb50..1f72b4fcf 100644 --- a/include/alp/reference/scalar.hpp +++ b/include/alp/reference/scalar.hpp @@ -37,7 +37,7 @@ namespace alp { namespace internal { template< typename T, typename Structure > - bool getInitialized( Scalar< T, Structure, reference > & ) noexcept; + bool getInitialized( const Scalar< T, Structure, reference > & ) noexcept; template< typename T, typename Structure > void setInitialized( Scalar< T, Structure, reference > &, bool ) noexcept; @@ -60,10 +60,12 @@ namespace alp { */ template< typename T, typename Structure > class Scalar< T, Structure, reference > { + private: + typedef Scalar< T, Structure, reference > self_type; - friend bool internal::getInitialized<>( self_type & ) noexcept; + friend bool internal::getInitialized<>( const self_type & ) noexcept; friend void internal::setInitialized<>( self_type &, bool ) noexcept; @@ -79,6 +81,7 @@ namespace alp { /** @see Vector::lambda_reference */ typedef T& lambda_reference; + typedef const T& const_lambda_reference; /** * The main ALP scalar constructor. @@ -167,12 +170,12 @@ namespace alp { /** \internal No implementation notes. */ lambda_reference operator*() noexcept { - assert( getInitialized( *this ) ); + assert( internal::getInitialized( *this ) ); return value; } /** \internal No implementation notes. */ - const lambda_reference operator*() const noexcept { + const_lambda_reference operator*() const noexcept { assert( getInitialized( *this ) ); return value; } @@ -185,13 +188,13 @@ namespace alp { namespace internal { template< typename T, typename Structure > - bool getInitialized( Scalar< T, Structure, reference > &s ) noexcept { + bool getInitialized( const Scalar< T, Structure, reference > &s ) noexcept { return s.initialized; } template< typename T, typename Structure > void setInitialized( Scalar< T, Structure, reference > &s, bool initialized ) noexcept { - s.initialized = s; + s.initialized = initialized; } } // end namespace ``alp::internal'' diff --git a/include/alp/reference/vector.hpp b/include/alp/reference/vector.hpp index 28c414728..e29a72b31 100644 --- a/include/alp/reference/vector.hpp +++ b/include/alp/reference/vector.hpp @@ -499,7 +499,7 @@ namespace alp { internal::requires_allocation< View >::value > * = nullptr > - Vector( bool initialized, const size_t length, LambdaType lambda ) : + Vector( std::function< bool() > initialized, const size_t length, LambdaType lambda ) : base_type( initialized, length, 1, lambda ) {} /** diff --git a/include/alp/structures.hpp b/include/alp/structures.hpp index f6b2b2ff5..c775f7091 100644 --- a/include/alp/structures.hpp +++ b/include/alp/structures.hpp @@ -112,14 +112,14 @@ namespace alp { namespace internal { /** * @internal Compile-time check if a tuple of intervals is sorted and non-overlapping. - * E.g., a pair ( [a, b) [c, d) ) with a < b <= c < d + * E.g., a pair ( [a, b) [c, d) ) with a < b < c < d */ template< typename IntervalTuple > struct is_tuple_sorted_non_overlapping; template< std::ptrdiff_t _left0, std::ptrdiff_t _right0, std::ptrdiff_t _left1, std::ptrdiff_t _right1, typename... Intervals > struct is_tuple_sorted_non_overlapping < std::tuple< Interval< _left0, _right0 >, Interval< _left1, _right1 >, Intervals... > > { - static constexpr bool value = ( _right0 <= _left1 ) && is_tuple_sorted_non_overlapping< std::tuple< Interval< _left1, _right1 >, Intervals... > >::value; + static constexpr bool value = ( _right0 < _left1 ) && is_tuple_sorted_non_overlapping< std::tuple< Interval< _left1, _right1 >, Intervals... > >::value; }; template< std::ptrdiff_t _left, std::ptrdiff_t _right > @@ -179,8 +179,48 @@ namespace alp { template< typename Structure, typename... Structures > struct is_in< Structure, std::tuple< Structure, Structures... > > : std::true_type {}; - namespace internal { + + + /** + * @brief interval_ge + * True if the intervals in the tuple on the left are larger or equal than + * the ones in the tuple on the right. Notice that a single interval on the + * left may include multiple ones on the right, e.g., + * tuple< [1, 10) > ">=" tuple< [1, 3), [4, 5) > + */ + template < typename LeftTuple, typename RightTuple > + struct interval_ge { + static constexpr bool value = false; + }; + + template < > + struct interval_ge< std::tuple< >, std::tuple< > > { + + static constexpr bool value = true; + + }; + + template < typename IntervalL, typename... IntervalsL > + struct interval_ge< std::tuple< IntervalL, IntervalsL... >, std::tuple< > > { + + static constexpr bool value = true; + + }; + + template < typename IntervalL, typename... IntervalsL, typename IntervalR, typename... IntervalsR > + struct interval_ge< std::tuple< IntervalL, IntervalsL... >, std::tuple< IntervalR, IntervalsR... > > { + + static constexpr bool value = + ( IntervalR::left >= IntervalL::left + && IntervalR::right <= IntervalL::right + && interval_ge< std::tuple, std::tuple< IntervalsR... > >::value ) + || ( IntervalR::left > IntervalL::right + && interval_ge< std::tuple< IntervalsL... >, std::tuple< IntervalR, IntervalsR... > >::value ); + + }; + + /** * @internal WIP interface. Symmetry may be extended so to describe the * direction of the symmetry. @@ -219,6 +259,19 @@ namespace alp { }; } // namespace internal + /** + * @brief band_ge + * If the band geometry of the left structure is larger or equal than + * the one of the right structure. + */ + + template < typename LeftStructure, typename RightStructure > + struct band_ge { + static constexpr bool value = internal::interval_ge< + typename LeftStructure::band_intervals, typename RightStructure::band_intervals + >::value; + }; + struct BaseStructure {}; struct UpperTriangular; @@ -329,7 +382,9 @@ namespace alp { * * @tparam Intervals One or more \a alp::Interval types specifying the * bands of the structure. These intervals should be - * non-overlapping and sorted according to the above + * non-overlapping, compact (at least distance of one + * band between two intervals), and sorted according + * to the above * assumption that all intervals are defined assuming * the main diagonal has position zero. * \a alp::LeftOpenInterval ( \a alp::RightOpenInterval) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index eceae4497..f2c0d33f7 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -221,6 +221,10 @@ add_grb_executables( dense_dot_norm2 dense_dot_norm2.cpp BACKENDS alp_reference ) +add_grb_executables( dense_fold dense_fold.cpp + BACKENDS alp_reference +) + add_grb_executables( dense_matrix_imf dense_matrix_imf.cpp BACKENDS alp_reference ) diff --git a/tests/unit/dense_structured_matrix.cpp b/tests/unit/dense_structured_matrix.cpp index 37cd887dd..02b697e55 100644 --- a/tests/unit/dense_structured_matrix.cpp +++ b/tests/unit/dense_structured_matrix.cpp @@ -21,9 +21,32 @@ #include #include #include +#include #include +template< typename Tp, std::size_t I > +typename std::enable_if< I == std::tuple_size< Tp >::value, void >::type +interval_print() { +} + +template< typename Tp, std::size_t I > +typename std::enable_if< I < std::tuple_size< Tp >::value, void >::type +interval_print() { + std::cout + << " <"<< std::tuple_element::type::left << ", " + << std::tuple_element::type::right << "> "; + interval_print(); +} + +template +void interval_print( std::string header ) { + std::cout << header; + interval_print(); + std::cout << "\n"; +} + + template< typename StructuredMat > void ask_questions( const StructuredMat & M, std::string name ) { @@ -73,6 +96,57 @@ void alp_program( const size_t & n, alp::RC & rc ) { //std::cout << "v_view2( " << alp::getLength( v_view2 ) << " )" << std::endl; // TODO: temporarily comented until containers are ready + typedef alp::structures::Band< alp::LeftOpenInterval< -2 >, alp::Interval<1, 4> > LB0; + typedef alp::structures::Band< alp::Interval< -5, -4 >, alp::Interval< -3, -2 > > RB0; + + typedef std::tuple< alp::LeftOpenInterval< -2 >, alp::Interval<1, 4> > LeftTuple1; + typedef std::tuple< alp::Interval< -5, -4 >, alp::Interval< -3, -2 >, alp::Interval< 2 > > RightTuple1; + + typedef std::tuple< alp::LeftOpenInterval< -2 >, alp::Interval<1, 4> > LeftTuple2; + typedef std::tuple< alp::Interval< -5, -4 >, alp::Interval< -3, -2 >, alp::Interval< 2 >, alp::Interval< 4 > > RightTuple2; + + typedef std::tuple< alp::OpenInterval > LeftTuple3; + typedef std::tuple< alp::Interval< -5, -4 >, alp::Interval< -3, -2 >, alp::Interval< 2 >, alp::Interval< 4 > > RightTuple3; + + typedef std::tuple< > LeftTuple4; + typedef std::tuple< alp::Interval< -5, -4 >, alp::Interval< -3, -2 >, alp::Interval< 2 >, alp::Interval< 4 > > RightTuple4; + + typedef std::tuple< > Tuple5; + + std::cout << "\n"; + interval_print( "LB0= " ); + interval_print( "RB0= " ); + std::cout << "Is LB0 >= RB0: " << alp::structures::band_ge< LB0, RB0 >::value << std::endl; + std::cout << "Is RB0 >= LB0: " << alp::structures::band_ge< RB0, LB0 >::value << std::endl; + + std::cout << "\n"; + interval_print( "LeftTuple1= " ); + interval_print( "RightTuple1= " ); + std::cout << "Is super set 1: " << alp::structures::internal::interval_ge< LeftTuple1, RightTuple1 >::value << std::endl; + std::cout << "Is super set 1 rev: " << alp::structures::internal::interval_ge< RightTuple1, LeftTuple1 >::value << std::endl; + + std::cout << "\n"; + interval_print( "LeftTuple2= " ); + interval_print( "RightTuple2= " ); + std::cout << "Is super set 2: " << alp::structures::internal::interval_ge< LeftTuple2, RightTuple2 >::value << std::endl; + std::cout << "Is super set 2 rev: " << alp::structures::internal::interval_ge< RightTuple2, LeftTuple2 >::value << std::endl; + + std::cout << "\n"; + interval_print( "LeftTuple3= " ); + interval_print( "RightTuple3= " ); + std::cout << "Is super set 3: " << alp::structures::internal::interval_ge< LeftTuple3, RightTuple3 >::value << std::endl; + std::cout << "Is super set 3 rev: " << alp::structures::internal::interval_ge< RightTuple3, LeftTuple3 >::value << std::endl; + + std::cout << "\n"; + interval_print( "LeftTuple4= " ); + interval_print( "RightTuple4= " ); + std::cout << "Is super set 4: " << alp::structures::internal::interval_ge< LeftTuple4, RightTuple4 >::value << std::endl; + std::cout << "Is super set 4 rev: " << alp::structures::internal::interval_ge< RightTuple4, LeftTuple4 >::value << std::endl; + + std::cout << "\n"; + interval_print( "Tuple5= " ); + std::cout << "Is super set 5: " << alp::structures::internal::interval_ge< Tuple5, Tuple5 >::value << std::endl; + //alp::Matrix< float, alp::structures::Band< alp::Interval<-2, 5> > > BM0( n, n ); //alp::Matrix< float, alp::structures::Band< alp::RightOpenInterval<-2> > > BM1( n, n ); //alp::Matrix< float, alp::structures::Band< alp::LeftOpenInterval<-2> > > BM2( n, n );