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 REMORA_ASSIGNMENT_HPP
34#define REMORA_ASSIGNMENT_HPP
35
36#include "kernels/matrix_assign.hpp"
37#include "kernels/vector_assign.hpp"
38#include "detail/traits.hpp"
39#include "detail/proxy_optimizers_fwd.hpp"
40
41#include <type_traits>
42namespace remora{
43
44//////////////////////////////////////////////////////////////////////
45/////Evaluate blockwise expressions
46//////////////////////////////////////////////////////////////////////
47
48///\brief conditionally evaluates a vector expression if it is a block expression
49///
50/// If the expression is a block expression, a temporary vector is created to which
51/// the expression is assigned, which is then returned, otherwise the expression itself
52/// is returned
53template<class E, class Device>
54typename std::conditional<
55 std::is_base_of<
56 blockwise_tag,
57 typename E::evaluation_category
58 >::value,
59 typename vector_temporary<E>::type,
60 E const&
61>::type
62eval_block(vector_expression<E, Device> const& e){
63 return e();//either casts to E const& or returns the copied expression
64}
65///\brief conditionally evaluates a matrix expression if it is a block expression
66///
67/// If the expression is a block expression, a temporary matrix is created to which
68/// the expression is assigned, which is then returned, otherwise the expression itself
69/// is returned
70template<class E, class Device>
71typename std::conditional<
72 std::is_base_of<
73 blockwise_tag,
74 typename E::evaluation_category
75 >::value,
76 typename matrix_temporary<E>::type,
77 E const&
78>::type
79eval_block(matrix_expression<E, Device> const& e){
80 return e();//either casts to E const& or returns the copied expression
81}
82
83
84///\brief Evaluates an expression if it does not have a standard storage layout
85///
86/// This function evaluates an expression to a temporary if it does not have
87/// a known storage type. i.e. proxy expressions and containers are not evaluated but passed
88/// through while everything else is evaluated.
89template<class E, class Device>
90typename std::conditional<
91 std::is_same<
92 unknown_storage,
93 typename E::storage_type
94 >::value,
95 typename vector_temporary<E>::type,
96 E const&
97>::type
98eval_expression(vector_expression<E, Device> const& e){
99 return e();//either casts to E const& or returns the evaluated expression
100}
101///\brief Evaluates an expression if it does not have a standard storage layout
102///
103/// This function evaluates an expression to a temporary if it does not have
104/// a known storage type. i.e. proxy expressions and containers are not evaluated but passed
105/// through while everything else is evaluated.
106template<class E, class Device>
107typename std::conditional<
108 std::is_same<
109 unknown_storage,
110 typename E::storage_type
111 >::value,
112 typename matrix_temporary<E>::type,
113 E const&
114>::type
115eval_expression(matrix_expression<E, Device> const& e){
116 return e();//either casts to E const& or returns the evaluated expression
117}
118
119/////////////////////////////////////////////////////////////////////////////////////
120////// Vector Assign
121////////////////////////////////////////////////////////////////////////////////////
122
123namespace detail{
124 template<class VecX, class VecV, class Device>
125 void assign(
126 vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v,
127 elementwise_tag
128 ){
129 kernels::assign(x, v);
130 }
131 template<class VecX, class VecV, class Device>
132 void assign(
133 vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v,
134 blockwise_tag
135 ){
136 v().assign_to(x);
137 }
138 template<class VecX, class VecV, class Device>
139 void plus_assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v,
140 elementwise_tag
141 ){
142 typename device_traits<Device>:: template add<typename common_value_type<VecX,VecV>::type> f;
143 kernels::assign(x, v, f);
144 }
145 template<class VecX, class VecV, class Device>
146 void plus_assign(
147 vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v,
148 blockwise_tag
149 ){
150 v().plus_assign_to(x);
151 }
152}
153
154/// \brief Dispatches vector assignment on an expression level
155///
156/// This dispatcher takes care for whether the blockwise evaluation
157/// or the elementwise evaluation is called
158template<class VecX, class VecV, class Device>
159VecX& assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
160 REMORA_SIZE_CHECK(x().size() == v().size());
161 detail::assign(x,v,typename VecV::evaluation_category());
162 return x();
163}
164
165/// \brief Dispatches vector plus-assignment on an expression level
166///
167/// This dispatcher takes care for whether the blockwise evaluation
168/// or the elementwise evaluation is called
169template<class VecX, class VecV, class Device>
170VecX& plus_assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
171 REMORA_SIZE_CHECK(x().size() == v().size());
172 detail::plus_assign(x,v,typename VecV::evaluation_category());
173 return x();
174}
175
176/// \brief Dispatches vector minus-assignment on an expression level
177///
178/// This dispatcher takes care for whether the blockwise evaluation
179/// or the elementwise evaluation is called
180template<class VecX, class VecV, class Device>
181VecX& minus_assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
182 REMORA_SIZE_CHECK(x().size() == v().size());
183
184 typedef typename VecV::value_type value_type;
185 auto minusV = detail::vector_scalar_multiply_optimizer<VecV>::create(v(),value_type(-1));
186 detail::plus_assign(x,minusV,typename VecV::evaluation_category());
187 return x();
188}
189
190/// \brief Dispatches vector multiply-assignment on an expression level
191///
192/// This dispatcher takes care for whether the blockwise evaluation
193/// or the elementwise evaluation is called
194template<class VecX, class VecV, class Device>
195VecX& multiply_assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
196 REMORA_SIZE_CHECK(x().size() == v().size());
197 auto&& veval = eval_block(v);
198 typedef typename device_traits<Device>:: template multiply<typename common_value_type<VecX,VecV>::type> F;
199 kernels::assign(x, veval, F());
200 return x();
201}
202
203/// \brief Dispatches vector multiply-assignment on an expression level
204///
205/// This dispatcher takes care for whether the blockwise evaluation
206/// or the elementwise evaluation is called
207template<class VecX, class VecV, class Device>
208VecX& divide_assign(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
209 REMORA_SIZE_CHECK(x().size() == v().size());
210 auto&& veval = eval_block(v);
211 typedef typename device_traits<Device>:: template divide<typename common_value_type<VecX,VecV>::type> F;
212 kernels::assign(x, veval, F());
213 return x();
214}
215
216/////////////////////////////////////////////////////////////////////////////////////
217////// Matrix Assign
218////////////////////////////////////////////////////////////////////////////////////
219
220namespace detail{
221 template<class MatA, class MatB, class Device>
222 void assign(
223 matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B,
224 elementwise_tag
225 ){
226 kernels::assign(A, B);
227 }
228 template<class MatA, class MatB, class Device>
229 void assign(
230 matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B,
231 blockwise_tag
232 ){
233 B().assign_to(A);
234 }
235 template<class MatA, class MatB, class Device>
236 void plus_assign(
237 matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B,
238 elementwise_tag
239 ){
240 typename device_traits<Device>:: template add<typename common_value_type<MatA,MatB>::type> f;
241 kernels::assign(A,B,f);
242 }
243 template<class MatA, class MatB, class Device>
244 void plus_assign(
245 matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B,
246 blockwise_tag
247 ){
248 B().plus_assign_to(A);
249 }
250}
251
252/// \brief Dispatches matrix assignment on an expression level
253///
254/// This dispatcher takes care for whether the blockwise evaluation
255/// or the elementwise evaluation is called
256template<class MatA, class MatB, class Device>
257MatA& assign(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
258 REMORA_SIZE_CHECK(A().size1() == B().size1());
259 REMORA_SIZE_CHECK(A().size2() == B().size2());
260 detail::assign(A,B, typename MatB::evaluation_category());
261 return A();
262}
263
264/// \brief Dispatches matrix plus-assignment on an expression level
265///
266/// This dispatcher takes care for whether the blockwise evaluation
267/// or the elementwise evaluation is called
268template<class MatA, class MatB, class Device>
269MatA& plus_assign(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
270 REMORA_SIZE_CHECK(A().size1() == B().size1());
271 REMORA_SIZE_CHECK(A().size2() == B().size2());
272 detail::plus_assign(A,B, typename MatB::evaluation_category());
273 return A();
274}
275
276/// \brief Dispatches matrix minus-assignment on an expression level
277///
278/// This dispatcher takes care for whether the blockwise evaluation
279/// or the elementwise evaluation is called
280template<class MatA, class MatB, class Device>
281MatA& minus_assign(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
282 REMORA_SIZE_CHECK(A().size1() == B().size1());
283 REMORA_SIZE_CHECK(A().size2() == B().size2());
284 typedef typename MatB::value_type value_type;
285 auto minusB = detail::matrix_scalar_multiply_optimizer<MatB>::create(B(),value_type(-1));
286 detail::plus_assign(A, minusB, typename MatB::evaluation_category());
287 return A();
288}
289
290/// \brief Dispatches matrix multiply-assignment on an expression level
291///
292/// This dispatcher takes care for whether the blockwise evaluation
293/// or the elementwise evaluation is called
294template<class MatA, class MatB, class Device>
295MatA& multiply_assign(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
296 REMORA_SIZE_CHECK(A().size1() == B().size1());
297 REMORA_SIZE_CHECK(A().size2() == B().size2());
298 auto&& Beval = eval_block(B);
299 typedef typename device_traits<Device>:: template multiply<typename common_value_type<MatA,MatB>::type> F;
300 kernels::assign(A, Beval, F());
301 return A();
302}
303
304/// \brief Dispatches matrix divide-assignment on an expression level
305///
306/// This dispatcher takes care for whether the blockwise evaluation
307/// or the elementwise evaluation is called
308template<class MatA, class MatB, class Device>
309MatA& divide_assign(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
310 REMORA_SIZE_CHECK(A().size1() == B().size1());
311 REMORA_SIZE_CHECK(A().size2() == B().size2());
312 auto&& Beval = eval_block(B);
313 typedef typename device_traits<Device>:: template divide<typename common_value_type<MatA,MatB>::type> F;
314 kernels::assign(A, Beval, F());
315 return A();
316}
317
318//////////////////////////////////////////////////////////////////////////////////////
319///// Vector Operators
320/////////////////////////////////////////////////////////////////////////////////////
321
322/// \brief Add-Assigns two vector expressions
323///
324/// Performs the operation x_i+=v_i for all elements.
325/// Assumes that the right and left hand side aliases and therefore
326/// performs a copy of the right hand side before assigning
327/// use noalias as in noalias(x)+=v to avoid this if A and B do not alias
328template<class VecX, class VecV, class Device>
329VecX& operator+=(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
330 REMORA_SIZE_CHECK(x().size() == v().size());
331 typename vector_temporary<VecX>::type temporary(v);
332 return plus_assign(x,temporary);
333}
334
335template<class VecX, class VecV, class Device>
336typename VecX::closure_type operator+=(vector_expression<VecX, Device>&& x, vector_expression<VecV, Device> const& v){
337 REMORA_SIZE_CHECK(x().size() == v().size());
338 typename vector_temporary<VecX>::type temporary(v);
339 return plus_assign(x,temporary);
340}
341
342/// \brief Subtract-Assigns two vector expressions
343///
344/// Performs the operation x_i-=v_i for all elements.
345/// Assumes that the right and left hand side aliases and therefore
346/// performs a copy of the right hand side before assigning
347/// use noalias as in noalias(x)-=v to avoid this if A and B do not alias
348template<class VecX, class VecV, class Device>
349VecX& operator-=(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
350 REMORA_SIZE_CHECK(x().size() == v().size());
351 typename vector_temporary<VecX>::type temporary(v);
352 return minus_assign(x,temporary);
353}
354
355template<class VecX, class VecV, class Device>
356typename VecX::closure_type operator-=(vector_expression<VecX, Device>&& x, vector_expression<VecV, Device> const& v){
357 REMORA_SIZE_CHECK(x().size() == v().size());
358 typename vector_temporary<VecX>::type temporary(v);
359 return minus_assign(x,temporary);
360}
361
362/// \brief Multiply-Assigns two vector expressions
363///
364/// Performs the operation x_i*=v_i for all elements.
365/// Assumes that the right and left hand side aliases and therefore
366/// performs a copy of the right hand side before assigning
367/// use noalias as in noalias(x)*=v to avoid this if A and B do not alias
368template<class VecX, class VecV, class Device>
369VecX& operator*=(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
370 REMORA_SIZE_CHECK(x().size() == v().size());
371 typename vector_temporary<VecX>::type temporary(v);
372 return multiply_assign(x,temporary);
373}
374
375template<class VecX, class VecV, class Device>
376typename VecX::closure_type operator*=(vector_expression<VecX, Device>&& x, vector_expression<VecV, Device> const& v){
377 REMORA_SIZE_CHECK(x().size() == v().size());
378 typename vector_temporary<VecX>::type temporary(v);
379 multiply_assign(x,temporary);
380}
381
382/// \brief Divide-Assigns two vector expressions
383///
384/// Performs the operation x_i/=v_i for all elements.
385/// Assumes that the right and left hand side aliases and therefore
386/// performs a copy of the right hand side before assigning
387/// use noalias as in noalias(x)/=v to avoid this if A and B do not alias
388template<class VecX, class VecV, class Device>
389VecX& operator/=(vector_expression<VecX, Device>& x, vector_expression<VecV, Device> const& v){
390 REMORA_SIZE_CHECK(x().size() == v().size());
391 typename vector_temporary<VecX>::type temporary(v);
392 return divide_assign(x,temporary);
393}
394
395template<class VecX, class VecV, class Device>
396typename VecX::closure_type operator/=(vector_expression<VecX, Device>&& x, vector_expression<VecV, Device> const& v){
397 REMORA_SIZE_CHECK(x().size() == v().size());
398 typename vector_temporary<VecX>::type temporary(v);
399 divide_assign(x,temporary);
400}
401
402/// \brief Adds a scalar to all elements of the vector
403///
404/// Performs the operation x_i += t for all elements.
405template<class VecX, class T, class Device>
406typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,VecX&>::type
407operator+=(vector_expression<VecX, Device>& x, T t){
408 kernels::assign<typename device_traits<Device>:: template add<typename VecX::value_type> > (x, t);
409 return x();
410}
411
412template<class VecX, class T, class Device>
413typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,typename VecX::closure_type>::type
414operator+=(vector_expression<VecX, Device>&& x, T t){
415 kernels::assign<typename device_traits<Device>:: template add<typename VecX::value_type> > (x, t);
416 return x();
417}
418
419/// \brief Subtracts a scalar from all elements of the vector
420///
421/// Performs the operation x_i += t for all elements.
422template<class VecX, class T, class Device>
423typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,VecX&>::type
424operator-=(vector_expression<VecX, Device>& x, T t){
425 kernels::assign<typename device_traits<Device>:: template subtract<typename VecX::value_type> > (x, t);
426 return x();
427}
428
429template<class VecX, class T, class Device>
430typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,typename VecX::closure_type>::type
431operator-=(vector_expression<VecX, Device>&& x, T t){
432 kernels::assign<typename device_traits<Device>:: template subtract<typename VecX::value_type> > (x, t);
433 return x();
434}
435
436/// \brief Multiplies a scalar with all elements of the vector
437///
438/// Performs the operation x_i *= t for all elements.
439template<class VecX, class T, class Device>
440typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,VecX&>::type
441operator*=(vector_expression<VecX, Device>& x, T t){
442 kernels::assign<typename device_traits<Device>:: template multiply<typename VecX::value_type> > (x, t);
443 return x();
444}
445
446template<class VecX, class T, class Device>
447typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,typename VecX::closure_type>::type
448operator*=(vector_expression<VecX, Device>&& x, T t){
449 kernels::assign<typename device_traits<Device>:: template multiply<typename VecX::value_type> > (x, t);
450 return x();
451}
452
453/// \brief Divides all elements of the vector by a scalar
454///
455/// Performs the operation x_i /= t for all elements.
456template<class VecX, class T, class Device>
457typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,VecX&>::type
458operator/=(vector_expression<VecX, Device>& x, T t){
459 kernels::assign<typename device_traits<Device>:: template divide<typename VecX::value_type> > (x, t);
460 return x();
461}
462
463template<class VecX, class T, class Device>
464typename std::enable_if<std::is_convertible<T, typename VecX::value_type>::value,typename VecX::closure_type>::type
465operator/=(vector_expression<VecX, Device>&& x, T t){
466 kernels::assign<typename device_traits<Device>:: template divide<typename VecX::value_type> > (x, t);
467 return x();
468}
469
470
471
472//////////////////////////////////////////////////////////////////////////////////////
473///// Matrix Operators
474/////////////////////////////////////////////////////////////////////////////////////
475
476/// \brief Add-Assigns two matrix expressions
477///
478/// Performs the operation A_ij+=B_ij for all elements.
479/// Assumes that the right and left hand side aliases and therefore
480/// performs a copy of the right hand side before assigning
481/// use noalias as in noalias(A)+=B to avoid this if A and B do not alias
482template<class MatA, class MatB, class Device>
483MatA& operator+=(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
484 REMORA_SIZE_CHECK(A().size1() == B().size1());
485 REMORA_SIZE_CHECK(A().size2() == B().size2());
486 typename matrix_temporary<MatA>::type temporary(B);
487 return plus_assign(A,temporary);
488}
489
490template<class MatA, class MatB, class Device>
491typename MatA::closure_type operator+=(matrix_expression<MatA, Device>&& A, matrix_expression<MatB, Device> const& B){
492 REMORA_SIZE_CHECK(A().size1() == B().size1());
493 REMORA_SIZE_CHECK(A().size2() == B().size2());
494 typename matrix_temporary<MatA>::type temporary(B);
495 return plus_assign(A,temporary);
496}
497
498/// \brief Subtract-Assigns two matrix expressions
499///
500/// Performs the operation A_ij-=B_ij for all elements.
501/// Assumes that the right and left hand side aliases and therefore
502/// performs a copy of the right hand side before assigning
503/// use noalias as in noalias(A)-=B to avoid this if A and B do not alias
504template<class MatA, class MatB, class Device>
505MatA& operator-=(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
506 REMORA_SIZE_CHECK(A().size1() == B().size1());
507 REMORA_SIZE_CHECK(A().size2() == B().size2());
508 typename matrix_temporary<MatA>::type temporary(B);
509 return minus_assign(A,temporary);
510}
511
512template<class MatA, class MatB, class Device>
513typename MatA::closure_type operator-=(matrix_expression<MatA, Device>&& A, matrix_expression<MatB, Device> const& B){
514 REMORA_SIZE_CHECK(A().size1() == B().size1());
515 REMORA_SIZE_CHECK(A().size2() == B().size2());
516 typename matrix_temporary<MatA>::type temporary(B);
517 return minus_assign(A,temporary);
518}
519
520/// \brief Multiply-Assigns two matrix expressions
521///
522/// Performs the operation A_ij*=B_ij for all elements.
523/// Assumes that the right and left hand side aliases and therefore
524/// performs a copy of the right hand side before assigning
525/// use noalias as in noalias(A)*=B to avoid this if A and B do not alias
526template<class MatA, class MatB, class Device>
527MatA& operator*=(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
528 REMORA_SIZE_CHECK(A().size1() == B().size1());
529 REMORA_SIZE_CHECK(A().size2() == B().size2());
530 typename matrix_temporary<MatA>::type temporary(B);
531 return multiply_assign(A,temporary);
532}
533
534template<class MatA, class MatB, class Device>
535typename MatA::closure_type operator*=(matrix_expression<MatA, Device>&& A, matrix_expression<MatB, Device> const& B){
536 REMORA_SIZE_CHECK(A().size1() == B().size1());
537 REMORA_SIZE_CHECK(A().size2() == B().size2());
538 typename matrix_temporary<MatA>::type temporary(B);
539 return multiply_assign(A,temporary);
540}
541
542/// \brief Divide-Assigns two matrix expressions
543///
544/// Performs the operation A_ij/=B_ij for all elements.
545/// Assumes that the right and left hand side aliases and therefore
546/// performs a copy of the right hand side before assigning
547/// use noalias as in noalias(A)/=B to avoid this if A and B do not alias
548template<class MatA, class MatB, class Device>
549MatA& operator/=(matrix_expression<MatA, Device>& A, matrix_expression<MatB, Device> const& B){
550 REMORA_SIZE_CHECK(A().size1() == B().size1());
551 REMORA_SIZE_CHECK(A().size2() == B().size2());
552 typename matrix_temporary<MatA>::type temporary(B);
553 return divide_assign(A,temporary);
554}
555
556template<class MatA, class MatB, class Device>
557typename MatA::closure_type operator/=(matrix_expression<MatA, Device>&& A, matrix_expression<MatB, Device> const& B){
558 REMORA_SIZE_CHECK(A().size1() == B().size1());
559 REMORA_SIZE_CHECK(A().size2() == B().size2());
560 typename matrix_temporary<MatA>::type temporary(B);
561 return divide_assign(A,temporary);
562}
563
564/// \brief Adds a scalar to all elements of the matrix
565///
566/// Performs the operation A_ij += t for all elements.
567template<class MatA, class T, class Device>
568typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,MatA&>::type
569operator+=(matrix_expression<MatA, Device>& A, T t){
570 kernels::assign<typename device_traits<Device>:: template add<typename MatA::value_type> > (A, t);
571 return A();
572}
573
574template<class MatA, class T, class Device>
575typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,typename MatA::closure_type>::type
576operator+=(matrix_expression<MatA, Device>&& A, T t){
577 kernels::assign<typename device_traits<Device>:: template add<typename MatA::value_type> > (A, t);
578 return A();
579}
580
581/// \brief Subtracts a scalar from all elements of the matrix
582///
583/// Performs the operation A_ij -= t for all elements.
584template<class MatA, class T, class Device>
585typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,MatA&>::type
586operator-=(matrix_expression<MatA, Device>& A, T t){
587 kernels::assign<typename device_traits<Device>:: template subtract<typename MatA::value_type> > (A, t);
588 return A();
589}
590
591template<class MatA, class T, class Device>
592typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,typename MatA::closure_type>::type
593operator-=(matrix_expression<MatA, Device>&& A, T t){
594 kernels::assign<typename device_traits<Device>:: template subtract<typename MatA::value_type> > (A, t);
595 return A();
596}
597
598/// \brief Multiplies a scalar to all elements of the matrix
599///
600/// Performs the operation A_ij *= t for all elements.
601template<class MatA, class T, class Device>
602typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,MatA&>::type
603operator*=(matrix_expression<MatA, Device>& A, T t){
604 kernels::assign<typename device_traits<Device>:: template multiply<typename MatA::value_type> > (A, t);
605 return A();
606}
607
608template<class MatA, class T, class Device>
609typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,typename MatA::closure_type>::type
610operator*=(matrix_expression<MatA, Device>&& A, T t){
611 kernels::assign<typename device_traits<Device>:: template multiply<typename MatA::value_type> > (A, t);
612 return A();
613}
614
615/// \brief Divides all elements of the matrix by a scalar
616///
617/// Performs the operation A_ij /= t for all elements.
618template<class MatA, class T, class Device>
619typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,MatA&>::type
620operator/=(matrix_expression<MatA, Device>& A, T t){
621 kernels::assign<typename device_traits<Device>:: template divide<typename MatA::value_type> > (A, t);
622 return A();
623}
624
625template<class MatA, class T, class Device>
626typename std::enable_if<std::is_convertible<T, typename MatA::value_type>::value,typename MatA::closure_type>::type
627operator/=(matrix_expression<MatA, Device>&& A, T t){
628 kernels::assign<typename device_traits<Device>:: template divide<typename MatA::value_type> > (A, t);
629 return A();
630}
631
632// Assignment proxy.
633// Provides temporary free assigment when LHS has no alias on RHS
634template<class C>
635class noalias_proxy{
636public:
637 typedef typename C::closure_type closure_type;
638 typedef typename C::value_type value_type;
639
640 noalias_proxy(C &lval): m_lval(lval) {}
641
642 noalias_proxy(const noalias_proxy &p):m_lval(p.m_lval) {}
643
644 template <class E>
645 closure_type &operator= (const E &e) {
646 return assign(m_lval, e);
647 }
648
649 template <class E>
650 closure_type &operator+= (const E &e) {
651 return plus_assign(m_lval, e);
652 }
653
654 template <class E>
655 closure_type &operator-= (const E &e) {
656 return minus_assign(m_lval, e);
657 }
658
659 template <class E>
660 closure_type &operator*= (const E &e) {
661 return multiply_assign(m_lval, e);
662 }
663
664 template <class E>
665 closure_type &operator/= (const E &e) {
666 return divide_assign(m_lval, e);
667 }
668
669 //this is not needed, but prevents errors when for example doing noalias(x)+=2;
670 closure_type &operator+= (value_type t) {
671 return m_lval += t;
672 }
673
674 //this is not needed, but prevents errors when for example doing noalias(x)-=2;
675 closure_type &operator-= (value_type t) {
676 return m_lval -= t;
677 }
678
679 //this is not needed, but prevents errors when for example doing noalias(x)*=2;
680 closure_type &operator*= (value_type t) {
681 return m_lval *= t;
682 }
683
684 //this is not needed, but prevents errors when for example doing noalias(x)/=2;
685 closure_type &operator/= (value_type t) {
686 return m_lval /= t;
687 }
688
689private:
690 closure_type m_lval;
691};
692
693// Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
694// noalias(lhs) = rhs_expression
695template <class C, class Device>
696noalias_proxy<C> noalias(matrix_expression<C, Device>& lvalue) {
697 return noalias_proxy<C> (lvalue());
698}
699template <class C, class Device>
700noalias_proxy<C> noalias(vector_expression<C, Device>& lvalue) {
701 return noalias_proxy<C> (lvalue());
702}
703template <class C, class Device>
704noalias_proxy<C> noalias(vector_expression<C, Device>&& rvalue) {
705 return noalias_proxy<C>(rvalue());
706}
707template <class C, class Device>
708noalias_proxy<C> noalias(matrix_expression<C, Device>&& lvalue) {
709 return noalias_proxy<C> (lvalue());
710}
711
712template <class C, class Device>
713noalias_proxy<C> noalias(vector_set_expression<C, Device>& lvalue) {
714 return noalias_proxy<C> (lvalue());
715}
716
717
718}
719#endif