Thu Nov 26 15:41:02 UTC 2015

This is an automated email from the git hooks/post-receive script.

ghisvail-guest pushed a commit to branch master
in repository shark.

commit 1b43f7c5998c17b8d397a2f651e72cf4129fe98a
Author: Oswin <oswin.krause at di.ku.dk>
Date:   Mon Oct 26 08:36:08 2015 +0100

    made the matrix-vector prod a block expression which internally uses kernels::gemv. this required a lot of changes in vector_expression to be able to write stuff like 2*prod(A,v)+b efficiently. namely two new expressions vector_scalar_multiply which multiplies a vector by a scalar and vector_addition which adds two vectors. This is because scalar multiplication can be carried out by gemv efficiently so it makes sense to have plus_assign_to or assign_to take an additional scalar argume [...]
---
 include/shark/LinAlg/BLAS/assignment.hpp        |   2 +-
 include/shark/LinAlg/BLAS/matrix_expression.hpp | 151 +++++--------
 include/shark/LinAlg/BLAS/vector_expression.hpp | 276 ++++++++++++++++++++++--
 3 files changed, 322 insertions(+), 107 deletions(-)

diff --git a/include/shark/LinAlg/BLAS/assignment.hpp b/include/shark/LinAlg/BLAS/assignment.hpp
index 30ccf97..eb95297 100644
--- a/include/shark/LinAlg/BLAS/assignment.hpp
+++ b/include/shark/LinAlg/BLAS/assignment.hpp
@@ -45,7 +45,7 @@ namespace blas {
 ////////////////////////////////////////////////////////////////////////////////////
 	
 namespace detail{
-	template<class VecX, class VecV>
+	template<class VecX, class VecV, class T>
 	void assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
 		kernels::assign(x,v);
 	}
diff --git a/include/shark/LinAlg/BLAS/matrix_expression.hpp b/include/shark/LinAlg/BLAS/matrix_expression.hpp
index 88b1593..5061838 100644
--- a/include/shark/LinAlg/BLAS/matrix_expression.hpp
+++ b/include/shark/LinAlg/BLAS/matrix_expression.hpp
@@ -513,120 +513,85 @@ safe_div(
 	return matrix_binary<E1, E2, functor_type>(e1,e2, functor_type(defaultValue));
 }
 
-template<class E1, class E2, class F>
-class matrix_vector_binary1:
-	public vector_expression<matrix_vector_binary1<E1, E2, F> > {
-
+template<class MatA, class VecV>
+class matrix_vector_prod:
+	public vector_expression<matrix_vector_prod<MatA, VecV> > {
 public:
-	typedef E1 expression1_type;
-	typedef E2 expression2_type;
-	typedef F functor_type;
-	typedef typename E1::const_closure_type expression1_closure_type;
-	typedef typename E2::const_closure_type expression2_closure_type;
-private:
-	typedef matrix_vector_binary1<E1, E2, F> self_type;
+	typedef typename MatA::const_closure_type matrix_closure_type;
+	typedef typename VecV::const_closure_type vector_closure_type;
 public:
 	typedef typename promote_traits<
-		typename E1::size_type, typename E2::size_type
-	>::promote_type size_type;
-	typedef typename promote_traits<
-		typename E1::difference_type, typename E2::difference_type
-	>::promote_type difference_type;
-	typedef typename F::result_type value_type;
-	typedef value_type scalar_type;
-	typedef value_type const_reference;
-	typedef const_reference reference;
-	typedef value_type const* const_pointer;
-	typedef const_pointer pointer;
-
-	typedef typename E1::index_type index_type;
-	typedef typename E1::const_index_pointer const_index_pointer;
-	typedef typename index_pointer<E1>::type index_pointer;
-
-	typedef const self_type const_closure_type;
+		typename MatA::scalar_type,
+		typename VecV::scalar_type
+	>::promote_type scalar_type;
+	typedef scalar_type value_type;
+	typedef typename MatA::size_type size_type;
+	typedef typename MatA::difference_type difference_type;
+	typedef typename MatA::index_type index_type;
+	typedef typename MatA::index_pointer index_pointer;
+	typedef typename MatA::const_index_pointer const_index_pointer;
+
+	typedef matrix_vector_prod<MatA, VecV> const_closure_type;
 	typedef const_closure_type closure_type;
 	typedef unknown_storage_tag storage_category;
-	typedef elementwise_tag evaluation_category;
+	typedef blockwise_tag evaluation_category;
+
+
+	//FIXME: This workaround is required to be able to generate
+	// temporary vectors
+	typedef typename MatA::const_row_iterator const_iterator;
+	typedef const_iterator iterator;
 
 	// Construction and destruction
-	matrix_vector_binary1(const expression1_type &e1, const expression2_type &e2):
-		m_expression1(e1), m_expression2(e2) {}
+	matrix_vector_prod(
+		matrix_closure_type const& matrix,
+		vector_closure_type  const& vector
+	):m_matrix(matrix), m_vector(vector) {}
 
-	// Accessors
 	size_type size() const {
-		return m_expression1.size1();
+		return m_matrix.size1();
 	}
-
-	// Expression accessors
-	const expression1_closure_type &expression1() const {
-		return m_expression1;
+	
+	matrix_closure_type const& matrix() const {
+		return m_matrix;
 	}
-	const expression2_closure_type &expression2() const {
-		return m_expression2;
+	vector_closure_type const& vector() const {
+		return m_vector;
 	}
-
-	// Element access
-	const_reference operator()(index_type i) const {
-		return functor_type::apply(m_expression1, m_expression2, i);
+	
+	//computation kernels
+	template<class VecX>
+	void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		x().clear();
+		plus_assign_to(x,alpha);
 	}
-
-	// Iterator types
-private:
-	typedef typename E1::const_row_iterator const_subrow_iterator_type;
-public:
-	typedef indexed_iterator<const_closure_type> const_iterator;
-	typedef const_iterator iterator;
-
-	const_iterator begin() const {
-		return const_iterator(*this,0);
+	template<class VecX>
+	void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		kernels::gemv(eval_block(m_matrix), eval_block(m_vector), x, alpha);
 	}
-	const_iterator end() const {
-		return const_iterator(*this,size());
+	
+	template<class VecX>
+	void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		kernels::gemv(eval_block(m_matrix), eval_block(m_vector), x, -alpha);
 	}
+	
 private:
-	expression1_closure_type m_expression1;
-	expression2_closure_type m_expression2;
+	matrix_closure_type m_matrix;
+	vector_closure_type m_vector;
 };
 
-template<class E1, class E2>
-struct matrix_vector_binary_traits {
-private:
-	template<class T>
-	struct matrix_vector_prod{
-		typedef T value_type;
-		typedef T result_type;
-
-		template<class U,class V>
-		static result_type apply(
-			const matrix_expression<U> &e1,
-			const vector_expression<V> &e2,
-			std::size_t i
-		) {
-			return inner_prod(row(e1,i),e2);
-		}
-	};
-public:
-	typedef typename E1::value_type T1;
-	typedef typename E2::value_type T2;
-	typedef unknown_storage_tag storage_category;
-	typedef row_major orientation;
-	typedef typename promote_traits<T1, T2>::promote_type promote_type;
-	typedef matrix_vector_binary1<E1, E2, matrix_vector_prod<promote_type> > expression_type;
-	typedef expression_type result_type;
-};
 
-template<class M, class V>
-typename matrix_vector_binary_traits<M,V>::result_type
-prod(matrix_expression<M> const& matrix,vector_expression<V> const& vector){
-	typedef typename matrix_vector_binary_traits<M,V>::expression_type expression;
-	return expression(matrix(),vector());
+/// \brief computes the matrix-vector product x+=Av
+template<class MatA, class VecV>
+matrix_vector_prod<MatA,VecV> prod(matrix_expression<MatA> const& A,vector_expression<VecV> const& v) {
+	return matrix_vector_prod<MatA,VecV>(A(),v());
 }
 
-template<class V, class M>
-typename matrix_vector_binary_traits<matrix_transpose<M const>,V>::result_type
-prod(vector_expression<V> const& vector,matrix_expression<M> const& matrix){
-	typedef typename matrix_vector_binary_traits<matrix_transpose<M const>,V>::expression_type expression;
-	return expression(trans(matrix),vector());
+/// \brief computes the matrix-vector product x+=v^TA
+template<class MatA, class VecV>
+matrix_vector_prod<matrix_transpose<MatA>,VecV> prod(vector_expression<VecV> const& v,matrix_expression<MatA> const& A) {
+	//compute it using the identity (v^TA)^T= A^Tv
+	return matrix_vector_prod<matrix_transpose<MatA>,VecV>(trans(A),v());
 }
 
 //matrix-matrix prod
diff --git a/include/shark/LinAlg/BLAS/vector_expression.hpp b/include/shark/LinAlg/BLAS/vector_expression.hpp
index ad6f39a..f0c0b00 100644
--- a/include/shark/LinAlg/BLAS/vector_expression.hpp
+++ b/include/shark/LinAlg/BLAS/vector_expression.hpp
@@ -9,6 +9,114 @@
 namespace shark {
 namespace blas {
 	
+///\brief Implements multiplications of a vector by a scalar
+template<class E>
+class vector_scalar_multiply:
+	public vector_expression<vector_scalar_multiply <E> > {
+	typedef vector_scalar_multiply<E> self_type;
+public:
+	typedef typename E::const_closure_type expression_closure_type;
+	typedef typename E::size_type size_type;
+	typedef typename E::difference_type difference_type;
+	typedef typename E::value_type value_type;
+	typedef typename E::scalar_type scalar_type;
+	typedef value_type const_reference;
+	typedef value_type reference;
+	typedef value_type const * const_pointer;
+	typedef value_type const*  pointer;
+
+	typedef typename E::index_type index_type;
+	typedef typename E::const_index_pointer const_index_pointer;
+	typedef typename index_pointer<E>::type index_pointer;
+
+	typedef self_type const const_closure_type;
+	typedef self_type closure_type;
+	typedef unknown_storage_tag storage_category;
+	typedef typename E::evaluation_category evaluation_category;
+
+	// Construction and destruction
+	// May be used as mutable expression.
+	vector_scalar_multiply(vector_expression<E> const &e, scalar_type scalar):
+		m_expression(e()), m_scalar(scalar) {}
+
+	// Accessors
+	size_type size() const {
+		return m_expression.size();
+	}
+
+	// Expression accessors
+	expression_closure_type const &expression() const {
+		return m_expression;
+	}
+
+public:
+	// Element access
+	const_reference operator()(index_type i) const {
+		return m_scalar * m_expression(i);
+	}
+
+	const_reference operator[](index_type i) const {
+		return m_scalar * m_expression(i);
+	}
+	
+	
+	//computation kernels
+	template<class VecX>
+	void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		m_expression.assign_to(x,m_scalar*alpha);
+	}
+	template<class VecX>
+	void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		m_expression.plus_assign_to(x,m_scalar*alpha);
+	}
+	
+	template<class VecX>
+	void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		m_expression.minus_assign_to(x,m_scalar*alpha);
+	}
+
+	
+	//iterators
+	typedef transform_iterator<typename E::const_iterator,scalar_multiply1<value_type, scalar_type> > const_iterator;
+	typedef const_iterator iterator;
+	
+	const_iterator begin() const {
+		return const_iterator(m_expression.begin(),scalar_multiply1<value_type, scalar_type>(m_scalar));
+	}
+	const_iterator end() const {
+		return const_iterator(m_expression.end(),scalar_multiply1<value_type, scalar_type>(m_scalar));
+	}
+private:
+	expression_closure_type m_expression;
+	scalar_type m_scalar;
+};
+
+
+template<class T, class E>
+typename boost::enable_if<
+	boost::is_convertible<T, typename E::scalar_type >,
+        vector_scalar_multiply<E>
+>::type
+operator* (vector_expression<E> const& e, T scalar){
+	typedef typename E::scalar_type scalar_type;
+	return vector_scalar_multiply<E>(e, scalar_type(scalar));
+}
+template<class T, class E>
+typename boost::enable_if<
+	boost::is_convertible<T, typename E::scalar_type >,
+        vector_scalar_multiply<E>
+>::type
+operator* (T scalar, vector_expression<E> const& e){
+	typedef typename E::scalar_type scalar_type;
+	return vector_scalar_multiply<E>(e, scalar_type(scalar));//explicit cast prevents warning, alternative would be to template vector_scalar_multiply on T as well
+}
+
+template<class E>
+vector_scalar_multiply<E> operator-(vector_expression<E> const& e){
+	typedef typename E::scalar_type scalar_type;
+	return vector_scalar_multiply<E>(e, scalar_type(-1));//explicit cast prevents warning, alternative would be to template vector_scalar_multiply on T as well
+}
+	
 /// \brief Vector expression representing a constant valued vector.
 template<class T>
 class scalar_vector:public vector_expression<scalar_vector<T> > {
@@ -96,8 +204,6 @@ template<class E, class F>
 class vector_unary:
 	public vector_expression<vector_unary<E, F> > {
 	typedef vector_unary<E, F> self_type;
-	typedef E const expression_type;
-
 public:
 	typedef F functor_type;
 	typedef typename E::const_closure_type expression_closure_type;
@@ -133,6 +239,33 @@ public:
 	expression_closure_type const &expression() const {
 		return m_expression;
 	}
+	
+	//computation kernels
+	template<class VecX>
+	void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		//compute this by first assigning the result of the argument and then applying
+		//the function to every element
+		assign(x,m_expression);
+		typename VecX::iterator end=x().end();
+		for(typename VecX::iterator pos =x().begin(); pos != end; ++pos){
+			*pos= alpha * m_functor(*pos);
+		}
+	}
+	template<class VecX>
+	void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		//First assign result of this expression to a temporary and then perform plus_assignment to x
+		typename vector_temporary<self_type>::type temporary(size());
+		assign_to(temporary,alpha);
+		plus_assign_to(x,temporary);
+	}
+	
+	template<class VecX>
+	void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		//First assign result of this expression to a temporary and then perform minus_assignment to x
+		typename vector_temporary<self_type>::type temporary(size());
+		assign_to(temporary,alpha);
+		minus_assign_to(x,temporary);
+	}
 
 public:
 	// Element access
@@ -168,7 +301,6 @@ name(vector_expression<E> const& e){\
 	typedef F<typename E::value_type> functor_type;\
 	return vector_unary<E, functor_type>(e, functor_type());\
 }
-SHARK_UNARY_VECTOR_TRANSFORMATION(operator-, scalar_negate)
 SHARK_UNARY_VECTOR_TRANSFORMATION(abs, scalar_abs)
 SHARK_UNARY_VECTOR_TRANSFORMATION(log, scalar_log)
 SHARK_UNARY_VECTOR_TRANSFORMATION(exp, scalar_exp)
@@ -198,7 +330,6 @@ name (vector_expression<E> const& e, T scalar){ \
 }
 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator+, scalar_add)
 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator-, scalar_subtract2)
-SHARK_VECTOR_SCALAR_TRANSFORMATION(operator*, scalar_multiply2)
 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator/, scalar_divide)
 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator<, scalar_less_than)
 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator<=, scalar_less_equal_than)
