assignment.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Assignment and evaluation of vector expressions
5  *
6  *
7  *
8  * \author O. Krause
9  * \date 2013
10  *
11  *
12  * \par Copyright 1995-2015 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://image.diku.dk/shark/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 
33 #ifndef SHARK_LINALG_BLAS_ASSIGNMENT_HPP
34 #define SHARK_LINALG_BLAS_ASSIGNMENT_HPP
35 
38 #include "detail/traits.hpp"
39 
40 namespace shark {
41 namespace blas {
42 
43 /////////////////////////////////////////////////////////////////////////////////////
44 ////// Vector Assign
45 ////////////////////////////////////////////////////////////////////////////////////
46 
47 namespace detail{
48  template<class VecX, class VecV>
49  void assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
50  kernels::assign(x,v);
51  }
52  template<class VecX, class VecV>
53  void assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
54  v().assign_to(x);
55  }
56  template<class VecX, class VecV>
57  void plus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
58  kernels::assign<scalar_plus_assign> (x, v);
59  }
60  template<class VecX, class VecV>
61  void plus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
62  v().plus_assign_to(x);
63  }
64  template<class VecX, class VecV>
65  void minus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
66  kernels::assign<scalar_minus_assign> (x, v);
67  }
68  template<class VecX, class VecV>
69  void minus_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
70  v().minus_assign_to(x);
71  }
72  template<class VecX, class VecV>
73  void multiply_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
74  kernels::assign<scalar_multiply_assign> (x, v);
75  }
76  template<class VecX, class VecV>
77  void multiply_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
78  typename vector_temporary<VecX>::type temporary(v);
79  kernels::assign<scalar_multiply_assign> (x, temporary);
80  }
81  template<class VecX, class VecV>
82  void divide_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,elementwise_tag){
83  kernels::assign<scalar_divide_assign> (x, v);
84  }
85  template<class VecX, class VecV>
86  void divide_assign(vector_expression<VecX>& x, vector_expression<VecV> const& v,blockwise_tag){
87  typename vector_temporary<VecX>::type temporary(v);
88  kernels::assign<scalar_divide_assign> (x, temporary);
89  }
90 }
91 
92 
93 /// \brief Dispatches vector assignment on an expression level
94 ///
95 /// This dispatcher takes care for whether the blockwise evaluation
96 /// or the elementwise evaluation is called
97 template<class VecX, class VecV>
99  SIZE_CHECK(x().size() == v().size());
100  detail::assign(x,v,typename VecV::evaluation_category());
101  return x();
102 }
103 
104 /// \brief Dispatches vector plus-assignment on an expression level
105 ///
106 /// This dispatcher takes care for whether the blockwise evaluation
107 /// or the elementwise evaluation is called
108 template<class VecX, class VecV>
110  SIZE_CHECK(x().size() == v().size());
111  detail::plus_assign(x,v,typename VecV::evaluation_category());
112  return x();
113 }
114 
115 /// \brief Dispatches vector minus-assignment on an expression level
116 ///
117 /// This dispatcher takes care for whether the blockwise evaluation
118 /// or the elementwise evaluation is called
119 template<class VecX, class VecV>
121  SIZE_CHECK(x().size() == v().size());
122  detail::minus_assign(x,v,typename VecV::evaluation_category());
123  return x();
124 }
125 
126 /// \brief Dispatches vector multiply-assignment on an expression level
127 ///
128 /// This dispatcher takes care for whether the blockwise evaluation
129 /// or the elementwise evaluation is called
130 template<class VecX, class VecV>
132  SIZE_CHECK(x().size() == v().size());
133  detail::multiply_assign(x,v,typename VecV::evaluation_category());
134  return x();
135 }
136 
137 /// \brief Dispatches vector multiply-assignment on an expression level
138 ///
139 /// This dispatcher takes care for whether the blockwise evaluation
140 /// or the elementwise evaluation is called
141 template<class VecX, class VecV>
143  SIZE_CHECK(x().size() == v().size());
144  detail::divide_assign(x,v,typename VecV::evaluation_category());
145  return x();
146 }
147 
148 /////////////////////////////////////////////////////////////////////////////////////
149 ////// Matrix Assign
150 ////////////////////////////////////////////////////////////////////////////////////
151 
152 namespace detail{
153  template<class MatA, class MatB>
154  void assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
155  kernels::assign(A,B);
156  }
157  template<class MatA, class MatB>
158  void assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
159  B().assign_to(A);
160  }
161  template<class MatA, class MatB>
162  void plus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
163  kernels::assign<scalar_plus_assign> (A, B);
164  }
165  template<class MatA, class MatB>
166  void plus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
167  B().plus_assign_to(A);
168  }
169  template<class MatA, class MatB>
170  void minus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
171  kernels::assign<scalar_minus_assign> (A, B);
172  }
173  template<class MatA, class MatB>
174  void minus_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
175  B().minus_assign_to(A);
176  }
177  template<class MatA, class MatB>
178  void multiply_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
179  kernels::assign<scalar_multiply_assign> (A, B);
180  }
181  template<class MatA, class MatB>
182  void multiply_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
183  typename matrix_temporary<MatA>::type temporary(B);
184  kernels::assign<scalar_multiply_assign> (A, B);
185  }
186  template<class MatA, class MatB>
187  void divide_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,elementwise_tag){
188  kernels::assign<scalar_divide_assign> (A, B);
189  }
190  template<class MatA, class MatB>
191  void divide_assign(matrix_expression<MatA>& A, matrix_expression<MatB> const& B,blockwise_tag){
192  typename matrix_temporary<MatA>::type temporary(B);
193  kernels::assign<scalar_divide_assign> (A, B);
194  }
195 }
196 
197 
198 /// \brief Dispatches matrix assignment on an expression level
199 ///
200 /// This dispatcher takes care for whether the blockwise evaluation
201 /// or the elementwise evaluation is called
202 template<class MatA, class MatB>
204  SIZE_CHECK(A().size1() == B().size1());
205  SIZE_CHECK(A().size2() == B().size2());
206  detail::assign(A,B, typename MatB::evaluation_category());
207  return A();
208 }
209 
210 /// \brief Dispatches matrix plus-assignment on an expression level
211 ///
212 /// This dispatcher takes care for whether the blockwise evaluation
213 /// or the elementwise evaluation is called
214 template<class MatA, class MatB>
216  SIZE_CHECK(A().size1() == B().size1());
217  SIZE_CHECK(A().size2() == B().size2());
218  detail::plus_assign(A,B, typename MatB::evaluation_category());
219  return A();
220 }
221 
222 /// \brief Dispatches matrix plus-assignment on an expression level
223 ///
224 /// This dispatcher takes care for whether the blockwise evaluation
225 /// or the elementwise evaluation is called
226 template<class MatA, class MatB>
228  SIZE_CHECK(A().size1() == B().size1());
229  SIZE_CHECK(A().size2() == B().size2());
230  detail::minus_assign(A,B, typename MatB::evaluation_category());
231  return A();
232 }
233 
234 /// \brief Dispatches matrix multiply-assignment on an expression level
235 ///
236 /// This dispatcher takes care for whether the blockwise evaluation
237 /// or the elementwise evaluation is called
238 template<class MatA, class MatB>
240  SIZE_CHECK(A().size1() == B().size1());
241  SIZE_CHECK(A().size2() == B().size2());
242  detail::multiply_assign(A,B, typename MatB::evaluation_category());
243  return A();
244 }
245 
246 /// \brief Dispatches matrix divide-assignment on an expression level
247 ///
248 /// This dispatcher takes care for whether the blockwise evaluation
249 /// or the elementwise evaluation is called
250 template<class MatA, class MatB>
252  SIZE_CHECK(A().size1() == B().size1());
253  SIZE_CHECK(A().size2() == B().size2());
254  detail::divide_assign(A,B, typename MatB::evaluation_category());
255  return A();
256 }
257 
258 //////////////////////////////////////////////////////////////////////////////////////
259 ///// Vector Operators
260 /////////////////////////////////////////////////////////////////////////////////////
261 
262 /// \brief Add-Assigns two vector expressions
263 ///
264 /// Performs the operation x_i+=v_i for all elements.
265 /// Assumes that the right and left hand side aliases and therefore
266 /// performs a copy of the right hand side before assigning
267 /// use noalias as in noalias(x)+=v to avoid this if A and B do not alias
268 template<class VecX, class VecV>
270  SIZE_CHECK(x().size() == v().size());
271  typename vector_temporary<VecX>::type temporary(v);
272  return plus_assign(x,temporary);
273 }
274 
275 /// \brief Subtract-Assigns two vector expressions
276 ///
277 /// Performs the operation x_i-=v_i for all elements.
278 /// Assumes that the right and left hand side aliases and therefore
279 /// performs a copy of the right hand side before assigning
280 /// use noalias as in noalias(x)-=v to avoid this if A and B do not alias
281 template<class VecX, class VecV>
283  SIZE_CHECK(x().size() == v().size());
284  typename vector_temporary<VecX>::type temporary(v);
285  return minus_assign(x,temporary);
286 }
287 
288 /// \brief Multiply-Assigns two vector expressions
289 ///
290 /// Performs the operation x_i*=v_i for all elements.
291 /// Assumes that the right and left hand side aliases and therefore
292 /// performs a copy of the right hand side before assigning
293 /// use noalias as in noalias(x)*=v to avoid this if A and B do not alias
294 template<class VecX, class VecV>
296  SIZE_CHECK(x().size() == v().size());
297  typename vector_temporary<VecX>::type temporary(v);
298  return multiply_assign(x,temporary);
299 }
300 
301 /// \brief Divide-Assigns two vector expressions
302 ///
303 /// Performs the operation x_i/=v_i for all elements.
304 /// Assumes that the right and left hand side aliases and therefore
305 /// performs a copy of the right hand side before assigning
306 /// use noalias as in noalias(x)/=v to avoid this if A and B do not alias
307 template<class VecX, class VecV>
309  SIZE_CHECK(x().size() == v().size());
310  typename vector_temporary<VecX>::type temporary(v);
311  return divide_assign(x,temporary);
312 }
313 
314 /// \brief Adds a scalar to all elements of the vector
315 ///
316 /// Performs the operation x_i += t for all elements.
317 template<class VecX>
318 VecX& operator+=(vector_expression<VecX>& x, typename VecX::scalar_type t){
319  kernels::assign<scalar_plus_assign> (x, t);
320  return x();
321 }
322 
323 /// \brief Subtracts a scalar from all elements of the vector
324 ///
325 /// Performs the operation x_i += t for all elements.
326 template<class VecX>
327 VecX& operator-=(vector_expression<VecX>& x, typename VecX::scalar_type t){
328  kernels::assign<scalar_minus_assign> (x, t);
329  return x();
330 }
331 
332 /// \brief Multiplies a scalar with all elements of the vector
333 ///
334 /// Performs the operation x_i *= t for all elements.
335 template<class VecX>
336 VecX& operator*=(vector_expression<VecX>& x, typename VecX::scalar_type t){
337  kernels::assign<scalar_multiply_assign> (x, t);
338  return x();
339 }
340 
341 /// \brief Divides all elements of the vector by a scalar
342 ///
343 /// Performs the operation x_i /= t for all elements.
344 template<class VecX>
345 VecX& operator/=(vector_expression<VecX>& x, typename VecX::scalar_type t){
346  kernels::assign<scalar_divide_assign> (x, t);
347  return x();
348 }
349 
350 
351 
352 //////////////////////////////////////////////////////////////////////////////////////
353 ///// Matrix Operators
354 /////////////////////////////////////////////////////////////////////////////////////
355 
356 /// \brief Add-Assigns two matrix expressions
357 ///
358 /// Performs the operation A_ij+=B_ij for all elements.
359 /// Assumes that the right and left hand side aliases and therefore
360 /// performs a copy of the right hand side before assigning
361 /// use noalias as in noalias(A)+=B to avoid this if A and B do not alias
362 template<class MatA, class MatB>
364  SIZE_CHECK(A().size1() == B().size1());
365  SIZE_CHECK(A().size2() == B().size2());
366  typename matrix_temporary<MatA>::type temporary(B);
367  return plus_assign(A,temporary);
368 }
369 
370 /// \brief Subtract-Assigns two matrix expressions
371 ///
372 /// Performs the operation A_ij-=B_ij for all elements.
373 /// Assumes that the right and left hand side aliases and therefore
374 /// performs a copy of the right hand side before assigning
375 /// use noalias as in noalias(A)-=B to avoid this if A and B do not alias
376 template<class MatA, class MatB>
378  SIZE_CHECK(A().size1() == B().size1());
379  SIZE_CHECK(A().size2() == B().size2());
380  typename matrix_temporary<MatA>::type temporary(B);
381  return minus_assign(A,temporary);
382 }
383 
384 /// \brief Multiply-Assigns two matrix expressions
385 ///
386 /// Performs the operation A_ij*=B_ij for all elements.
387 /// Assumes that the right and left hand side aliases and therefore
388 /// performs a copy of the right hand side before assigning
389 /// use noalias as in noalias(A)*=B to avoid this if A and B do not alias
390 template<class MatA, class MatB>
392  SIZE_CHECK(A().size1() == B().size1());
393  SIZE_CHECK(A().size2() == B().size2());
394  typename matrix_temporary<MatA>::type temporary(B);
395  return multiply_assign(A,temporary);
396 }
397 
398 /// \brief Divide-Assigns two matrix expressions
399 ///
400 /// Performs the operation A_ij/=B_ij for all elements.
401 /// Assumes that the right and left hand side aliases and therefore
402 /// performs a copy of the right hand side before assigning
403 /// use noalias as in noalias(A)/=B to avoid this if A and B do not alias
404 template<class MatA, class MatB>
406  SIZE_CHECK(A().size1() == B().size1());
407  SIZE_CHECK(A().size2() == B().size2());
408  typename matrix_temporary<MatA>::type temporary(B);
409  return divide_assign(A,temporary);
410 }
411 
412 /// \brief Adds a scalar to all elements of the matrix
413 ///
414 /// Performs the operation A_ij += t for all elements.
415 template<class MatA>
416 MatA& operator+=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
417  kernels::assign<scalar_plus_assign> (A, t);
418  return A();
419 }
420 
421 /// \brief Subtracts a scalar from all elements of the matrix
422 ///
423 /// Performs the operation A_ij -= t for all elements.
424 template<class MatA>
425 MatA& operator-=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
426  kernels::assign<scalar_minus_assign> (A, t);
427  return A();
428 }
429 
430 /// \brief Multiplies a scalar to all elements of the matrix
431 ///
432 /// Performs the operation A_ij *= t for all elements.
433 template<class MatA>
434 MatA& operator*=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
435  kernels::assign<scalar_multiply_assign> (A, t);
436  return A();
437 }
438 
439 /// \brief Divides all elements of the matrix by a scalar
440 ///
441 /// Performs the operation A_ij /= t for all elements.
442 template<class MatA>
443 MatA& operator /=(matrix_expression<MatA>& A, typename MatA::scalar_type t){
444  kernels::assign<scalar_divide_assign> (A, t);
445  return A();
446 }
447 
448 
449 
450 //////////////////////////////////////////////////////////////////////////////////////
451 ///// Temporary Proxy Operators
452 /////////////////////////////////////////////////////////////////////////////////////
453 
454 template<class T, class U>
456  static_cast<T&>(x) += arg;
457  return x;
458 }
459 template<class T, class U>
461  static_cast<T&>(x) -= arg;
462  return x;
463 }
464 template<class T, class U>
466  static_cast<T&>(x) *= arg;
467  return x;
468 }
469 template<class T, class U>
471  static_cast<T&>(x) /= arg;
472  return x;
473 }
474 
475 
476 
477 
478 // Assignment proxy.
479 // Provides temporary free assigment when LHS has no alias on RHS
480 template<class C>
482 public:
483  typedef typename C::closure_type closure_type;
484  typedef typename C::scalar_type scalar_type;
485 
486  noalias_proxy(C &lval): m_lval(lval) {}
487 
488  noalias_proxy(const noalias_proxy &p):m_lval(p.m_lval) {}
489 
490  template <class E>
491  closure_type &operator= (const E &e) {
492  return assign(m_lval, e);
493  }
494 
495  template <class E>
496  closure_type &operator+= (const E &e) {
497  return plus_assign(m_lval, e);
498  }
499 
500  template <class E>
501  closure_type &operator-= (const E &e) {
502  return minus_assign(m_lval, e);
503  }
504 
505  template <class E>
506  closure_type &operator*= (const E &e) {
507  return multiply_assign(m_lval, e);
508  }
509 
510  template <class E>
511  closure_type &operator/= (const E &e) {
512  return divide_assign(m_lval, e);
513  }
514 
515  //this is not needed, but prevents errors when for example doing noalias(x)+=2;
516  closure_type &operator+= (scalar_type t) {
517  return m_lval += t;
518  }
519 
520  //this is not needed, but prevents errors when for example doing noalias(x)-=2;
521  closure_type &operator-= (scalar_type t) {
522  return m_lval -= t;
523  }
524 
525  //this is not needed, but prevents errors when for example doing noalias(x)*=2;
526  closure_type &operator*= (scalar_type t) {
527  return m_lval *= t;
528  }
529 
530  //this is not needed, but prevents errors when for example doing noalias(x)/=2;
531  closure_type &operator/= (scalar_type t) {
532  return m_lval /= t;
533  }
534 
535 private:
536  closure_type m_lval;
537 };
538 
539 // Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
540 // noalias(lhs) = rhs_expression
541 template <class C>
543  return noalias_proxy<C> (lvalue());
544 }
545 template <class C>
547  return noalias_proxy<C> (lvalue());
548 }
549 
550 template <class C>
552  return noalias_proxy<C> (lvalue());
553 }
554 template <class C>
556  return noalias_proxy<C> (lvalue());
557 }
558 template <class C>
560  return noalias_proxy<C> (static_cast<C&>(lvalue));
561 }
562 
563 
564 
565 
566 //////////////////////////////////////////////////////////////////////
567 /////Evaluate blockwise expressions
568 //////////////////////////////////////////////////////////////////////
569 
570 ///\brief conditionally evaluates a vector expression if it is a block expression
571 ///
572 /// If the expression is a block expression, a temporary vector is created to which
573 /// the expression is assigned, which is then returned, otherwise the expression itself
574 /// is returned
575 template<class E>
576 typename boost::mpl::eval_if<
577  boost::is_same<
578  typename E::evaluation_category,
579  blockwise_tag
580  >,
581  vector_temporary<E>,
582  boost::mpl::identity<E const&>
583 >::type
585  return e();//either casts to E const& or returns the copied expression
586 }
587 ///\brief conditionally evaluates a matrix expression if it is a block expression
588 ///
589 /// If the expression is a block expression, a temporary matrix is created to which
590 /// the expression is assigned, which is then returned, otherwise the expression itself
591 /// is returned
592 template<class E>
593 typename boost::mpl::eval_if<
594  boost::is_same<
595  typename E::evaluation_category,
596  blockwise_tag
597  >,
598  matrix_temporary<E>,
599  boost::mpl::identity<E const&>
600 >::type
602  return e();//either casts to E const& or returns the copied expression
603 }
604 
605 }}
606 #endif