510 foldl foldr dense reference implementation

include/alp/reference/blas1.hpp

@@ -23,11 +23,17 @@
 #ifndef _H_ALP_REFERENCE_BLAS1
 #define _H_ALP_REFERENCE_BLAS1
 
+#include <functional>
+
 #include <alp/backends.hpp>
 #include <alp/config.hpp>
 #include <alp/rc.hpp>
 #include <alp/scalar.hpp>
 #include <alp/density.hpp>
+#include <alp/matrix.hpp>
+#include <alp/vector.hpp>
+#include <alp/blas0.hpp>
+
 
 #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;
 	}
 

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''

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''
 

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 ) {}
 
 			/**

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<IntervalL, IntervalsL... >, 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) 

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
 )