@@ -224,11 +355,130 @@ name (T scalar, vector_expression<E> const& e){ \
 }
 SHARK_VECTOR_SCALAR_TRANSFORMATION_2(operator+, scalar_add)
 SHARK_VECTOR_SCALAR_TRANSFORMATION_2(operator-, scalar_subtract1)
-SHARK_VECTOR_SCALAR_TRANSFORMATION_2(operator*, scalar_multiply1)
 SHARK_VECTOR_SCALAR_TRANSFORMATION_2(min, scalar_min)
 SHARK_VECTOR_SCALAR_TRANSFORMATION_2(max, scalar_max)
 #undef SHARK_VECTOR_SCALAR_TRANSFORMATION_2
 
+template<class E1, class E2>
+class vector_addition: public vector_expression<vector_addition<E1,E2> > {
+private:
+	typedef scalar_binary_plus<
+		typename E1::value_type,
+		typename E2::value_type
+	> functor_type;
+public:
+	typedef std::size_t size_type;
+	typedef std::ptrdiff_t difference_type;
+	typedef typename functor_type::result_type value_type;
+	typedef value_type scalar_type;
+	typedef value_type const_reference;
+	typedef value_type reference;
+	typedef value_type const * const_pointer;
+	typedef value_type const*  pointer;
+
+	typedef typename E1::index_type index_type;
+	typedef typename E1::const_index_pointer const_index_pointer;
+	typedef typename index_pointer<E1>::type index_pointer;
+
+	typedef typename E1::const_closure_type expression_closure1_type;
+	typedef typename E2::const_closure_type expression_closure2_type;
+	
+	typedef vector_addition<E1,E2> const_closure_type;
+	typedef const_closure_type closure_type;
+	typedef unknown_storage_tag storage_category;
+	typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
+
+	// Construction and destruction
+	explicit vector_addition (
+		expression_closure1_type e1, 
+		expression_closure2_type e2
+	):m_expression1(e1),m_expression2(e2){
+		SIZE_CHECK(e1.size() == e2.size());
+	}
+
+	// Accessors
+	size_type size() const {
+		return m_expression1.size();
+	}
+
+	// Expression accessors
+	expression_closure1_type const& expression1() const {
+		return m_expression1;
+	}
+	expression_closure2_type const& expression2() const {
+		return m_expression2;
+	}
+
+	// Element access
+	const_reference operator() (index_type i) const {
+		SIZE_CHECK(i < size());
+		return m_expression1(i) + m_expression2(i);
+	}
+
+	const_reference operator[] (index_type i) const {
+		SIZE_CHECK(i < size());
+		return m_expression1(i) + m_expression2(i);
+	}
+	
+	//computation kernels
+	template<class VecX>
+	void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		assign(x,alpha*m_expression1);
+		plus_assign(x,alpha*m_expression2);
+	}
+	template<class VecX>
+	void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		plus_assign(x,alpha*m_expression1);
+		plus_assign(x,alpha*m_expression2);
+	}
+	
+	template<class VecX>
+	void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
+		minus_assign(x,alpha*m_expression1);
+		minus_assign(x,alpha*m_expression2);
+	}
+
+	// Iterator types
+	typedef binary_transform_iterator<
+		typename E1::const_iterator,
+		typename E2::const_iterator,
+		functor_type
+	> const_iterator;
+	typedef const_iterator iterator;
+
+	const_iterator begin () const {
+		return const_iterator(functor_type(),
+			m_expression1.begin(),m_expression1.end(),
+			m_expression2.begin(),m_expression2.end()
+		);
+	}
+	const_iterator end() const {
+		return const_iterator(functor_type(),
+			m_expression1.end(),m_expression1.end(),
+			m_expression2.end(),m_expression2.end()
+		);
+	}
+
+private:
+	expression_closure1_type m_expression1;
+	expression_closure2_type m_expression2;
+};
+
+template<class E1, class E2>
+vector_addition<E1, E2 > operator+ (
+	vector_expression<E1> const& e1,
+	vector_expression<E2> const& e2
+){
+	return vector_addition<E1, E2>(e1(),e2());
+}
+template<class E1, class E2>
+vector_addition<E1, vector_scalar_multiply<E2> > operator- (
+	vector_expression<E1> const& e1,
+	vector_expression<E2> const& e2
+){
+	return vector_addition<E1, vector_scalar_multiply<E2> >(e1(),-e2());
+}
+
 
 template<class E1, class E2, class F>
 class vector_binary:
