516 implement eWiseLambda

include/alp/reference/blas1.hpp

@@ -27,12 +27,12 @@
 #include <alp/backends.hpp>
 #include <alp/config.hpp>
 #include <alp/rc.hpp>
-#include <alp/scalar.hpp>
-#include <alp/vector.hpp>
 #include <alp/density.hpp>
-#include <alp/matrix.hpp>
-#include <alp/vector.hpp>
-#include <alp/blas0.hpp>
+#include "scalar.hpp"
+#include "matrix.hpp"
+#include "vector.hpp"
+#include "blas0.hpp"
+#include "blas2.hpp"
 
 #ifndef NO_CAST_ASSERT
 #define NO_CAST_ASSERT( x, y, z )                                              \
@@ -3952,16 +3952,20 @@ namespace alp {
 	 * @see Vector::operator[]()
 	 * @see Vector::lambda_reference
 	 */
-	template< typename Func,
+	template<
+		typename Func,
 		typename DataType1, typename DataStructure1, typename DataView1, typename InputImfR1, typename InputImfC1,
 		typename DataType2, typename DataStructure2, typename DataView2, typename InputImfR2, typename InputImfC2,
-		typename... Args >
-	RC eWiseLambda( const Func f,
-		const Vector< DataType1, DataStructure1, Density::Dense, DataView1, InputImfR1, InputImfC1, reference > & x,
-		const Vector< DataType2, DataStructure2, Density::Dense, DataView2, InputImfR2, InputImfC2, reference > & y,
-		Args const &... args ) {
+		typename... Args
+	>
+	RC eWiseLambda(
+		const Func f,
+		Vector< DataType1, DataStructure1, Density::Dense, DataView1, InputImfR1, InputImfC1, reference > &x,
+		const Vector< DataType2, DataStructure2, Density::Dense, DataView2, InputImfR2, InputImfC2, reference > &y,
+		Args const &... args
+	) {
 		// catch mismatch
-		if( size( x ) != size( y ) ) {
+		if( getLength( x ) != getLength( y ) ) {
 			return MISMATCH;
 		}
 		// continue
@@ -3975,15 +3979,22 @@ namespace alp {
 	 * @see Vector::operator[]()
 	 * @see Vector::lambda_reference
 	 */
-	template< typename Func,
+	template<
+		typename Func,
 		typename DataType, typename DataStructure, typename DataView, typename DataImfR, typename DataImfC
 	>
-	RC eWiseLambda( const Func f, const Vector< DataType, DataStructure, Density::Dense, DataView, DataImfR, DataImfC, reference > & x ) {
-	#ifdef _DEBUG
+	RC eWiseLambda( const Func f, Vector< DataType, DataStructure, Density::Dense, DataView, DataImfR, DataImfC, reference > &x ) {
+#ifdef _DEBUG
 		std::cout << "Info: entering eWiseLambda function on vectors.\n";
-	#endif
-		throw std::runtime_error( "Needs an implementation." );
-		return SUCCESS;
+#endif
+		auto x_as_matrix = get_view< view::matrix >( x );
+		return eWiseLambda(
+			[ &f ]( const size_t i, const size_t j, DataType &val ) {
+				(void)j;
+				f( i, val );
+			},
+			x_as_matrix
+		);
 	}
 
 	/**

include/alp/reference/blas2.hpp

@@ -28,6 +28,7 @@
 #include <alp/backends.hpp>
 #include <alp/config.hpp>
 #include <alp/rc.hpp>
+#include <alp/matrix.hpp>
 
 namespace alp {
 
@@ -330,23 +331,111 @@ namespace alp {
 		return SUCCESS;
 	}
 
+	namespace internal {
+
+		/**
+		 * Applies the provided function to each element of the given band.
+		 * This function is called by the public eWiseLambda variant.
+		 * Forward declaration. Specializations handle bound checking.
+		 */
+		template<
+			size_t BandIndex, typename Func,
+			typename DataType, typename Structure, typename View, typename ImfR, typename ImfC,
+			typename std::enable_if_t<
+				BandIndex >= std::tuple_size< typename Structure::band_intervals >::value
+			> * = nullptr
+		>
+		RC eWiseLambda(
+			const Func f,
+			alp::Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A
+		);
+
+		/** Specialization for an out-of-bounds band index */
+		template<
+			size_t BandIndex, typename Func,
+			typename DataType, typename Structure, typename View, typename ImfR, typename ImfC,
+			typename std::enable_if_t<
+				BandIndex >= std::tuple_size< typename Structure::band_intervals >::value
+			> * = nullptr
+		>
+		RC eWiseLambda(
+			const Func f,
+			alp::Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A
+		) {
+			(void)f;
+			(void)A;
+			// nothing to do
+			return SUCCESS;
+		}
+
+		/**
+		 * Specialization for a within-the-range band index.
+		 * Applies the provided function to each element of the given band.
+		 * Upon completion, calls itself for the next band.
+		 */
+		template<
+			size_t BandIndex, typename Func,
+			typename DataType, typename Structure, typename View, typename ImfR, typename ImfC,
+			typename std::enable_if_t<
+				BandIndex < std::tuple_size< typename Structure::band_intervals >::value
+			> * = nullptr
+		>
+		RC eWiseLambda(
+			const Func f,
+			alp::Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A
+		) {
+			const size_t M = nrows( A );
+			const size_t N = ncols( A );
+
+			const std::ptrdiff_t l = structures::get_lower_bandwidth< BandIndex >( A );
+			const std::ptrdiff_t u = structures::get_upper_bandwidth< BandIndex >( A );
+
+			// In case of symmetry the iteration domain intersects the the upper
+			// (or lower) domain of A
+			constexpr bool is_sym_a = structures::is_a< Structure, structures::Symmetric >::value;
+
+			// Temporary until adding multiple symmetry directions
+			constexpr bool sym_up_a = is_sym_a;
+
+			/** i-coordinate lower and upper limits considering matrix size and band limits */
+			const std::ptrdiff_t i_l_lim = std::max( static_cast< std::ptrdiff_t >( 0 ), -u );
+			const std::ptrdiff_t i_u_lim = std::min( M, -l + N );
+
+			for( size_t i = static_cast< size_t >( i_l_lim ); i < static_cast< size_t >( i_u_lim ); ++i ) {
+				/** j-coordinate lower and upper limits considering matrix size and symmetry */
+				const std::ptrdiff_t j_sym_l_lim = is_sym_a && sym_up_a ? i : 0;
+				const std::ptrdiff_t j_sym_u_lim = is_sym_a && !sym_up_a ? i + 1 : N;
+				/** j-coordinate lower and upper limits, also considering the band limits in addition to the factors above */
+				const std::ptrdiff_t j_l_lim = std::max( j_sym_l_lim, l );
+				const std::ptrdiff_t j_u_lim = std::min( j_sym_u_lim, u );
+
+				for( size_t j = static_cast< size_t >( j_l_lim ); j < static_cast< size_t >( j_u_lim ); ++j ) {
+					auto &a_val = internal::access( A, internal::getStorageIndex( A, i, j ) );
+					f( i, j, a_val );
+				}
+			}
+			return eWiseLambda< BandIndex + 1 >( f, A );
+		}
+
+	} // namespace internal
+
 	/**
-	 * Straightforward implementation using the column-major layout.
+	 * Delegates to single-band variant.
 	 *
 	 * @see alp::eWiseLambda for the user-level specification.
 	 */
-	template< class ActiveDistribution, typename Func,
+	template<
+		typename Func,
 		typename DataType, typename Structure, typename View, typename ImfR, typename ImfC
 	>
-	RC eWiseLambda( const Func f,
-		const Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > & A,
-		const size_t s,
-		const size_t P ) {
+	RC eWiseLambda(
+		const Func f,
+		Matrix< DataType, Structure, Density::Dense, View, ImfR, ImfC, reference > &A
+	) {
 #ifdef _DEBUG
 		std::cout << "entering alp::eWiseLambda (matrices, reference ). A is " << alp::nrows( A ) << " by " << alp::ncols( A ) << " and holds " << alp::nnz( A ) << " nonzeroes.\n";
 #endif
-		throw std::runtime_error( "Needs an implementation." );
-		return SUCCESS;
+		return internal::eWiseLambda< 0 >( f, A );
 	}
 
 	/**
@@ -355,21 +444,25 @@ namespace alp {
 	 *
 	 * @see alp::eWiseLambda for the user-level specification.
 	 */
-	template< typename Func,
+	template<
+		typename Func,
 		typename DataType1, typename DataStructure1, typename DataView1, typename DataImfR1, typename DataImfC1,
 		typename DataType2, typename DataStructure2, typename DataView2, typename DataImfR2, typename DataImfC2,
 		typename... Args
 	>
-	RC eWiseLambda( const Func f,
-		const Matrix< DataType1, DataStructure1, Density::Dense, DataView1, DataImfR1, DataImfC1, reference > & A,
-		const Vector< DataType2, DataStructure2, Density::Dense, DataView2, DataImfR2, DataImfC2, reference > x,
-		Args... args ) {
+	RC eWiseLambda(
+		const Func f,
+		Matrix< DataType1, DataStructure1, Density::Dense, DataView1, DataImfR1, DataImfC1, reference > &A,
+		const Vector< DataType2, DataStructure2, Density::Dense, DataView2, DataImfR2, DataImfC2, reference > &x,
+		Args const &... args
+	) {
 		// do size checking
-		if( ! ( size( x ) == nrows( A ) || size( x ) == ncols( A ) ) ) {
-			std::cerr << "Mismatching dimensions: given vector of size " << size( x ) << " has nothing to do with either matrix dimension (" << nrows( A ) << " nor " << ncols( A ) << ").\n";
+		if( !( getLength( x ) == nrows( A ) || getLength( x ) == ncols( A ) ) ) {
+			std::cerr << "Mismatching dimensions: given vector of size " << size( x )
+				<< " has nothing to do with either matrix dimension (" << nrows( A ) << " nor " << ncols( A ) << ").\n";
 			return MISMATCH;
 		}
-		// no need for synchronisation, everything is local in reference implementation
+
 		return eWiseLambda( f, A, args... );
 	}
 

tests/unit/CMakeLists.txt

@@ -217,6 +217,10 @@ add_grb_executables( dense_matrix_access dense_matrix_access.cpp
 #	BACKENDS alp_reference
 #)
 
+add_grb_executables( dense_eWiseLambda dense_eWiseLambda.cpp
+	BACKENDS alp_reference
+)
+
 add_grb_executables( dense_dot_norm2 dense_dot_norm2.cpp
 	BACKENDS alp_reference
 )

tests/unit/dense_eWiseLambda.cpp

@@ -0,0 +1,151 @@
+
+/*
+ *   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 <iostream>
+#include <sstream>
+
+#include <alp.hpp>
+
+using namespace alp;
+
+static const int data1[ 15 ] = { 4, 7, 4, 6, 4, 7, 1, 7, 3, 6, 7, 5, 1, 8, 7 };
+static const size_t I[ 15 ] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 9, 9, 8, 7, 6 };
+static const size_t J[ 15 ] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 5, 7, 5, 1 };
+static const double data2[ 6 ] = { 1, 1, 1, 1, 1, 1 };
+static const size_t I2[ 6 ] = { 0, 1, 0, 2, 1, 2 };
+static const size_t J2[ 6 ] = { 1, 0, 2, 0, 2, 1 };
+static const double testv[ 3 ] = { 0.1, 2.1, -2.3 };
+
+void alp_program( const size_t &n, alp::RC &rc ) {
+	// initialize test
+	typedef int T;
+	alp::Matrix< T, structures::General > A( n, n );
+	alp::Vector< T > u( n );
+	alp::Vector< T > v( n );
+
+	internal::setInitialized( A, true );
+	internal::setInitialized( u, true );
+	internal::setInitialized( v, true );
+
+	// test eWiseLambda on matrix
+	rc = alp::eWiseLambda(
+		[]( const size_t i, const size_t j, T &val ) {
+			(void)i;
+			(void)j;
+			val = 1;
+		},
+		A
+	);
+	if( rc != alp::SUCCESS ){
+		std::cerr << "\talp::eWiseLambda (matrix, no vectors) FAILED\n";
+		return;
+	}
+	assert( internal::access( A, internal::getStorageIndex( A, 0, 0 ) ) == 1 );
+
+	// test eWiseLambda on vector
+	rc = alp::eWiseLambda(
+		[]( const size_t i, T &val ) {
+			(void)i;
+			val = 2;
+		},
+		v
+	);
+
+	if( rc != SUCCESS ) {
+		std::cerr << "\talp::eWiseLambda (vector) FAILED\n";
+		return;
+	}
+	assert( v[ 0 ] == 2 );
+
+	// test eWiseLambda on vector, consuming from another vector
+	rc = alp::eWiseLambda(
+		[ &v ]( const size_t i, T &val ) {
+			val = v[ i ];
+		},
+		u, v
+	);
+
+	if( rc != SUCCESS ) {
+		std::cerr << "\talp::eWiseLambda (vector, vectors...) FAILED\n";
+		return;
+	}
+	assert( v[ 0 ] = u[ 0 ] );
+
+	// test eWiseLambda on Matrix, consuming two other vectors
+	rc = alp::eWiseLambda(
+		[ &u, &v ]( const size_t i, const size_t j, T &val ) {
+			val = val + u[ i ] * v[ j ];
+		},
+		A, u, v
+	);
+	assert( internal::access( A, internal::getStorageIndex( A, 0, 0 ) ) == 1 + u[ 0 ] + v[ 0 ] );
+
+	if( rc != SUCCESS ) {
+		std::cerr << "\talp::eWiseLambda (matrix, vectors...) FAILED\n";
+		return;
+	}
+
+}
+
+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 % 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< AUTOMATIC > launcher;
+	alp::RC out;
+	if( launcher.exec( &alp_program, in, out, true ) != SUCCESS ) {
+		std::cerr << "Launching test FAILED\n";
+		return 255;
+	}
+	if( out != SUCCESS ) {
+		std::cout << "Test FAILED (" << alp::toString( out ) << ")" << std::endl;
+		return out;
+	} else {
+		std::cout << "Test OK" << std::endl;
+		return 0;
+	}
+}