[shark] 54/79: reorganized asignment operators by centralizing them through freestanding operators. also theoretical upport for blockwie expression is added albeit no expression is declared blockwise right now
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Thu Nov 26 15:41:01 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 2b9bd3ae659d246e6a68360f3c92eee147e182e2
Author: Oswin <oswin.krause at di.ku.dk>
Date: Fri Oct 23 15:51:54 2015 +0200
reorganized asignment operators by centralizing them through freestanding operators. also theoretical upport for blockwie expression is added albeit no expression is declared blockwise right now
---
include/shark/LinAlg/BLAS/assignment.hpp | 564 +++++++++++++++++++++
include/shark/LinAlg/BLAS/detail/traits.hpp | 28 +-
include/shark/LinAlg/BLAS/expression_types.hpp | 94 ----
.../shark/LinAlg/BLAS/kernels/matrix_assign.hpp | 13 +-
.../shark/LinAlg/BLAS/kernels/vector_assign.hpp | 2 +-
include/shark/LinAlg/BLAS/lu.hpp | 4 +-
include/shark/LinAlg/BLAS/matrix.hpp | 124 +----
include/shark/LinAlg/BLAS/matrix_expression.hpp | 10 +-
include/shark/LinAlg/BLAS/matrix_proxy.hpp | 546 +-------------------
include/shark/LinAlg/BLAS/matrix_sparse.hpp | 53 +-
include/shark/LinAlg/BLAS/triangular_matrix.hpp | 109 +---
include/shark/LinAlg/BLAS/vector.hpp | 216 +-------
include/shark/LinAlg/BLAS/vector_expression.hpp | 3 +
include/shark/LinAlg/BLAS/vector_proxy.hpp | 221 +-------
include/shark/LinAlg/BLAS/vector_sparse.hpp | 88 +---
15 files changed, 676 insertions(+), 1399 deletions(-)
diff --git a/include/shark/LinAlg/BLAS/assignment.hpp b/include/shark/LinAlg/BLAS/assignment.hpp
new file mode 100644
index 0000000..d64f5f0
--- /dev/null
+++ b/include/shark/LinAlg/BLAS/assignment.hpp
@@ -0,0 +1,564 @@
+/*!
+ *
+ *
+ * \brief Matrix proxy classes.
+ *
+ *
+ *
+ * \author O. Krause
+ * \date 2013
+ *
+ *
+ * \par Copyright 1995-2015 Shark Development Team
+ *
+ * <BR><HR>
+ * This file is part of Shark.
+ * <http://image.diku.dk/shark/>
+ *
+ * Shark is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Shark is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Shark. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef SHARK_LINALG_BLAS_ASSIGNMENT_HPP
+#define SHARK_LINALG_BLAS_ASSIGNMENT_HPP
+
+#include "kernels/matrix_assign.hpp"
+#include "kernels/vector_assign.hpp"
+#include "detail/traits.hpp"
+
+namespace shark {
+namespace blas {
+
+/////////////////////////////////////////////////////////////////////////////////////
+////// Vector Assign
+////////////////////////////////////////////////////////////////////////////////////
+
+namespace detail{
+ template<class VecX, class VecV>
+ void assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
+ kernels::assign(x,v);
+ }
+ template<class VecX, class VecV>
+ void assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
+ v().assign_to(x);
+ }
+ template<class VecX, class VecV>
+ void plus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
+ kernels::assign<scalar_plus_assign> (x, v);
+ }
+ template<class VecX, class VecV>
+ void plus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
+ x().plus_assign_to(v);
+ }
+ template<class VecX, class VecV>
+ void minus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
+ kernels::assign<scalar_minus_assign> (x, v);
+ }
+ template<class VecX, class VecV>
+ void minus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
+ x().minus_assign_to(v);
+ }
+ template<class VecX, class VecV>
+ void multiply_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
+ kernels::assign<scalar_multiply_assign> (x, v);
+ }
+ template<class VecX, class VecV>
+ void multiply_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
+ typename vector_temporary<VecX>::type temporary(v);
+ kernels::assign<scalar_multiply_assign> (x, temporary);
+ }
+ template<class VecX, class VecV>
+ void divide_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
+ kernels::assign<scalar_divide_assign> (x, v);
+ }
+ template<class VecX, class VecV>
+ void divide_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
+ typename vector_temporary<VecX>::type temporary(v);
+ kernels::assign<scalar_divide_assign> (x, temporary);
+ }
+}
+
+
+/// \brief Dispatches vector assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class VecX, class VecV>
+VecX& assign(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ detail::assign(x,v,typename VecV::evaluation_category());
+ return x();
+}
+
+/// \brief Dispatches vector plus-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class VecX, class VecV>
+VecX& plus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ detail::plus_assign(x,v,typename VecV::evaluation_category());
+ return x();
+}
+
+/// \brief Dispatches vector minus-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class VecX, class VecV>
+VecX& minus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ detail::minus_assign(x,v,typename VecV::evaluation_category());
+ return x();
+}
+
+/// \brief Dispatches vector multiply-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class VecX, class VecV>
+VecX& multiply_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ detail::multiply_assign(x,v,typename VecV::evaluation_category());
+ return x();
+}
+
+/// \brief Dispatches vector multiply-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class VecX, class VecV>
+VecX& divide_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ detail::divide_assign(x,v,typename VecV::evaluation_category());
+ return x();
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+////// Matrix Assign
+////////////////////////////////////////////////////////////////////////////////////
+
+namespace detail{
+ template<class MatA, class MatB>
+ void assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
+ kernels::assign(A,B);
+ }
+ template<class MatA, class MatB>
+ void assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
+ B().assign_to(A);
+ }
+ template<class MatA, class MatB>
+ void plus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
+ kernels::assign<scalar_plus_assign> (A, B);
+ }
+ template<class MatA, class MatB>
+ void plus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
+ B().plus_assign_to(A);
+ }
+ template<class MatA, class MatB>
+ void minus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
+ kernels::assign<scalar_minus_assign> (A, B);
+ }
+ template<class MatA, class MatB>
+ void minus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
+ B().minus_assign_to(A);
+ }
+ template<class MatA, class MatB>
+ void multiply_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
+ kernels::assign<scalar_multiply_assign> (A, B);
+ }
+ template<class MatA, class MatB>
+ void multiply_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
+ typename matrix_temporary<MatA>::type temporary(B);
+ kernels::assign<scalar_multiply_assign> (A, B);
+ }
+ template<class MatA, class MatB>
+ void divide_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
+ kernels::assign<scalar_divide_assign> (A, B);
+ }
+ template<class MatA, class MatB>
+ void divide_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
+ typename matrix_temporary<MatA>::type temporary(B);
+ kernels::assign<scalar_divide_assign> (A, B);
+ }
+}
+
+
+/// \brief Dispatches matrix assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class MatA, class MatB>
+MatA& assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ detail::assign(A,B, typename MatB::evaluation_category());
+ return A();
+}
+
+/// \brief Dispatches matrix plus-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class MatA, class MatB>
+MatA& plus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ detail::plus_assign(A,B, typename MatB::evaluation_category());
+ return A();
+}
+
+/// \brief Dispatches matrix plus-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class MatA, class MatB>
+MatA& minus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ detail::minus_assign(A,B, typename MatB::evaluation_category());
+ return A();
+}
+
+/// \brief Dispatches matrix multiply-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class MatA, class MatB>
+MatA& multiply_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ detail::multiply_assign(A,B, typename MatB::evaluation_category());
+ return A();
+}
+
+/// \brief Dispatches matrix divide-assignment on an expression level
+///
+/// This dispatcher takes care for whether the blockwise evaluation
+/// or the elementwise evaluation is called
+template<class MatA, class MatB>
+MatA& divide_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ detail::divide_assign(A,B, typename MatB::evaluation_category());
+ return A();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+///// Vector Operators
+/////////////////////////////////////////////////////////////////////////////////////
+
+/// \brief Add-Assigns two vector expressions
+///
+/// Performs the operation x_i+=v_i for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(x)+=v to avoid this if A and B do not alias
+template<class VecX, class VecV>
+VecX& operator+=(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ typename vector_temporary<VecX>::type temporary(v);
+ return plus_assign(x,temporary);
+}
+
+/// \brief Subtract-Assigns two vector expressions
+///
+/// Performs the operation x_i-=v_i for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(x)-=v to avoid this if A and B do not alias
+template<class VecX, class VecV>
+VecX& operator-=(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ typename vector_temporary<VecX>::type temporary(v);
+ return minus_assign(x,temporary);
+}
+
+/// \brief Multiply-Assigns two vector expressions
+///
+/// Performs the operation x_i*=v_i for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(x)*=v to avoid this if A and B do not alias
+template<class VecX, class VecV>
+VecX& operator*=(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ typename vector_temporary<VecX>::type temporary(v);
+ return multiply_assign(x,temporary);
+}
+
+/// \brief Divide-Assigns two vector expressions
+///
+/// Performs the operation x_i/=v_i for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(x)/=v to avoid this if A and B do not alias
+template<class VecX, class VecV>
+VecX& operator/=(vector_expression<VecX>& x, vector_expression<VecV> const& v){
+ SIZE_CHECK(x().size() == v().size());
+ typename vector_temporary<VecX>::type temporary(v);
+ return divide_assign(x,temporary);
+}
+
+/// \brief Adds a scalar to all elements of the vector
+///
+/// Performs the operation x_i += t for all elements.
+template<class VecX>
+VecX& operator+=(vector_expression<VecX>& x, typename VecX::scalar_type t){
+ kernels::assign<scalar_plus_assign> (x, t);
+ return x();
+}
+
+/// \brief Subtracts a scalar from all elements of the vector
+///
+/// Performs the operation x_i += t for all elements.
+template<class VecX>
+VecX& operator-=(vector_expression<VecX>& x, typename VecX::scalar_type t){
+ kernels::assign<scalar_minus_assign> (x, t);
+ return x();
+}
+
+/// \brief Multiplies a scalar with all elements of the vector
+///
+/// Performs the operation x_i *= t for all elements.
+template<class VecX>
+VecX& operator*=(vector_expression<VecX>& x, typename VecX::scalar_type t){
+ kernels::assign<scalar_multiply_assign> (x, t);
+ return x();
+}
+
+/// \brief Divides all elements of the vector by a scalar
+///
+/// Performs the operation x_i /= t for all elements.
+template<class VecX>
+VecX& operator/=(vector_expression<VecX>& x, typename VecX::scalar_type t){
+ kernels::assign<scalar_divide_assign> (x, t);
+ return x();
+}
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+///// Matrix Operators
+/////////////////////////////////////////////////////////////////////////////////////
+
+/// \brief Add-Assigns two matrix expressions
+///
+/// Performs the operation A_ij+=B_ij for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(A)+=B to avoid this if A and B do not alias
+template<class MatA, class MatB>
+MatA& operator+=(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ typename matrix_temporary<MatA>::type temporary(B);
+ return plus_assign(A,temporary);
+}
+
+/// \brief Subtract-Assigns two matrix expressions
+///
+/// Performs the operation A_ij-=B_ij for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(A)-=B to avoid this if A and B do not alias
+template<class MatA, class MatB>
+MatA& operator-=(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ typename matrix_temporary<MatA>::type temporary(B);
+ return minus_assign(A,temporary);
+}
+
+/// \brief Multiply-Assigns two matrix expressions
+///
+/// Performs the operation A_ij*=B_ij for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(A)*=B to avoid this if A and B do not alias
+template<class MatA, class MatB>
+MatA& operator*=(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ typename matrix_temporary<MatA>::type temporary(B);
+ return multiply_assign(A,temporary);
+}
+
+/// \brief Divide-Assigns two matrix expressions
+///
+/// Performs the operation A_ij/=B_ij for all elements.
+/// Assumes that the right and left hand side aliases and therefore
+/// performs a copy of the right hand side before assigning
+/// use noalias as in noalias(A)/=B to avoid this if A and B do not alias
+template<class MatA, class MatB>
+MatA& operator/=(matrix_expression<MatA>& A, matrix_expression<MatB> const& B){
+ SIZE_CHECK(A().size1() == B().size1());
+ SIZE_CHECK(A().size2() == B().size2());
+ typename matrix_temporary<MatA>::type temporary(B);
+ return divide_assign(A,temporary);
+}
+
+/// \brief Adds a scalar to all elements of the matrix
+///
+/// Performs the operation A_ij += t for all elements.
+template<class MatA>
+MatA& operator+=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
+ kernels::assign<scalar_plus_assign> (A, t);
+ return A();
+}
+
+/// \brief Subtracts a scalar from all elements of the matrix
+///
+/// Performs the operation A_ij -= t for all elements.
+template<class MatA>
+MatA& operator-=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
+ kernels::assign<scalar_minus_assign> (A, t);
+ return A();
+}
+
+/// \brief Multiplies a scalar to all elements of the matrix
+///
+/// Performs the operation A_ij *= t for all elements.
+template<class MatA>
+MatA& operator*=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
+ kernels::assign<scalar_multiply_assign> (A, t);
+ return A();
+}
+
+/// \brief Divides all elements of the matrix by a scalar
+///
+/// Performs the operation A_ij /= t for all elements.
+template<class MatA>
+MatA& operator /=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
+ kernels::assign<scalar_divide_assign> (A, t);
+ return A();
+}
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+///// Temporary Proxy Operators
+/////////////////////////////////////////////////////////////////////////////////////
+
+template<class T, class U>
+temporary_proxy<T> operator+=(temporary_proxy<T> x, U const& arg){
+ static_cast<T&>(x) += arg;
+ return x;
+}
+template<class T, class U>
+temporary_proxy<T> operator-=(temporary_proxy<T> x, U const& arg){
+ static_cast<T&>(x) -= arg;
+ return x;
+}
+template<class T, class U>
+temporary_proxy<T> operator*=(temporary_proxy<T> x, U const& arg){
+ static_cast<T&>(x) *= arg;
+ return x;
+}
+template<class T, class U>
+temporary_proxy<T> operator/=(temporary_proxy<T> x, U const& arg){
+ static_cast<T&>(x) /= arg;
+ return x;
+}
+
+
+
+
+// Assignment proxy.
+// Provides temporary free assigment when LHS has no alias on RHS
+template<class C>
+class noalias_proxy{
+public:
+ typedef typename C::closure_type closure_type;
+ typedef typename C::scalar_type scalar_type;
+
+ noalias_proxy(C &lval): m_lval(lval) {}
+
+ noalias_proxy(const noalias_proxy &p):m_lval(p.m_lval) {}
+
+ template <class E>
+ closure_type &operator= (const E &e) {
+ return assign(m_lval, e);
+ }
+
+ template <class E>
+ closure_type &operator+= (const E &e) {
+ return plus_assign(m_lval, e);
+ }
+
+ template <class E>
+ closure_type &operator-= (const E &e) {
+ return minus_assign(m_lval, e);
+ }
+
+ template <class E>
+ closure_type &operator*= (const E &e) {
+ return multiply_assign(m_lval, e);
+ }
+
+ template <class E>
+ closure_type &operator/= (const E &e) {
+ return divide_assign(m_lval, e);
+ }
+
+ //this is not needed, but prevents errors when for example doing noalias(x)+=2;
+ closure_type &operator+= (scalar_type t) {
+ return m_lval += t;
+ }
+
+ //this is not needed, but prevents errors when for example doing noalias(x)-=2;
+ closure_type &operator-= (scalar_type t) {
+ return m_lval -= t;
+ }
+
+ //this is not needed, but prevents errors when for example doing noalias(x)*=2;
+ closure_type &operator*= (scalar_type t) {
+ return m_lval *= t;
+ }
+
+ //this is not needed, but prevents errors when for example doing noalias(x)/=2;
+ closure_type &operator/= (scalar_type t) {
+ return m_lval /= t;
+ }
+
+private:
+ closure_type m_lval;
+};
+
+// Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
+// noalias(lhs) = rhs_expression
+template <class C>
+noalias_proxy<C> noalias(matrix_expression<C>& lvalue) {
+ return noalias_proxy<C> (lvalue());
+}
+template <class C>
+noalias_proxy<C> noalias(vector_expression<C>& lvalue) {
+ return noalias_proxy<C> (lvalue());
+}
+
+template <class C>
+noalias_proxy<C> noalias(matrix_set_expression<C>& lvalue) {
+ return noalias_proxy<C> (lvalue());
+}
+template <class C>
+noalias_proxy<C> noalias(vector_set_expression<C>& lvalue) {
+ return noalias_proxy<C> (lvalue());
+}
+template <class C>
+noalias_proxy<C> noalias(temporary_proxy<C> lvalue) {
+ return noalias_proxy<C> (static_cast<C&>(lvalue));
+}
+
+}}
+#endif
\ No newline at end of file
diff --git a/include/shark/LinAlg/BLAS/detail/traits.hpp b/include/shark/LinAlg/BLAS/detail/traits.hpp
index 958a17b..d81af9f 100644
--- a/include/shark/LinAlg/BLAS/detail/traits.hpp
+++ b/include/shark/LinAlg/BLAS/detail/traits.hpp
@@ -298,15 +298,27 @@ struct sparse_tag:public unknown_storage_tag{};
struct dense_tag: public unknown_storage_tag{};
struct packed_tag: public unknown_storage_tag{};
+//evaluation tags
+struct elementwise_tag{};
+struct blockwise_tag{};
+
+namespace detail{
+ template<class S1, class S2>
+ struct evaluation_restrict_traits {
+ typedef S1 type;
+ };
+ template<>
+ struct evaluation_restrict_traits<elementwise_tag, blockwise_tag> {
+ typedef blockwise_tag type;
+ };
+}
+
+template<class E1, class E2>
+struct evaluation_restrict_traits: public detail::evaluation_restrict_traits<
+ typename E1::evaluation_category,
+ typename E2::evaluation_category
+>{};
-template<class S1, class S2>
-struct storage_restrict_traits {
- typedef S1 storage_category;
-};
-template<>
-struct storage_restrict_traits<dense_tag, sparse_tag> {
- typedef sparse_tag storage_category;
-};
template<class E>
struct closure: public boost::mpl::if_<
diff --git a/include/shark/LinAlg/BLAS/expression_types.hpp b/include/shark/LinAlg/BLAS/expression_types.hpp
index 5c6e323..2492f2d 100644
--- a/include/shark/LinAlg/BLAS/expression_types.hpp
+++ b/include/shark/LinAlg/BLAS/expression_types.hpp
@@ -149,100 +149,6 @@ struct temporary_proxy:public P{
}
};
-// Assignment proxy.
-// Provides temporary free assigment when LHS has no alias on RHS
-template<class C>
-class noalias_proxy{
-public:
- typedef typename C::closure_type closure_type;
- typedef typename C::scalar_type scalar_type;
-
- noalias_proxy(C &lval): m_lval(lval) {}
-
- noalias_proxy(const noalias_proxy &p):m_lval(p.m_lval) {}
-
- template <class E>
- closure_type &operator= (const E &e) {
- m_lval.assign(e);
- return m_lval;
- }
-
- template <class E>
- closure_type &operator+= (const E &e) {
- m_lval.plus_assign(e);
- return m_lval;
- }
-
- template <class E>
- closure_type &operator-= (const E &e) {
- m_lval.minus_assign(e);
- return m_lval;
- }
-
- template <class E>
- closure_type &operator*= (const E &e) {
- m_lval.multiply_assign(e);
- return m_lval;
- }
-
- template <class E>
- closure_type &operator/= (const E &e) {
- m_lval.divide_assign(e);
- return m_lval;
- }
-
- //this is not needed, but prevents errors when for example doing noalias(x)+=2;
- closure_type &operator+= (scalar_type t) {
- m_lval += t;
- return m_lval;
- }
-
- //this is not needed, but prevents errors when for example doing noalias(x)-=2;
- closure_type &operator-= (scalar_type t) {
- m_lval -=t;
- return m_lval;
- }
-
- //this is not needed, but prevents errors when for example doing noalias(x)*=2;
- closure_type &operator*= (scalar_type t) {
- m_lval *= t;
- return m_lval;
- }
-
- //this is not needed, but prevents errors when for example doing noalias(x)/=2;
- closure_type &operator/= (scalar_type t) {
- m_lval *=t;
- return m_lval;
- }
-
-private:
- closure_type m_lval;
-};
-
-// Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
-// noalias(lhs) = rhs_expression
-template <class C>
-noalias_proxy<C> noalias(matrix_expression<C>& lvalue) {
- return noalias_proxy<C> (lvalue());
-}
-template <class C>
-noalias_proxy<C> noalias(vector_expression<C>& lvalue) {
- return noalias_proxy<C> (lvalue());
-}
-
-template <class C>
-noalias_proxy<C> noalias(matrix_set_expression<C>& lvalue) {
- return noalias_proxy<C> (lvalue());
-}
-template <class C>
-noalias_proxy<C> noalias(vector_set_expression<C>& lvalue) {
- return noalias_proxy<C> (lvalue());
-}
-template <class C>
-noalias_proxy<C> noalias(temporary_proxy<C> lvalue) {
- return noalias_proxy<C> (static_cast<C&>(lvalue));
-}
-
}
}
diff --git a/include/shark/LinAlg/BLAS/kernels/matrix_assign.hpp b/include/shark/LinAlg/BLAS/kernels/matrix_assign.hpp
index fcf31d7..df6c458 100644
--- a/include/shark/LinAlg/BLAS/kernels/matrix_assign.hpp
+++ b/include/shark/LinAlg/BLAS/kernels/matrix_assign.hpp
@@ -32,6 +32,7 @@ public:
typedef internal_transpose_proxy<M> closure_type;
typedef typename M::orientation::transposed_orientation orientation;
typedef typename M::storage_category storage_category;
+ typedef typename M::evaluation_category evaluation_category;
// Construction and destruction
explicit internal_transpose_proxy(matrix_closure_type m):
@@ -189,7 +190,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size1(); ++i){
matrix_row<M> rowM(m(),i);
- assign(rowM,row(e,i));
+ kernels::assign(rowM,row(e,i));
}
}
@@ -240,7 +241,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size2(); ++i){
matrix_column<M> columnM(m(),i);
- assign(columnM,column(e,i));
+ kernels::assign(columnM,column(e,i));
}
}
@@ -254,7 +255,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size1(); ++i){
matrix_column<M> rowM(m(),i);
- assign(rowM,row(e,i));
+ kernels::assign(rowM,row(e,i));
}
}
@@ -419,7 +420,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size1(); ++i){
matrix_row<M> rowM(m(),i);
- assign<F>(rowM,row(e,i));
+ kernels::assign<F>(rowM,row(e,i));
}
}
@@ -471,7 +472,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size2(); ++i){
matrix_column<M> columnM(m(),i);
- assign<F>(columnM,column(e,i));
+ kernels::assign<F>(columnM,column(e,i));
}
}
@@ -484,7 +485,7 @@ void assign(
) {
for(std::size_t i = 0; i != m().size1(); ++i){
matrix_column<M> rowM(m(),i);
- assign<F>(rowM,row(e,i));
+ kernels::assign<F>(rowM,row(e,i));
}
}
diff --git a/include/shark/LinAlg/BLAS/kernels/vector_assign.hpp b/include/shark/LinAlg/BLAS/kernels/vector_assign.hpp
index e784b73..e9ae3a8 100644
--- a/include/shark/LinAlg/BLAS/kernels/vector_assign.hpp
+++ b/include/shark/LinAlg/BLAS/kernels/vector_assign.hpp
@@ -379,7 +379,7 @@ void assign(
typename vector_temporary<V>::type temporary(v());
assign_sparse(temporary,e, f);
v().clear();
- assign(v, temporary);
+ kernels::assign(v, temporary);
}
// Dispatcher
diff --git a/include/shark/LinAlg/BLAS/lu.hpp b/include/shark/LinAlg/BLAS/lu.hpp
index 70df930..719b7e3 100644
--- a/include/shark/LinAlg/BLAS/lu.hpp
+++ b/include/shark/LinAlg/BLAS/lu.hpp
@@ -27,7 +27,7 @@ typename M::size_type lu_factorize(M &m) {
} else if (singular == 0) {
singular = i + 1;
}
- subrange(m, i + 1, size1, i + 1, size2).minus_assign(
+ noalias(subrange(m, i + 1, size1, i + 1, size2))-=(
outer_prod(subrange(mci, i + 1, size1),
subrange(mri, i + 1, size2)));
}
@@ -61,7 +61,7 @@ typename M::size_type lu_factorize(M &m, PM &pm) {
} else if (singular == 0) {
singular = i + 1;
}
- subrange(m,i + 1, size1, i + 1, size2).minus_assign(
+ noalias(subrange(m,i + 1, size1, i + 1, size2))-=(
outer_prod(subrange(mci, i + 1, size1),
subrange(mri, i + 1, size2)));
}
diff --git a/include/shark/LinAlg/BLAS/matrix.hpp b/include/shark/LinAlg/BLAS/matrix.hpp
index 70f0979..8d1b2ba 100644
--- a/include/shark/LinAlg/BLAS/matrix.hpp
+++ b/include/shark/LinAlg/BLAS/matrix.hpp
@@ -46,6 +46,7 @@ public:
typedef const matrix_reference<const self_type> const_closure_type;
typedef matrix_reference<self_type> closure_type;
typedef dense_tag storage_category;
+ typedef elementwise_tag evaluation_category;
typedef L orientation;
// Construction and destruction
@@ -88,7 +89,7 @@ public:
template<class E>
matrix(matrix_expression<E> const& e):
m_size1(e().size1()), m_size2(e().size2()), m_data(m_size1 * m_size2) {
- assign(e);
+ assign(*this,e);
}
// ---------
@@ -161,51 +162,15 @@ public:
}
// Assignment
-
- template<class E>
- matrix& assign(matrix_expression<E> const& e) {
- kernels::assign(*this,e);
- return *this;
- }
- template<class E>
- matrix& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
-
- template<class E>
- matrix& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
-
/*! @note "pass by value" the key idea to enable move semantics */
matrix& operator = (matrix m) {
swap(m);
return *this;
}
template<class C> // Container assignment without temporary
- matrix& operator = (const matrix_container<C>& m) {
+ matrix& operator = (matrix_container<C> const& m) {
resize(m().size1(), m().size2());
- assign(m);
+ assign(*this, m);
return *this;
}
template<class E>
@@ -214,83 +179,6 @@ public:
swap(temporary);
return *this;
}
- template<class E>
- matrix& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this + e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- matrix& operator += (const matrix_container<C>& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return plus_assign(e);
- }
-
- template<class E>
- matrix& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this - e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- matrix& operator -= (matrix_container<C> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return minus_assign(e);
- }
-
- template<class E>
- matrix& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this * e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- matrix& operator *= (const matrix_container<C>& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return multiply_assign(e);
- }
-
- template<class E>
- matrix& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this / e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- matrix& operator /= (matrix_container<C> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return divide_assign(e);
- }
-
- matrix& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- matrix& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- matrix& operator += (scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- matrix& operator -= (scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
- }
// Swapping
void swap(matrix& m) {
@@ -383,7 +271,6 @@ public:
}
void reserve(size_type non_zeros) {}
-
void reserve_row(std::size_t, std::size_t){}
void reserve_column(std::size_t, std::size_t){}
@@ -442,6 +329,7 @@ public:
typedef const matrix_reference<const self_type> const_closure_type;
typedef matrix_reference<self_type> closure_type;
typedef sparse_tag storage_category;
+ typedef elementwise_tag evaluation_category;
typedef row_major orientation;
// Construction and destruction
@@ -470,7 +358,7 @@ public:
}
// Assignment
- diagonal_matrix& operator = (const diagonal_matrix& m) {
+ diagonal_matrix& operator = (diagonal_matrix const& m) {
m_diagonal = m.m_diagonal;
return *this;
}
diff --git a/include/shark/LinAlg/BLAS/matrix_expression.hpp b/include/shark/LinAlg/BLAS/matrix_expression.hpp
index 91bc864..88b1593 100644
--- a/include/shark/LinAlg/BLAS/matrix_expression.hpp
+++ b/include/shark/LinAlg/BLAS/matrix_expression.hpp
@@ -66,6 +66,7 @@ public:
typedef const_closure_type closure_type;
typedef unknown_orientation orientation;
typedef unknown_storage_tag storage_category;
+ typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
// Construction and destruction
@@ -177,6 +178,7 @@ public:
typedef const_closure_type closure_type;
typedef blas::row_major orientation;
typedef blas::unknown_storage_tag storage_category;
+ typedef typename V::evaluation_category evaluation_category;
// Construction and destruction
explicit vector_repeater (expression_type const& e, std::size_t rows):
@@ -267,6 +269,7 @@ public:
typedef const_closure_type closure_type;
typedef typename E::orientation orientation;
typedef blas::unknown_storage_tag storage_category;
+ typedef typename E::evaluation_category evaluation_category;
// Construction and destruction
matrix_unary(blas::matrix_expression<E> const &e, F const &functor):
@@ -407,6 +410,7 @@ public:
typedef const_closure_type closure_type;
typedef typename E1::orientation orientation;
typedef blas::unknown_storage_tag storage_category;
+ typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
typedef F functor_type;
@@ -542,6 +546,7 @@ public:
typedef const self_type const_closure_type;
typedef const_closure_type closure_type;
typedef unknown_storage_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
matrix_vector_binary1(const expression1_type &e1, const expression2_type &e2):
@@ -880,11 +885,7 @@ public:
/** \brief A matrix with all values of type \c T equal to the same value
*
- * Changing one value has the effect of changing all the values. Assigning it to a normal matrix will copy
- * the same value everywhere in this matrix. All accesses are constant time, due to the trivial value.
- *
* \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
- * \tparam ALLOC an allocator for storing the unique value. By default, a standar allocator is used.
*/
template<class T>
class scalar_matrix:
@@ -909,6 +910,7 @@ public:
typedef matrix_reference<self_type> closure_type;
typedef dense_tag storage_category;
typedef unknown_orientation orientation;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
scalar_matrix():
diff --git a/include/shark/LinAlg/BLAS/matrix_proxy.hpp b/include/shark/LinAlg/BLAS/matrix_proxy.hpp
index 086c60a..ad08b93 100644
--- a/include/shark/LinAlg/BLAS/matrix_proxy.hpp
+++ b/include/shark/LinAlg/BLAS/matrix_proxy.hpp
@@ -33,7 +33,7 @@
#ifndef SHARK_LINALG_BLAS_MATRIX_PROXY_HPP
#define SHARK_LINALG_BLAS_MATRIX_PROXY_HPP
-#include "kernels/matrix_assign.hpp"
+#include "assignment.hpp"
namespace shark {
namespace blas {
@@ -60,6 +60,8 @@ public:
typedef self_type closure_type;
typedef typename M::orientation orientation;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
+
// Construction and destruction
matrix_reference(M& m):
@@ -165,95 +167,11 @@ public:
// Assignment
-
- template<class E>
- matrix_reference& assign(matrix_expression<E> const& e) {
- expression().assign(e);
- return *this;
- }
- template<class E>
- matrix_reference& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().plus_assign(e);
- return *this;
- }
- template<class E>
- matrix_reference& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().minus_assign(e);
- return *this;
- }
- template<class E>
- matrix_reference& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().multiply_assign(e);
- return *this;
- }
- template<class E>
- matrix_reference& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().divide_assign(e);
- return *this;
- }
-
- matrix_reference& operator = (const matrix_reference& m) {
- expression() = m();
- return *this;
- }
template<class E>
matrix_reference& operator = (matrix_expression<E> const& e) {
expression() = e();
return *this;
}
- template<class E>
- matrix_reference& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() +=e();
- return *this;
- }
- template<class E>
- matrix_reference& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() -= e();
- return *this;
- }
- template<class E>
- matrix_reference& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() *= e();
- return *this;
- }
- template<class E>
- matrix_reference& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() /= e();
- return *this;
- }
- matrix_reference& operator *= (scalar_type t) {
- expression() *= t;
- return *this;
- }
- matrix_reference& operator /= (scalar_type t) {
- expression() /= t;
- return *this;
- }
-
- matrix_reference& operator += (scalar_type t) {
- expression() += t;
- return *this;
- }
- matrix_reference& operator -= (scalar_type t) {
- expression() -= t;
- return *this;
- }
// Iterator types
typedef typename row_iterator<M>::type row_iterator;
@@ -358,6 +276,7 @@ public:
typedef matrix_transpose<M> closure_type;
typedef typename M::orientation::transposed_orientation orientation;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
explicit matrix_transpose(matrix_closure_type m):
@@ -521,41 +440,6 @@ public:
// Assignment
//we implement it by using the identity A^T op= B <=> A op= B^T where op= is one of =,-=,+=
-
- template<class E>
- matrix_transpose& assign(matrix_expression<E> const& e) {
- expression().assign(matrix_transpose<E const>(e()));
- return *this;
- }
- template<class E>
- matrix_transpose& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().plus_assign(matrix_transpose<E const>(e()));
- return *this;
- }
- template<class E>
- matrix_transpose& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().minus_assign(matrix_transpose<E const>(e()));
- return *this;
- }
- template<class E>
- matrix_transpose& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().multiply_assign(matrix_transpose<E const>(e()));
- return *this;
- }
- template<class E>
- matrix_transpose& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression().divide_assign(matrix_transpose<E const>(e()));
- return *this;
- }
-
matrix_transpose& operator = (matrix_transpose const& m) {
expression() = m.expression();
return *this;
@@ -565,52 +449,6 @@ public:
expression() = matrix_transpose<E const>(e());
return *this;
}
- template<class E>
- matrix_transpose& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() += matrix_transpose<E const>(e());
- return *this;
- }
- template<class E>
- matrix_transpose& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() -= matrix_transpose<E const>(e());
- return *this;
- }
- template<class E>
- matrix_transpose& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() *= matrix_transpose<E const>(e());
- return *this;
- }
- template<class E>
- matrix_transpose& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- expression() /= matrix_transpose<E const>(e());
- return *this;
- }
-
- matrix_transpose& operator *= (scalar_type t) {
- expression() *= t;
- return *this;
- }
- matrix_transpose& operator /= (scalar_type t) {
- expression() /= t;
- return *this;
- }
-
- matrix_transpose& operator += (scalar_type t) {
- expression() += t;
- return *this;
- }
- matrix_transpose& operator -= (scalar_type t) {
- expression() -= t;
- return *this;
- }
private:
matrix_closure_type m_expression;
};
@@ -653,6 +491,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
matrix_row(matrix_closure_type const& expression, index_type i):m_expression(expression), m_i(i) {
@@ -727,82 +566,15 @@ public:
}
// Assignment
- template<class E>
- matrix_row& assign(vector_expression<E> const& e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- matrix_row& plus_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_row& minus_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_row& multiply_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_row& divide_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
template<class E>
matrix_row& operator = (vector_expression<E> const& e) {
- return assign(typename vector_temporary<M>::type(e));
+ return assign(*this, typename vector_temporary<M>::type(e));
}
matrix_row& operator = (matrix_row const& e) {
- return assign(typename vector_temporary<M>::type(e));
- }
- template<class E>
- matrix_row& operator += (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- return plus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_row& operator -= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_row& operator *= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_row& operator /= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- return divide_assign(typename vector_temporary<E>::type(e));
+ return assign(*this, typename vector_temporary<M>::type(e));
}
- matrix_row& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- matrix_row& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- matrix_row& operator += (scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- matrix_row& operator -= (scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
- }
-
// Iterator types
typedef typename M::const_row_iterator const_iterator;
typedef typename row_iterator<M>::type iterator;
@@ -941,6 +713,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
matrix_column(matrix_closure_type expression, index_type j)
@@ -1017,36 +790,6 @@ public:
// Assignment
template<class E>
- matrix_column& assign(vector_expression<E> const& e) {
- m_wrapper.assign(e());
- return *this;
- }
- template<class E>
- matrix_column& plus_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper.plus_assign(e());
- return *this;
- }
- template<class E>
- matrix_column& minus_assign( vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper.minus_assign(e());
- return *this;
- }
- template<class E>
- matrix_column& multiply_assign( vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper.multiply_assign(e());
- return *this;
- }
- template<class E>
- matrix_column& divide_assign( vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper.divide_assign(e());
- return *this;
- }
-
- template<class E>
matrix_column& operator = (vector_expression<E> const& e) {
m_wrapper = e();
return *this;
@@ -1055,48 +798,6 @@ public:
m_wrapper = e;
return *this;
}
- template<class E>
- matrix_column& operator += (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper += e();
- return *this;
- }
- template<class E>
- matrix_column& operator -= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper -= e();
- return *this;
- }
- template<class E>
- matrix_column& operator *= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper *= e();
- return *this;
- }
- template<class E>
- matrix_column& operator /= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- m_wrapper /= e();
- return *this;
- }
-
- matrix_column& operator *= (scalar_type at) {
- m_wrapper *= at;
- return *this;
- }
- matrix_column& operator /= (scalar_type at) {
- m_wrapper /= at;
- return *this;
- }
-
- matrix_column& operator += (scalar_type at) {
- m_wrapper += at;
- return *this;
- }
- matrix_column& operator -= (scalar_type at) {
- m_wrapper -= at;
- return *this;
- }
// Iterator types
typedef typename wrapper_type::const_iterator const_iterator;
@@ -1178,6 +879,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
matrix_vector_range(matrix_type& expression, range const&r1, range const&r2):
@@ -1244,77 +946,10 @@ public:
}
// Assignment
- template<class E>
- matrix_vector_range& assign(vector_expression<E> const& e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- matrix_vector_range& plus_assign(vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_vector_range& minus_assign(vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_vector_range& multiply_assign(vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- matrix_vector_range& divide_assign(vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
template<class E>
matrix_vector_range& operator = (vector_expression<E> const& e) {
- return assign(typename vector_temporary<M>::type(e));
- }
- template<class E>
- matrix_vector_range& operator += (vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- return plus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_vector_range& operator -= (vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_vector_range& operator *= (vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- matrix_vector_range& operator /= (vector_expression<E>const& e) {
- SIZE_CHECK(size() == e().size());
- return divide_assign(typename vector_temporary<E>::type(e));
- }
-
- matrix_vector_range& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- matrix_vector_range& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- matrix_vector_range& operator += (scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- matrix_vector_range& operator -= (scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
+ return assign(*this, typename vector_temporary<M>::type(e));
}
typedef indexed_iterator<closure_type> iterator;
@@ -1405,6 +1040,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename M::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
typedef typename M::orientation orientation;
// Construction and destruction
@@ -1473,88 +1109,13 @@ public:
}
// Assignment
- template<class E>
- self_type& assign(matrix_expression<E> const& e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- self_type& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
self_type& operator = (self_type const& e) {
- return assign(typename matrix_temporary<self_type>::type(e));
+ return assign(*this, typename matrix_temporary<self_type>::type(e));
}
template<class E>
self_type& operator = (matrix_expression<E> const& e) {
- return assign(typename matrix_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return plus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return divide_assign(typename vector_temporary<E>::type(e));
- }
-
- self_type& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- self_type& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- self_type& operator += (scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- self_type& operator -= (scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
+ return assign(*this, typename matrix_temporary<E>::type(e));
}
// Iterator types
@@ -1755,6 +1316,7 @@ public:
typedef matrix_reference<self_type const> const const_closure_type;
typedef matrix_reference<self_type> closure_type;
typedef dense_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
@@ -1856,93 +1418,17 @@ public:
// -------
// ASSIGNING
// -------
-
- template<class E>
- self_type& assign(matrix_expression<E> const& e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- self_type& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- self_type& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
self_type& operator = (self_type const& e) {
SIZE_CHECK(size1() == e().size1());
SIZE_CHECK(size2() == e().size2());
- return assign(typename matrix_temporary<self_type>::type(e));
+ return assign(*this, typename matrix_temporary<self_type>::type(e));
}
template<class E>
self_type& operator = (matrix_expression<E> const& e) {
SIZE_CHECK(size1() == e().size1());
SIZE_CHECK(size2() == e().size2());
- return assign(typename matrix_temporary<self_type>::type(e));
- }
- template<class E>
- self_type& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return plus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- self_type& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return divide_assign(typename vector_temporary<E>::type(e));
- }
-
- self_type& operator*=(scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- self_type& operator/=(scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- self_type& operator+=(scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- self_type& operator-=(scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
+ return assign(*this, typename matrix_temporary<self_type>::type(e));
}
// --------------
diff --git a/include/shark/LinAlg/BLAS/matrix_sparse.hpp b/include/shark/LinAlg/BLAS/matrix_sparse.hpp
index bc74b14..d29d7f1 100644
--- a/include/shark/LinAlg/BLAS/matrix_sparse.hpp
+++ b/include/shark/LinAlg/BLAS/matrix_sparse.hpp
@@ -86,6 +86,7 @@ public:
typedef matrix_reference<self_type> closure_type;
typedef sparse_tag storage_category;
+ typedef elementwise_tag evaluation_category;
typedef row_major orientation;
// Construction and destruction
@@ -105,7 +106,7 @@ public:
, m_rowStart(e().size1() + 1, 0)
, m_rowEnd(e().size1(), 0)
, m_indices(non_zeros), m_values(non_zeros), m_zero(0) {
- kernels::assign(*this, e);
+ assign(*this, e);
}
// Accessors
@@ -259,59 +260,13 @@ public:
template<class C> // Container assignment without temporary
compressed_matrix &operator = (const matrix_container<C> &m) {
resize(m().size1(), m().size2());
- assign(m);
- return *this;
- }
- compressed_matrix &assign_temporary(compressed_matrix &m) {
- swap(m);
+ assign(*this, m);
return *this;
}
template<class E>
compressed_matrix &operator = (const matrix_expression<E> &e) {
self_type temporary(e, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class E>
- compressed_matrix &assign(const matrix_expression<E> &e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- compressed_matrix &operator += (const matrix_expression<E> &e) {
- self_type temporary(*this + e, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class C> // Container assignment without temporary
- compressed_matrix &operator += (const matrix_container<C> &m) {
- plus_assign(m);
- return *this;
- }
- template<class E>
- compressed_matrix &plus_assign(const matrix_expression<E> &e) {
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- compressed_matrix &operator -= (const matrix_expression<E> &e) {
- self_type temporary(*this - e, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class C> // Container assignment without temporary
- compressed_matrix &operator -= (const matrix_container<C> &m) {
- minus_assign(m);
- return *this;
- }
- template<class E>
- compressed_matrix &minus_assign(const matrix_expression<E> &e) {
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- compressed_matrix &operator *= (value_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- compressed_matrix &operator /= (value_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
+ swap(temporary);
return *this;
}
diff --git a/include/shark/LinAlg/BLAS/triangular_matrix.hpp b/include/shark/LinAlg/BLAS/triangular_matrix.hpp
index 1838405..db5ba1e 100644
--- a/include/shark/LinAlg/BLAS/triangular_matrix.hpp
+++ b/include/shark/LinAlg/BLAS/triangular_matrix.hpp
@@ -34,6 +34,7 @@ public:
typedef const matrix_reference<const self_type> const_closure_type;
typedef matrix_reference<self_type> closure_type;
typedef packed_tag storage_category;
+ typedef elementwise_tag evaluation_category;
typedef packed<Orientation,TriangularType> orientation;
// Construction and destruction
@@ -64,7 +65,7 @@ public:
triangular_matrix(matrix_expression<E> const& e)
:m_size(e().size1()), m_data(m_size * (m_size+1)/2)
{
- assign(e);
+ assign(*this, e);
}
// ---------
@@ -145,42 +146,6 @@ public:
return orientation::non_zero(i,j);
}
- // Assignment
- template<class E>
- triangular_matrix& assign(matrix_expression<E> const& e) {
- kernels::assign(*this,e);
- return *this;
- }
- template<class E>
- triangular_matrix& plus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- triangular_matrix& minus_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
-
- template<class E>
- triangular_matrix& multiply_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- triangular_matrix& divide_assign(matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
-
/*! @note "pass by value" the key idea to enable move semantics */
triangular_matrix& operator = (triangular_matrix m) {
swap(m);
@@ -190,7 +155,7 @@ public:
triangular_matrix& operator = (const matrix_container<C>& m) {
SIZE_CHECK(m().size1()==m().size2());
resize(m().size1());
- assign(m);
+ assign(*this, m);
return *this;
}
template<class E>
@@ -199,74 +164,6 @@ public:
swap(temporary);
return *this;
}
- template<class E>
- triangular_matrix& operator += (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this + e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- triangular_matrix& operator += (const matrix_container<C>& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return plus_assign(e);
- }
-
- template<class E>
- triangular_matrix& operator -= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this - e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- triangular_matrix& operator -= (matrix_container<C> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return minus_assign(e);
- }
-
- template<class E>
- triangular_matrix& operator *= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this * e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- triangular_matrix& operator *= (const matrix_container<C>& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return multiply_assign(e);
- }
-
- template<class E>
- triangular_matrix& operator /= (matrix_expression<E> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- self_type temporary(*this / e);
- swap(temporary);
- return *this;
- }
- template<class C> // Container assignment without temporary
- triangular_matrix& operator /= (matrix_container<C> const& e) {
- SIZE_CHECK(size1() == e().size1());
- SIZE_CHECK(size2() == e().size2());
- return divide_assign(e);
- }
-
- triangular_matrix& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- triangular_matrix& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
// Swapping
void swap(triangular_matrix& m) {
diff --git a/include/shark/LinAlg/BLAS/vector.hpp b/include/shark/LinAlg/BLAS/vector.hpp
index 16f082d..fdd3e01 100644
--- a/include/shark/LinAlg/BLAS/vector.hpp
+++ b/include/shark/LinAlg/BLAS/vector.hpp
@@ -44,6 +44,7 @@ public:
typedef vector_reference<self_type> closure_type;
typedef self_type vector_temporary_type;
typedef dense_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
@@ -72,17 +73,16 @@ public:
/// \brief Copy-constructor of a vector
/// \param v is the vector to be duplicated
- vector(const vector& v):m_storage(v.m_storage) {}
+ vector(vector const& v):m_storage(v.m_storage) {}
/// \brief Copy-constructor of a vector from a vector_expression
/// Depending on the vector_expression, this constructor can have the cost of the computations
/// of the expression (trivial to say it, but it is to take into account in your complexity calculations).
/// \param e the vector_expression which values will be duplicated into the vector
template<class E>
- vector(vector_expression<E> const& e):
- vector_container<self_type> (),
- m_storage(e().size()) {
- kernels::assign (*this, e);
+ vector(vector_expression<E> const& e)
+ :m_storage(e().size()) {
+ assign(*this, e);
}
// ---------
@@ -195,73 +195,20 @@ public:
void clear() {
std::fill(m_storage.begin(), m_storage.end(), value_type/*zero*/());
}
-
+
// -------------------
- // Assignment Functions
+ // Assignment operators
// -------------------
-
- /// \brief Assign the result of a vector_expression to the vector
- /// Assign the result of a vector_expression to the vector.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& assign(vector_expression<E> const& e) {
- kernels::assign (*this, e);
- return *this;
- }
-
- /// \brief Assign the sum of the vector and a vector_expression to the vector
- /// Assign the sum of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& plus_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- /// \brief Assign the difference of the vector and a vector_expression to the vector
- /// Assign the difference of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& minus_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
-
- /// \brief Assign the elementwise product of the vector and a vector_expression to the vector
- /// Assign the difference of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& multiply_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_multiply_assign> (*this, e);
+ /// \brief Assign a full vector (\e RHS-vector) to the current vector (\e LHS-vector)
+ /// Assign a full vector (\e RHS-vector) to the current vector (\e LHS-vector). This method does not create any temporary.
+ /// \param v is the source vector container
+ /// \return a reference to a vector (i.e. the destination vector)
+ vector& operator = (vector const& v) {
+ m_storage = v.m_storage;
return *this;
}
- /// \brief Assign the elementwise division of the vector and a vector_expression to the vector
- /// Assign the difference of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& divide_assign(vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
-
- // -------------------
- // Assignment operators
- // -------------------
-
/// \brief Assign a full vector (\e RHS-vector) to the current vector (\e LHS-vector)
/// Assign a full vector (\e RHS-vector) to the current vector (\e LHS-vector). This method does not create any temporary.
/// \param v is the source vector container
@@ -269,7 +216,7 @@ public:
template<class C> // Container assignment without temporary
vector& operator = (vector_container<C> const& v) {
resize(v().size());
- return assign(v);
+ return assign(*this, v);
}
/// \brief Assign the result of a vector_expression to the vector
@@ -282,141 +229,6 @@ public:
swap(*this,temporary);
return *this;
}
-
- // Assignment
- /// \brief Assign a full vector (\e RHS-vector) to the current vector (\e LHS-vector)
- /// \param v is the source vector
- /// \return a reference to a vector (i.e. the destination vector)
- vector& operator = (vector v) {
- swap(*this,v);
- return *this;
- }
-
-
- /// \brief Assign the sum of the vector and a vector_expression to the vector
- /// Assign the sum of the vector and a vector_expression to the vector.
- /// A temporary is created for the computations.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& operator += (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- self_type temporary(e);
- return plus_assign(temporary);
- }
-
- /// \brief Assign the sum of the vector and a vector_expression to the vector
- /// Assign the sum of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param v is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class C> // Container assignment without temporary
- vector& operator += (vector_container<C> const& v) {
- SIZE_CHECK(size() == v().size());
- return plus_assign(v);
- }
-
- /// \brief Assign the difference of the vector and a vector_expression to the vector
- /// Assign the difference of the vector and a vector_expression to the vector.
- /// A temporary is created for the computations.
- /// \param e is a const reference to the vector_expression
- template<class E>
- vector& operator -= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- self_type temporary(e);
- return minus_assign(temporary);
- }
-
- /// \brief Assign the difference of the vector and a vector_expression to the vector
- /// Assign the difference of the vector and a vector_expression to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param v is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class C> // Container assignment without temporary
- vector& operator -= (vector_container<C> const& v) {
- SIZE_CHECK(size() == v().size());
- return minus_assign(v);
- }
-
- /// \brief Assign the elementwise product of the vector and a vector_expression to the vector
- /// A temporary is created for the computations.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& operator *= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- self_type temporary(e);
- return multiply_assign(temporary);
- }
-
- /// \brief Assign the elementwise product of the vector and a vector_expression to the vector
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param v is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class C> // Container assignment without temporary
- vector& operator *= (vector_container<C> const& v) {
- SIZE_CHECK(size() == v().size());
- return multiply_assign(v);
- }
-
- /// \brief Assign the elementwise division of the vector and a vector_expression to the vector
- /// A temporary is created for the computations.
- /// \param e is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class E>
- vector& operator /= (vector_expression<E> const& e) {
- SIZE_CHECK(size() == e().size());
- self_type temporary(e);
- return divide_assign(temporary);
- }
-
- /// \brief Assign the elementwise product of the vector and a vector_expression to the vector
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param v is a const reference to the vector_expression
- /// \return a reference to the resulting vector
- template<class C> // Container assignment without temporary
- vector& operator /= (vector_container<C> const& v) {
- SIZE_CHECK(size() == v().size());
- return divide_assign(v);
- }
-
- /// \brief Assign the product of the vector and a scalar to the vector
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param t is a const reference to the scalar
- /// \return a reference to the resulting vector
- vector& operator *= (scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
-
- /// \brief Assign the division of the vector by a scalar to the vector
- /// Assign the division of the vector by a scalar to the vector.
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \tparam E is the type of the vector_expression
- /// \param t is a const reference to the scalar
- /// \return a reference to the resulting vector
- vector& operator /= (scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- /// \brief Add the scalar value to every element of the vector
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param t the scalar to add
- /// \return a reference to the resulting vector
- vector& operator += (scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
-
- /// \brief Add the scalar value to every element of the vector
- /// No temporary is created. Computations are done and stored directly into the resulting vector.
- /// \param t the scalar to add
- /// \return a reference to the resulting vector
- vector& operator -= (scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
- }
// Iterator types
typedef dense_storage_iterator<value_type> iterator;
diff --git a/include/shark/LinAlg/BLAS/vector_expression.hpp b/include/shark/LinAlg/BLAS/vector_expression.hpp
index 2ce6a56..ad6f39a 100644
--- a/include/shark/LinAlg/BLAS/vector_expression.hpp
+++ b/include/shark/LinAlg/BLAS/vector_expression.hpp
@@ -31,6 +31,7 @@ public:
typedef const vector_reference<const self_type> const_closure_type;
typedef vector_reference<self_type> closure_type;
typedef unknown_storage_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
scalar_vector()
@@ -116,6 +117,7 @@ public:
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.
@@ -255,6 +257,7 @@ public:
typedef self_type const const_closure_type;
typedef self_type closure_type;
typedef unknown_storage_tag storage_category;
+ typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
// Construction and destruction
explicit vector_binary (
diff --git a/include/shark/LinAlg/BLAS/vector_proxy.hpp b/include/shark/LinAlg/BLAS/vector_proxy.hpp
index e85a2a6..6d58283 100644
--- a/include/shark/LinAlg/BLAS/vector_proxy.hpp
+++ b/include/shark/LinAlg/BLAS/vector_proxy.hpp
@@ -34,7 +34,7 @@
#ifndef SHARK_LINALG_BLAS_VECTOR_PROXY_HPP
#define SHARK_LINALG_BLAS_VECTOR_PROXY_HPP
-#include "kernels/vector_assign.hpp"
+#include "assignment.hpp"
#include "detail/iterator.hpp"
namespace shark{
@@ -62,6 +62,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename V::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
vector_reference(referred_type& v):m_expression(&v){}
@@ -122,86 +123,14 @@ public:
reference operator [](index_type i) const{
return expression() [i];
}
-
- // Assignment
- template<class E>
- vector_reference& assign(vector_expression<E> const& e){
- expression().assign(e);
- return *this;
- }
- template<class E>
- vector_reference& minus_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression().minus_assign(e);
- return *this;
- }
- template<class E>
- vector_reference& plus_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression().plus_assign(e);
- return *this;
- }
- template<class E>
- vector_reference& multiply_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression().multiply_assign(e);
- return *this;
- }
- template<class E>
- vector_reference& divide_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression().divide_assign(e);
- return *this;
- }
vector_reference& operator = (vector_reference const& v){
- expression() = v;
+ expression() = v.expression();
return *this;
}
template<class E>
vector_reference& operator = (vector_expression<E> const& e){
- expression() = e;
- return *this;
- }
- template<class E>
- vector_reference& operator -= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression() -= e();//op() allows for optimization when e is vector_container
- return *this;
- }
- template<class E>
- vector_reference& operator += (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression() += e();
- return *this;
- }
- template<class E>
- vector_reference& operator *= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression() *= e();
- return *this;
- }
- template<class E>
- vector_reference& operator /= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- expression() /= e();
- return *this;
- }
- vector_reference& operator *= (scalar_type t){
- expression() *= t;
- return *this;
- }
- vector_reference& operator /= (scalar_type t){
- expression() /= t;
- return *this;
- }
-
- vector_reference& operator += (scalar_type t){
- expression() += t;
- return *this;
- }
- vector_reference& operator -= (scalar_type t){
- expression() -= t;
+ expression() = e();
return *this;
}
@@ -276,6 +205,7 @@ public:
typedef const self_type const_closure_type;
typedef self_type closure_type;
typedef typename V::storage_category storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
vector_range(vector_closure_type const& data, range const& r):
@@ -331,81 +261,14 @@ public:
return m_expression(m_range(i));
}
- // Assignment
- template<class E>
- vector_range& assign(vector_expression<E> const& e){
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- vector_range& plus_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- vector_range& minus_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- vector_range& multiply_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- vector_range& divide_assign(vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
-
+ // Assignment operators
vector_range& operator = (vector_range const& vr){
- return assign (typename vector_temporary<V>::type(vr));
+ return assign(*this, typename vector_temporary<V>::type(vr));
}
template<class E>
vector_range& operator = (vector_expression<E> const& e){
- return assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- vector_range& operator += (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- return plus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- vector_range& operator -= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- vector_range& operator *= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- vector_range& operator /= (vector_expression<E> const& e){
- SIZE_CHECK(size() == e().size());
- return divide_assign(typename vector_temporary<E>::type(e));
- }
-
- vector_range& operator += (scalar_type t){
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- vector_range& operator -= (scalar_type t){
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
- }
- vector_range& operator *= ( scalar_type t){
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- vector_range& operator /= ( scalar_type t){
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
+ return assign(*this, typename vector_temporary<E>::type(e));
}
typedef subrange_iterator< typename vector_closure_type::iterator> iterator;
@@ -511,10 +374,10 @@ public:
typedef index_type const* const_index_pointer;
typedef index_type* index_pointer;
- //ublas types
typedef vector_reference<self_type const> const const_closure_type;
typedef vector_reference<self_type> closure_type;
typedef dense_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
@@ -624,35 +487,6 @@ public:
(*this)[i] = value_type/*zero*/();
}
- // -------
- // ASSIGNING
- // -------
-
- template<class E>
- dense_vector_adaptor& assign(const vector_expression<E>& e) {
- kernels::assign(*this, e);
- return *this;
- }
- template<class E>
- dense_vector_adaptor& plus_assign(const vector_expression<E>& e) {
- kernels::assign<scalar_plus_assign> (*this, e);
- return *this;
- }
- template<class E>
- dense_vector_adaptor& minus_assign(const vector_expression<E>& e) {
- kernels::assign<scalar_minus_assign> (*this, e);
- return *this;
- }
- template<class E>
- dense_vector_adaptor& multiply_assign(const vector_expression<E>& e) {
- kernels::assign<scalar_multiply_assign> (*this, e);
- return *this;
- }
- template<class E>
- dense_vector_adaptor& divide_assign(const vector_expression<E>& e) {
- kernels::assign<scalar_divide_assign> (*this, e);
- return *this;
- }
dense_vector_adaptor& operator = (dense_vector_adaptor const& e) {
return assign(typename vector_temporary<self_type>::type(e));
@@ -662,41 +496,6 @@ public:
return assign(typename vector_temporary<E>::type(e));
}
- template<class E>
- dense_vector_adaptor& operator += (const vector_expression<E>& e) {
- return plus_assign(typename vector_temporary<self_type>::type(e));
- }
- template<class E>
- dense_vector_adaptor& operator -= (const vector_expression<E>& e) {
- return minus_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- dense_vector_adaptor& operator *= (const vector_expression<E>& e) {
- return multiply_assign(typename vector_temporary<E>::type(e));
- }
- template<class E>
- dense_vector_adaptor& operator /= (const vector_expression<E>& e) {
- return divide_assign(typename vector_temporary<E>::type(e));
- }
-
- dense_vector_adaptor& operator *= ( scalar_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- dense_vector_adaptor& operator /= ( scalar_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
- return *this;
- }
-
- dense_vector_adaptor& operator += ( scalar_type t) {
- kernels::assign<scalar_plus_assign> (*this, t);
- return *this;
- }
- dense_vector_adaptor& operator -= ( scalar_type t) {
- kernels::assign<scalar_minus_assign> (*this, t);
- return *this;
- }
-
// --------------
// ITERATORS
// --------------
@@ -784,8 +583,8 @@ public:
typedef index_type const* const_index_pointer;
typedef const_index_pointer index_pointer;
- //ublas types
typedef sparse_tag storage_category;
+ typedef elementwise_tag evaluation_category;
typedef vector_reference<self_type const> const const_closure_type;
typedef vector_reference<self_type> closure_type;
diff --git a/include/shark/LinAlg/BLAS/vector_sparse.hpp b/include/shark/LinAlg/BLAS/vector_sparse.hpp
index fb1838f..3f9cf21 100644
--- a/include/shark/LinAlg/BLAS/vector_sparse.hpp
+++ b/include/shark/LinAlg/BLAS/vector_sparse.hpp
@@ -29,7 +29,7 @@ namespace blas {
template<class T, class I>
class compressed_vector:public vector_container<compressed_vector<T, I> > {
- typedef T &true_reference;
+ typedef T& true_reference;
typedef compressed_vector<T, I> self_type;
public:
@@ -68,7 +68,7 @@ public:
public:
// Construction and destruction
- reference(self_type &m, size_type i):
+ reference(self_type& m, size_type i):
m_vector(m), m_i(i) {}
// Assignment
@@ -76,7 +76,7 @@ public:
return ref()=d;
}
- value_type& operator=(reference const &v ){
+ value_type& operator=(reference const& v ){
return ref() = v.value();
}
@@ -112,6 +112,7 @@ public:
typedef const vector_reference<const self_type> const_closure_type;
typedef vector_reference<self_type> closure_type;
typedef sparse_tag storage_category;
+ typedef elementwise_tag evaluation_category;
// Construction and destruction
compressed_vector():m_size(0), m_nnz(0),m_indices(1,0),m_zero(0){}
@@ -121,7 +122,7 @@ public:
compressed_vector(vector_expression<AE> const& ae, size_type non_zeros = 0)
:m_size(ae().size()), m_nnz(0), m_indices(non_zeros,0), m_values(non_zeros),m_zero(0)
{
- kernels::assign(*this, ae);
+ assign(*this, ae);
}
// Accessors
@@ -144,22 +145,22 @@ public:
index_type const* indices() const{
if(nnz_capacity() == 0)
return 0;
- return &m_indices[0];
+ return& m_indices[0];
}
index_type* indices(){
if(nnz_capacity() == 0)
return 0;
- return &m_indices[0];
+ return& m_indices[0];
}
value_type const* values() const {
if(nnz_capacity() == 0)
return 0;
- return &m_values[0];
+ return& m_values[0];
}
value_type* values(){
if(nnz_capacity() == 0)
return 0;
- return &m_values[0];
+ return& m_values[0];
}
void resize(size_type size) {
@@ -199,7 +200,7 @@ public:
}
// Assignment
- compressed_vector &operator = (compressed_vector const& v) {
+ compressed_vector& operator = (compressed_vector const& v) {
m_size = v.m_size;
m_nnz = v.m_nnz;
m_indices = v.m_indices;
@@ -207,69 +208,20 @@ public:
return *this;
}
template<class C> // Container assignment without temporary
- compressed_vector &operator = (vector_container<C> const& v) {
+ compressed_vector& operator = (vector_container<C> const& v) {
resize(v().size(), false);
- assign(v);
- return *this;
- }
-
- compressed_vector &assign_temporary(compressed_vector& v) {
- swap(v);
+ assign(*this, v);
return *this;
}
template<class AE>
- compressed_vector &operator = (vector_expression<AE> const& ae) {
+ compressed_vector& operator = (vector_expression<AE> const& ae) {
self_type temporary(ae, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class AE>
- compressed_vector &assign(vector_expression<AE> const& ae) {
- kernels::assign(*this, ae);
- return *this;
- }
-
- // Computed assignment
- template<class AE>
- compressed_vector &operator += (vector_expression<AE> const& ae) {
- self_type temporary(*this + ae, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class C> // Container assignment without temporary
- compressed_vector &operator += (vector_container<C> const& v) {
- plus_assign(v);
- return *this;
- }
- template<class AE>
- compressed_vector &plus_assign(vector_expression<AE> const& ae) {
- kernels::assign<scalar_plus_assign> (*this, ae);
- return *this;
- }
- template<class AE>
- compressed_vector &operator -= (vector_expression<AE> const& ae) {
- self_type temporary(*this - ae, nnz_capacity());
- return assign_temporary(temporary);
- }
- template<class C> // Container assignment without temporary
- compressed_vector &operator -= (vector_container<C> const& v) {
- minus_assign(v);
- return *this;
- }
- template<class AE>
- compressed_vector &minus_assign(vector_expression<AE> const& ae) {
- kernels::assign<scalar_minus_assign> (*this, ae);
- return *this;
- }
- compressed_vector &operator *= (value_type t) {
- kernels::assign<scalar_multiply_assign> (*this, t);
- return *this;
- }
- compressed_vector &operator /= (value_type t) {
- kernels::assign<scalar_divide_assign> (*this, t);
+ swap(temporary);
return *this;
}
// Swapping
- void swap(compressed_vector &v) {
+ void swap(compressed_vector& v) {
std::swap(m_size, v.m_size);
std::swap(m_nnz, v.m_nnz);
m_indices.swap(v.m_indices);
@@ -367,16 +319,16 @@ public:
// Serialization
template<class Archive>
- void serialize(Archive &ar, const unsigned int /* file_version */) {
+ void serialize(Archive& ar, const unsigned int /* file_version */) {
boost::serialization::collection_size_type s(m_size);
- ar &boost::serialization::make_nvp("size",s);
+ ar & boost::serialization::make_nvp("size",s);
if (Archive::is_loading::value) {
m_size = s;
}
// ISSUE: m_indices and m_values are undefined between m_nnz and capacity (trouble with 'nan'-values)
- ar &boost::serialization::make_nvp("nnz", m_nnz);
- ar &boost::serialization::make_nvp("indices", m_indices);
- ar &boost::serialization::make_nvp("values", m_values);
+ ar & boost::serialization::make_nvp("nnz", m_nnz);
+ ar & boost::serialization::make_nvp("indices", m_indices);
+ ar & boost::serialization::make_nvp("values", m_values);
}
private:
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/shark.git
More information about the debian-science-commits
mailing list