@@ -297,7 +547,9 @@ public:
 	// Iterator enhances the iterator of the referenced expressions
 	// with the unary functor.
 	typedef binary_transform_iterator<
-		typename E1::const_iterator,typename E2::const_iterator,functor_type
+		typename E1::const_iterator,
+		typename E2::const_iterator,
+		functor_type
 	> const_iterator;
 	typedef const_iterator iterator;
 
@@ -327,8 +579,6 @@ name(vector_expression<E1> const& e1, vector_expression<E2> const& e2){\
 	typedef F<typename E1::value_type, typename E2::value_type> functor_type;\
 	return vector_binary<E1, E2, functor_type>(e1(),e2(), functor_type());\
 }
-SHARK_BINARY_VECTOR_EXPRESSION(operator+, scalar_binary_plus)
-SHARK_BINARY_VECTOR_EXPRESSION(operator-, scalar_binary_minus)
 SHARK_BINARY_VECTOR_EXPRESSION(operator*, scalar_binary_multiply)
 SHARK_BINARY_VECTOR_EXPRESSION(element_prod, scalar_binary_multiply)
 SHARK_BINARY_VECTOR_EXPRESSION(operator/, scalar_binary_divide)
@@ -417,14 +667,14 @@ soft_max(const vector_expression<E> &e) {
 template<class E>
 typename real_traits<typename E::value_type >::type
 norm_1(const vector_expression<E> &e) {
-	return sum(abs(e));
+	return sum(abs(eval_block(e)));
 }
 
 /// \brief norm_2 v = sum_i |v_i|^2
 template<class E>
 typename real_traits<typename E::value_type >::type
 norm_sqr(const vector_expression<E> &e) {
-	return sum(abs_sqr(e));
+	return sum(abs_sqr(eval_block(e)));
 }
 
 /// \brief norm_2 v = sqrt (sum_i |v_i|^2 )
@@ -439,13 +689,13 @@ norm_2(const vector_expression<E> &e) {
 template<class E>
 typename real_traits<typename E::value_type >::type
 norm_inf(vector_expression<E> const &e){
-	return max(abs(e));
+	return max(abs(eval_block(e)));
 }
 
 /// \brief index_norm_inf v = arg max_i |v_i|
 template<class E>
 std::size_t index_norm_inf(vector_expression<E> const &e){
-	return arg_max(abs(e));
+	return arg_max(abs(eval_block(e)));
 }
 
 // inner_prod (v1, v2) = sum_i v1_i * v2_i
@@ -463,7 +713,7 @@ inner_prod(
 		typename E2::value_type
 	>::promote_type value_type;
 	value_type result = value_type();
-	kernels::dot(e1,e2,result);
+	kernels::dot(eval_block(e1),eval_block(e2),result);
 	return result;
 }
 

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/shark.git

