Fixed vector assignment with operations: +=, -=, *=, /= (cd12eadf) · Commits · TNL / tnl-dev

src/TNL/Containers/Algorithms/VectorAssignment.h

+90 −122

Original line number	Diff line number	Diff line
		@@ -27,36 +27,13 @@ template< typename Vector,
		struct VectorAssignment{};

		/**
		* \brief Vector addition
		* \brief Vector assignment with an operation: +=, -=, *=, /=
		*/
		template< typename Vector,
		typename T,
		bool hasSubscriptOperator = HasSubscriptOperator< T >::value >
		struct VectorAddition{};

		/**
		* \brief Vector subtraction
		*/
		template< typename Vector,
		typename T,
		bool hasSubscriptOperator = HasSubscriptOperator< T >::value >
		struct VectorSubtraction{};

		/**
		* \brief Vector multiplication
		*/
		template< typename Vector,
		typename T,
		bool hasSubscriptOperator = HasSubscriptOperator< T >::value >
		struct VectorMultiplication{};

		/**
		* \brief Vector division
		*/
		template< typename Vector,
		typename T,
		bool hasSubscriptOperator = HasSubscriptOperator< T >::value >
		struct VectorDivision{};
		bool hasSubscriptOperator = HasSubscriptOperator< T >::value,
		bool hasSetSizeMethod = HasSetSizeMethod< T >::value >
		struct VectorAssignmentWithOperation{};

		/**
		* \brief Specialization of ASSIGNEMENT with subscript operator
		@@ -80,6 +57,8 @@ struct VectorAssignment< Vector, T, true >

		static void assign( Vector& v, const T& t )
		{
		static_assert( std::is_same< typename Vector::DeviceType, typename T::DeviceType >::value,
		"Cannot assign an expression to a vector allocated on a different device." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		@@ -130,54 +109,58 @@ struct VectorAssignment< Vector, T, false >
		};

		/**
		* \brief Specialization of ADDITION with subscript operator
		* \brief Specialization for types with subscript operator and setSize
		* method - i.e. for Vectors.
		*
		* This is necessary, because Vectors cannot be passed-by-value to CUDA kernels
		* so we have to do it with views.
		*/
		template< typename Vector,
		typename T >
		struct VectorAddition< Vector, T, true >
		struct VectorAssignmentWithOperation< Vector, T, true, true >
		{
		__cuda_callable__
		static void additionStatic( Vector& v, const T& t )
		static void addition( Vector& v, const T& t )
		{
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] += t[ i ];
		VectorAssignmentWithOperation< Vector, typename Vector::ConstViewType >::addition( v, t.getConstView() );
		}

		static void addition( Vector& v, const T& t )
		static void subtraction( Vector& v, const T& t )
		{
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;
		VectorAssignmentWithOperation< Vector, typename Vector::ConstViewType >::subtraction( v, t.getConstView() );
		}

		RealType* data = v.getData();
		auto add = [=] __cuda_callable__ ( IndexType i )
		static void multiplication( Vector& v, const T& t )
		{
		data[ i ] += t[ i ];
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), add );
		VectorAssignmentWithOperation< Vector, typename Vector::ConstViewType >::multiplication( v, t.getConstView() );
		}

		static void division( Vector& v, const T& t )
		{
		VectorAssignmentWithOperation< Vector, typename Vector::ConstViewType >::subtraction( v, t.getConstView() );
		}
		};

		/**
		* \brief Specialization of ADDITION for array-value assignment for other types. We assume
		* that T is convertible to Vector::ValueType.
		* \brief Specialization for types with subscript operator, but without setSize
		* method - i.e. for expressions, views and static vectors.
		*/
		template< typename Vector,
		typename T >
		struct VectorAddition< Vector, T, false >
		struct VectorAssignmentWithOperation< Vector, T, true, false >
		{
		__cuda_callable__
		static void additionStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_GT( v.getSize(), 0, "Cannot assign value to empty vector." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] += t;
		v[ i ] += t[ i ];
		}

		static void addition( Vector& v, const T& t )
		{
		static_assert( std::is_same< typename Vector::DeviceType, typename T::DeviceType >::value,
		"Cannot assign an expression to a vector allocated on a different device." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;
		@@ -185,19 +168,11 @@ struct VectorAddition< Vector, T, false >
		RealType* data = v.getData();
		auto add = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] += t;
		data[ i ] += t[ i ];
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), add );
		}
		};

		/**
		* \brief Specialization of SUBTRACTION with subscript operator
		*/
		template< typename Vector,
		typename T >
		struct VectorSubtraction< Vector, T, true >
		{
		__cuda_callable__
		static void subtractionStatic( Vector& v, const T& t )
		{
		@@ -208,6 +183,8 @@ struct VectorSubtraction< Vector, T, true >

		static void subtraction( Vector& v, const T& t )
		{
		static_assert( std::is_same< typename Vector::DeviceType, typename T::DeviceType >::value,
		"Cannot assign an expression to a vector allocated on a different device." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		@@ -220,141 +197,132 @@ struct VectorSubtraction< Vector, T, true >
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), subtract );
		}
		};

		/**
		* \brief Specialization of SUBTRACTION for array-value assignment for other types. We assume
		* that T is convertible to Vector::ValueType.
		*/
		template< typename Vector,
		typename T >
		struct VectorSubtraction< Vector, T, false >
		{
		__cuda_callable__
		static void subtractionStatic( Vector& v, const T& t )
		static void multiplicationStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_GT( v.getSize(), 0, "Cannot assign value to empty vector." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] -= t;
		v[ i ] *= t[ i ];
		}

		static void subtraction( Vector& v, const T& t )
		static void multiplication( Vector& v, const T& t )
		{
		static_assert( std::is_same< typename Vector::DeviceType, typename T::DeviceType >::value,
		"Cannot assign an expression to a vector allocated on a different device." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;

		RealType* data = v.getData();
		auto subtract = [=] __cuda_callable__ ( IndexType i )
		auto multiply = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] -= t;
		data[ i ] *= t[ i ];
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), subtract );
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), multiply );
		}
		};

		/**
		* \brief Specialization of MULTIPLICATION with subscript operator
		*/
		template< typename Vector,
		typename T >
		struct VectorMultiplication< Vector, T, true >
		{
		__cuda_callable__
		static void multiplicationStatic( Vector& v, const T& t )
		static void divisionStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] *= t[ i ];
		v[ i ] /= t[ i ];
		}

		static void multiplication( Vector& v, const T& t )
		static void division( Vector& v, const T& t )
		{
		static_assert( std::is_same< typename Vector::DeviceType, typename T::DeviceType >::value,
		"Cannot assign an expression to a vector allocated on a different device." );
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;

		RealType* data = v.getData();
		auto multiply = [=] __cuda_callable__ ( IndexType i )
		auto divide = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] *= t[ i ];
		data[ i ] /= t[ i ];
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), multiply );
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), divide );
		}
		};

		/**
		* \brief Specialization of MULTIPLICATION for array-value assignment for other types. We assume
		* \brief Specialization for array-value assignment for other types. We assume
		* that T is convertible to Vector::ValueType.
		*/
		template< typename Vector,
		typename T >
		struct VectorMultiplication< Vector, T, false >
		struct VectorAssignmentWithOperation< Vector, T, false, false >
		{
		__cuda_callable__
		static void multiplicationStatic( Vector& v, const T& t )
		static void additionStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_GT( v.getSize(), 0, "Cannot assign value to empty vector." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] *= t;
		v[ i ] += t;
		}

		static void multiplication( Vector& v, const T& t )
		static void addition( Vector& v, const T& t )
		{
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;

		RealType* data = v.getData();
		auto multiply = [=] __cuda_callable__ ( IndexType i )
		auto add = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] *= t;
		data[ i ] += t;
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), multiply );
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), add );
		}
		};

		__cuda_callable__
		static void subtractionStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_GT( v.getSize(), 0, "Cannot assign value to empty vector." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] -= t;
		}

		/**
		* \brief Specialization of DIVISION with subscript operator
		*/
		template< typename Vector,
		typename T >
		struct VectorDivision< Vector, T, true >
		static void subtraction( Vector& v, const T& t )
		{
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;

		RealType* data = v.getData();
		auto subtract = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] -= t;
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), subtract );
		}

		__cuda_callable__
		static void divisionStatic( Vector& v, const T& t )
		static void multiplicationStatic( Vector& v, const T& t )
		{
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		TNL_ASSERT_GT( v.getSize(), 0, "Cannot assign value to empty vector." );
		for( decltype( v.getSize() ) i = 0; i < v.getSize(); i++ )
		v[ i ] /= t[ i ];
		v[ i ] *= t;
		}

		static void division( Vector& v, const T& t )
		static void multiplication( Vector& v, const T& t )
		{
		TNL_ASSERT_EQ( v.getSize(), t.getSize(), "The sizes of the vectors must be equal." );
		using RealType = typename Vector::RealType;
		using DeviceType = typename Vector::DeviceType;
		using IndexType = typename Vector::IndexType;

		RealType* data = v.getData();
		auto divide = [=] __cuda_callable__ ( IndexType i )
		auto multiply = [=] __cuda_callable__ ( IndexType i )
		{
		data[ i ] /= t[ i ];
		data[ i ] *= t;
		};
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), divide );
		ParallelFor< DeviceType >::exec( ( IndexType ) 0, v.getSize(), multiply );
		}
		};

		/**
		* \brief Specialization of DIVISION for array-value assignment for other types. We assume
		* that T is convertible to Vector::ValueType.
		*/
		template< typename Vector,
		typename T >
		struct VectorDivision< Vector, T, false >
		{
		__cuda_callable__
		static void divisionStatic( Vector& v, const T& t )
		{

src/TNL/Containers/Vector.hpp

+4 −4

Original line number	Diff line number	Diff line
		@@ -160,7 +160,7 @@ Vector< Real, Device, Index, Allocator >&
		Vector< Real, Device, Index, Allocator >::
		operator-=( const VectorExpression& expression )
		{
		Algorithms::VectorSubtraction< Vector, VectorExpression >::subtraction( *this, expression );
		Algorithms::VectorAssignmentWithOperation< Vector, VectorExpression >::subtraction( *this, expression );
		return *this;
		}

		@@ -173,7 +173,7 @@ Vector< Real, Device, Index, Allocator >&
		Vector< Real, Device, Index, Allocator >::
		operator+=( const VectorExpression& expression )
		{
		Algorithms::VectorAddition< Vector, VectorExpression >::addition( *this, expression );
		Algorithms::VectorAssignmentWithOperation< Vector, VectorExpression >::addition( *this, expression );
		return *this;
		}

		@@ -186,7 +186,7 @@ Vector< Real, Device, Index, Allocator >&
		Vector< Real, Device, Index, Allocator >::
		operator*=( const VectorExpression& expression )
		{
		Algorithms::VectorMultiplication< Vector, VectorExpression >::multiplication( *this, expression );
		Algorithms::VectorAssignmentWithOperation< Vector, VectorExpression >::multiplication( *this, expression );
		return *this;
		}

		@@ -199,7 +199,7 @@ Vector< Real, Device, Index, Allocator >&
		Vector< Real, Device, Index, Allocator >::
		operator/=( const VectorExpression& expression )
		{
		Algorithms::VectorDivision< Vector, VectorExpression >::division( *this, expression );
		Algorithms::VectorAssignmentWithOperation< Vector, VectorExpression >::division( *this, expression );
		return *this;
		}

src/TNL/Containers/VectorView.hpp

+4 −4

Original line number	Diff line number	Diff line
		@@ -110,7 +110,7 @@ VectorView< Real, Device, Index >&
		VectorView< Real, Device, Index >::
		operator-=( const VectorExpression& expression )
		{
		Algorithms::VectorSubtraction< VectorView, VectorExpression >::subtraction( *this, expression );
		Algorithms::VectorAssignmentWithOperation< VectorView, VectorExpression >::subtraction( *this, expression );
		return *this;
		}

		@@ -122,7 +122,7 @@ VectorView< Real, Device, Index >&
		VectorView< Real, Device, Index >::
		operator+=( const VectorExpression& expression )
		{
		Algorithms::VectorAddition< VectorView, VectorExpression >::addition( *this, expression );
		Algorithms::VectorAssignmentWithOperation< VectorView, VectorExpression >::addition( *this, expression );
		return *this;
		}

		@@ -134,7 +134,7 @@ VectorView< Real, Device, Index >&
		VectorView< Real, Device, Index >::
		operator*=( const VectorExpression& expression )
		{
		Algorithms::VectorMultiplication< VectorView, VectorExpression >::multiplication( *this, expression );
		Algorithms::VectorAssignmentWithOperation< VectorView, VectorExpression >::multiplication( *this, expression );
		return *this;
		}

		@@ -146,7 +146,7 @@ VectorView< Real, Device, Index >&
		VectorView< Real, Device, Index >::
		operator/=( const VectorExpression& expression )
		{
		Algorithms::VectorDivision< VectorView, VectorExpression >::division( *this, expression );
		Algorithms::VectorAssignmentWithOperation< VectorView, VectorExpression >::division( *this, expression );
		return *this;
		}

src/TNL/TypeTraits.h

+28 −0

Original line number	Diff line number	Diff line
		@@ -49,6 +49,34 @@ public:
		static constexpr bool value = ( sizeof( test< T >(0) ) == sizeof( YesType ) );
		};

		/**
		* \brief Type trait for checking if T has setSize method.
		*/
		template< typename T >
		class HasSetSizeMethod
		{
		private:
		template< typename U >
		static constexpr auto check(U*)
		-> typename
		std::enable_if_t<
		std::is_same<
		decltype( std::declval<U>().setSize(0) ),
		void
		>::value,
		std::true_type
		>;

		template< typename >
		static constexpr std::false_type check(...);

		using type = decltype(check<T>(0));

		public:
		static constexpr bool value = type::value;
		};


		/**
		* \brief Type trait for checking if T has operator[] taking one index argument.
		*/

src/UnitTests/Containers/VectorBinaryOperationsTest.h

+118 −0

Original line number	Diff line number	Diff line
		@@ -272,6 +272,124 @@ TYPED_TEST( VectorBinaryOperationsTest, division )
		EXPECT_EQ( (L1 + L1) / (L1 + L1), L1 );
		}

		template< typename Left, typename Right, std::enable_if_t< std::is_const<typename Left::RealType>::value, bool > = true >
		void test_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{}
		template< typename Left, typename Right, std::enable_if_t< ! std::is_const<typename Left::RealType>::value, bool > = true >
		void test_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{
		// with vector or vector view
		L1 = R2;
		EXPECT_EQ( L1, R2 );
		// with scalar
		L1 = 1;
		EXPECT_EQ( L1, 1 );
		// with expression
		L1 = R1 + R1;
		EXPECT_EQ( L1, R1 + R1 );
		}
		TYPED_TEST( VectorBinaryOperationsTest, assignment )
		{
		SETUP_BINARY_VECTOR_TEST( VECTOR_TEST_SIZE );
		test_assignment( L1, L2, R1, R2 );
		}

		template< typename Left, typename Right, std::enable_if_t< std::is_const<typename Left::RealType>::value, bool > = true >
		void test_add_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{}
		template< typename Left, typename Right, std::enable_if_t< ! std::is_const<typename Left::RealType>::value, bool > = true >
		void test_add_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{
		// with vector or vector view
		L1 += R2;
		EXPECT_EQ( L1, R1 + R2 );
		// with scalar
		L1 = 1;
		L1 += 2;
		EXPECT_EQ( L1, 3 );
		// with expression
		L1 = 1;
		L1 += R1 + R1;
		EXPECT_EQ( L1, R1 + R1 + R1 );
		}
		TYPED_TEST( VectorBinaryOperationsTest, add_assignment )
		{
		SETUP_BINARY_VECTOR_TEST( VECTOR_TEST_SIZE );
		test_add_assignment( L1, L2, R1, R2 );
		}

		template< typename Left, typename Right, std::enable_if_t< std::is_const<typename Left::RealType>::value, bool > = true >
		void test_subtract_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{}
		template< typename Left, typename Right, std::enable_if_t< ! std::is_const<typename Left::RealType>::value, bool > = true >
		void test_subtract_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{
		// with vector or vector view
		L1 -= R2;
		EXPECT_EQ( L1, R1 - R2 );
		// with scalar
		L1 = 1;
		L1 -= 2;
		EXPECT_EQ( L1, -1 );
		// with expression
		L1 = 1;
		L1 -= R1 + R1;
		EXPECT_EQ( L1, -R1 );
		}
		TYPED_TEST( VectorBinaryOperationsTest, subtract_assignment )
		{
		SETUP_BINARY_VECTOR_TEST( VECTOR_TEST_SIZE );
		test_subtract_assignment( L1, L2, R1, R2 );
		}

		template< typename Left, typename Right, std::enable_if_t< std::is_const<typename Left::RealType>::value, bool > = true >
		void test_multiply_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{}
		template< typename Left, typename Right, std::enable_if_t< ! std::is_const<typename Left::RealType>::value, bool > = true >
		void test_multiply_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{
		// with vector or vector view
		L1 *= R2;
		EXPECT_EQ( L1, R2 );
		// with scalar
		L1 = 1;
		L1 *= 2;
		EXPECT_EQ( L1, 2 );
		// with expression
		L1 = 1;
		L1 *= R1 + R1;
		EXPECT_EQ( L1, R1 + R1 );
		}
		TYPED_TEST( VectorBinaryOperationsTest, multiply_assignment )
		{
		SETUP_BINARY_VECTOR_TEST( VECTOR_TEST_SIZE );
		test_multiply_assignment( L1, L2, R1, R2 );
		}

		template< typename Left, typename Right, std::enable_if_t< std::is_const<typename Left::RealType>::value, bool > = true >
		void test_divide_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{}
		template< typename Left, typename Right, std::enable_if_t< ! std::is_const<typename Left::RealType>::value, bool > = true >
		void test_divide_assignment( Left& L1, Left& L2, Right& R1, Right& R2 )
		{
		// with vector or vector view
		L2 /= R2;
		EXPECT_EQ( L1, R1 );
		// with scalar
		L2 = 2;
		L2 /= 2;
		EXPECT_EQ( L1, 1 );
		// with expression
		L2 = 2;
		L2 /= R1 + R1;
		EXPECT_EQ( L1, R1 );
		}
		TYPED_TEST( VectorBinaryOperationsTest, divide_assignment )
		{
		SETUP_BINARY_VECTOR_TEST( VECTOR_TEST_SIZE );
		test_divide_assignment( L1, L2, R1, R2 );
		}

		TYPED_TEST( VectorBinaryOperationsTest, scalarProduct )
		{
		// this test expects an odd size