Rheolef  7.1
an efficient C++ finite element environment
vec_expr_v2.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_VEC_EXPR_v2_H
2 #define _RHEOLEF_VEC_EXPR_v2_H
23 // vec-valued affine expressions:
24 //
25 // vec w = 2*u - v + 1:
26 //
27 // author: Pierre.Saramito@imag.fr
28 //
29 // date: 14 september 2015
30 //
31 // Notes; implementation uses template expressions and SFINAE techiques
32 // * template expressions allows one to eliminate temporaries
33 // * SFINAE reduces the code size
34 //
35 // OVERVIEW:
36 // 1. linear/affine algebra operators
37 // 1.1. vec_expr, a type concept for expression types
38 // 1.2. unary operations
39 // 1.3. binary operations
40 // 2. completing the interface of the vec<T,M> class
41 // 2.1. cstor & assignment
42 // 2.2. computed assignment
43 // 2.3. scalar product
44 //
45 #include "rheolef/vec.h"
46 #include "rheolef/dis_inner_product.h"
47 #include "rheolef/dis_accumulate.h"
48 #include "rheolef/expr_utilities.h"
49 
50 namespace rheolef {
51 
52 // ===================================================================
53 // 1. linear/affine algebra operators
54 // ===================================================================
55 // -------------------------------------------------------------------
56 // 1.1. vec_expr, a type concept for expression types
57 // -------------------------------------------------------------------
58 namespace details {
59 
60 // Define a trait type for detecting vec expression valid arguments
61 template<class T> struct is_vec : std::false_type {};
62 template<class T, class M> struct is_vec<vec<T, M> > : std::true_type {};
63 
64 // Define a trait type for detecting vec expression valid arguments
65 template<class T> struct is_vec_expr_v2_arg : std::false_type {};
66 template<class T, class M> struct is_vec_expr_v2_arg <vec<T, M> > : std::true_type {};
67 
68 } // namespace details
69 // -------------------------------------------
70 // 1.2. unary operations
71 // -------------------------------------------
72 namespace details {
73 
74 template <class Op, class Expr>
76 
77 // typedefs:
78 
79  typedef typename Expr::size_type size_type;
80  typedef typename Expr::value_type value_type;
81  typedef typename Expr::memory_type memory_type;
82  typedef typename Expr::const_iterator expr_const_iterator;
84 
85 // allocatos:
86 
87  vec_expr_v2_unary (const Op& op, const Expr& expr)
88  : _op(op), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}
89 
90  template <class BinaryOp, class Constant>
91  vec_expr_v2_unary (const BinaryOp& binop, const Constant& c, const Expr& expr)
92  : _op(binop,c), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}
93 
94  template <class BinaryOp, class Constant>
95  vec_expr_v2_unary (const BinaryOp& binop, const Expr& expr, const Constant& c)
96  : _op(binop,c), _expr_iter(expr.begin()), _ownership(expr.ownership()) {}
97 
98 // accessors:
99 
100  const distributor& ownership() const { return _ownership; }
101  size_type size() const { return _ownership.size(); }
102 
103 // minimal forward iterator interface:
104 
105  struct const_iterator {
106  typedef std::forward_iterator_tag iterator_category;
107  typedef typename Expr::value_type value_type;
109  typedef value_type* pointer;
110  typedef std::ptrdiff_t difference_type;
112  : _op(op), _expr_iter (expr_iter) {}
113  const_iterator& operator++ () { ++_expr_iter; return *this; }
114  value_type operator* () const { return _op (*_expr_iter); }
115  protected:
116  const Op _op;
118  };
120 protected:
121  const Op _op;
124 };
125 template<class Op, class Expr> struct is_vec_expr_v2_arg <vec_expr_v2_unary<Op,Expr> > : std::true_type {};
126 
127 } // namespace details
128 
129 #define _RHEOLEF_vec_expr_v2_unary_operator(OP, FUNCTOR) \
130 template <class Expr> \
131 inline \
132 typename \
133 std::enable_if< \
134  details::is_vec_expr_v2_arg<Expr>::value, \
135  details::vec_expr_v2_unary< \
136  FUNCTOR, \
137  Expr \
138  > \
139 >::type \
140 operator OP (const Expr& expr) \
141 { \
142  typedef details::vec_expr_v2_unary <FUNCTOR, Expr> expr_t; \
143  return expr_t (FUNCTOR(), expr); \
144 }
145 
146 _RHEOLEF_vec_expr_v2_unary_operator(+, details::generic_unary_plus<>)
147 _RHEOLEF_vec_expr_v2_unary_operator(-, details::generic_negate<>)
148 #undef _RHEOLEF_vec_expr_v2_unary_operator
149 
150 // -------------------------------------------
151 // 1.3. binary operations
152 // -------------------------------------------
153 namespace details {
154 
155 template <class Op, class Expr1, class Expr2>
157 
158  typedef typename Expr1::size_type size_type;
159  typedef typename promote<
160  typename Expr1::value_type,
161  typename Expr2::value_type>::type value_type;
162  typedef typename Expr1::memory_type memory_type; // TODO: check Expr2::memory_type
163  typedef typename Expr1::const_iterator expr1_const_iterator;
164  typedef typename Expr2::const_iterator expr2_const_iterator;
166 
167 // allocators:
168 
169  vec_expr_v2_binary (const Op& op, const Expr1& expr1, const Expr2& expr2)
170  : _op (op),
171  _iter1 (expr1.begin()),
172  _iter2 (expr2.begin()),
173  _ownership (expr1.ownership())
174  {
175  check_macro (expr1.size() == expr2.size(), "linear binary vec expression: incompatible sizes "
176  << expr1.size() << " and " << expr2.size());
177  }
178 
179 // accessors:
180 
181  const distributor& ownership() const { return _ownership; }
182  size_type size() const { return _ownership.size(); }
183 
184 // minimal forward iterator interface:
185 
186  struct const_iterator {
187  typedef std::forward_iterator_tag iterator_category;
188  typedef typename promote<
189  typename Expr1::value_type,
190  typename Expr2::value_type>::type value_type;
192  typedef value_type* pointer;
193  typedef std::ptrdiff_t difference_type;
195  : _op(op), _iter1 (iter1), _iter2 (iter2) {}
196  const_iterator& operator++ () { ++_iter1; ++_iter2; return *this; }
197  value_type operator* () const { return _op (*_iter1, *_iter2); }
198  protected:
199  const Op _op;
202  };
204 protected:
205  const Op _op;
209 };
210 template<class Op, class Expr1, class Expr2> struct is_vec_expr_v2_arg <vec_expr_v2_binary<Op,Expr1,Expr2> > : std::true_type {};
211 
212 template <class Op, class Expr1, class Expr2, class Sfinae = void>
213 struct vec_expr_v2_binary_traits { /* catch-all case */ };
214 
215 // vec_expr +- vec_expr
216 template <class Op, class Expr1, class Expr2>
217 struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
218  typename std::enable_if<
219  details::is_vec_expr_v2_arg<Expr1>::value &&
220  details::is_vec_expr_v2_arg<Expr2>::value>::type>
221 {
223 };
224 // constant +-* vec_expr
225 template <class Op, class Expr1, class Expr2>
226 struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
227  typename std::enable_if<
228  details::is_rheolef_arithmetic<Expr1>::value &&
229  details::is_vec_expr_v2_arg<Expr2>::value>::type>
230 {
233 };
234 // vec_expr +-* constant
235 template <class Op, class Expr1, class Expr2>
236 struct vec_expr_v2_binary_traits <Op, Expr1, Expr2,
237  typename std::enable_if<
238  details::is_vec_expr_v2_arg<Expr1>::value &&
239  details::is_rheolef_arithmetic<Expr2>::value>::type>
240 {
243 };
244 
245 } // namespace details
246 
247 // x+y; x+c ; c+x
248 #define _RHEOLEF_vec_expr_v2_binary_operator(OP, FUNCTOR) \
249 template <class Expr1, class Expr2> \
250 inline \
251 typename \
252 std::enable_if< \
253  (details::is_vec_expr_v2_arg<Expr1>::value && \
254  details::is_vec_expr_v2_arg<Expr2>::value) || \
255  (details::is_rheolef_arithmetic<Expr1>::value && \
256  details::is_vec_expr_v2_arg<Expr2>::value) || \
257  (details::is_vec_expr_v2_arg<Expr1>::value && \
258  details::is_rheolef_arithmetic<Expr2>::value), \
259  typename \
260  details::vec_expr_v2_binary_traits< \
261  FUNCTOR, \
262  Expr1, Expr2 \
263  >::type \
264 >::type \
265 operator OP (const Expr1& expr1, const Expr2& expr2) \
266 { \
267  typedef typename details::vec_expr_v2_binary_traits <FUNCTOR, Expr1, Expr2>::type expr_t; \
268  return expr_t (FUNCTOR(), expr1, expr2); \
269 }
270 _RHEOLEF_vec_expr_v2_binary_operator(+, details::generic_plus<>)
271 _RHEOLEF_vec_expr_v2_binary_operator(-, details::generic_minus<>)
272 #undef _RHEOLEF_vec_expr_v2_binary_operator
273 
274 // c*x ; x*c
275 template <class Expr1, class Expr2>
276 inline
277 typename
278 std::enable_if<
279  (details::is_rheolef_arithmetic<Expr1>::value &&
282  details::is_rheolef_arithmetic<Expr2>::value),
283  typename
285  details::generic_multiplies<>,
286  Expr1, Expr2
287  >::type
288 >::type
289 operator* (const Expr1& expr1, const Expr2& expr2)
290 {
291  typedef details::generic_multiplies<> fun_t;
293  return expr_t (fun_t(), expr1, expr2);
294 }
295 
296 // x/c
297 template <class Expr1, class Expr2>
298 inline
299 typename
300 std::enable_if<
301  (details::is_vec_expr_v2_arg<Expr1>::value &&
302  details::is_rheolef_arithmetic<Expr2>::value),
303  typename
304  details::vec_expr_v2_binary_traits<
305  details::generic_divides<>,
306  Expr1, Expr2
307  >::type
308 >::type
309 operator/ (const Expr1& expr1, const Expr2& expr2)
310 {
311  typedef details::generic_divides<> fun_t;
313  return expr_t (fun_t(), expr1, expr2);
314 }
315 // ===================================================================
316 // 2. completing the interface of the vec<T,M> class
317 // ===================================================================
318 // ---------------------------------------------------------------------------
319 // 2.1. cstor & assignment
320 // ---------------------------------------------------------------------------
321 // vec<double> x = expr;
322 template<class T, class M>
323 template<class Expr, class Sfinae>
324 inline
325 vec<T,M>::vec (const Expr& expr)
326  : disarray<T,M>()
327 {
328  operator= (expr);
329 }
330 // x = expr;
331 template< class T, class M>
332 template <class Expr, class Sfinae>
333 inline
334 vec<T, M>&
336 {
337  if (disarray<T,M>::dis_size() == 0) {
338  resize (expr.ownership());
339  } else {
340  std::size_t n = disarray<T,M>::size();
341  check_macro (n == expr.size(),
342  "vec = vec_expression : incompatible size "
343  << n << " and " << expr.size());
344  }
346  return *this;
347 }
348 // ---------------------------------------------------------------------------
349 // 2.2) computed assignment
350 // ---------------------------------------------------------------------------
351 // x +-= expr;
352 #define _RHEOLEF_vec_expr_v2_op_assign(OP, FUNCTOR) \
353 template<class T, class M, class Expr> \
354 inline \
355 typename std::enable_if< \
356  details::is_vec_expr_v2_arg<Expr>::value, \
357  vec<T,M>&>::type \
358 operator OP (vec<T,M>& x, const Expr& expr) \
359 { \
360  check_macro (x.size() == expr.size(), "vec " << #OP << " vec_expression : incompatible spaces " \
361  << x.size() << " and " << expr.size()); \
362  details::assign_with_operator (x.begin(), x.end(), expr.begin(), FUNCTOR()); \
363  return x; \
364 }
367 #undef _RHEOLEF_vec_expr_v2_op_assign
368 
369 // x -+*/= c
370 #define _RHEOLEF_vec_expr_v2_op_assign_constant(OP, FUNCTOR) \
371 template<class T, class M, class Expr> \
372 inline \
373 typename std::enable_if< \
374  details::is_rheolef_arithmetic<Expr>::value, \
375  vec<T,M>&>::type \
376 operator OP (vec<T,M>& x, const Expr& expr) \
377 { \
378  details::assign_with_operator (x.begin(), x.end(), details::iterator_on_constant<Expr>(expr), FUNCTOR()); \
379  return x; \
380 }
384 _RHEOLEF_vec_expr_v2_op_assign_constant(/=, details::divides_assign)
385 #undef _RHEOLEF_vec_expr_v2_op_assign_constant
386 
387 // ---------------------------------------------------------------------------
388 // 2.3. scalar product
389 // ---------------------------------------------------------------------------
391 template <class Expr1, class Expr2>
392 inline
393 typename
394 std::enable_if<
397  typename promote<
398  typename Expr1::float_type,
399  typename Expr2::float_type>::type
400 >::type
401 dot (const Expr1& expr1, const Expr2& expr2)
402 {
403  typedef typename Expr1::memory_type M;
404  return dis_inner_product (expr1.begin(), expr2.begin(), expr1.size(), expr1.ownership().comm(), M());
405 }
407 template <class Expr1, class Expr2>
408 inline
409 typename
410 std::enable_if<
411  details::is_rheolef_arithmetic<Expr1>::value &&
413  typename Expr2::float_type
414 >::type
415 dot (const Expr1& expr1, const Expr2& expr2)
416 {
417  typedef typename Expr2::memory_type M;
418  return expr1*dis_accumulate (expr2.begin(), expr2.size(), expr2.ownership().comm(), M());
419 }
421 template <class Expr1, class Expr2>
422 inline
423 typename
424 std::enable_if<
425  details::is_vec_expr_v2_arg<Expr1>::value &&
426  details::is_rheolef_arithmetic<Expr2>::value,
427  typename Expr1::float_type
428 >::type
429 dot (const Expr1& expr1, const Expr2& expr2)
430 {
431  typedef typename Expr1::memory_type M;
432  return dis_accumulate (expr1.begin(), expr1.size(), expr1.ownership().comm(), M())*expr2;
433 }
434 
435 } // namespace rheolef
436 #endif // _RHEOLEF_VEC_EXPR_v2_H
field::size_type size_type
Definition: branch.cc:425
see the disarray page for the full documentation
Definition: disarray.h:459
see the distributor page for the full documentation
Definition: distributor.h:62
size_type size(size_type iproc) const
Definition: distributor.h:163
see the vec page for the full documentation
Definition: vec.h:79
vec< T, M > & operator=(const vec< T, M > &x)
Definition: vec.h:175
vec(const vec< T, M > &)
Definition: vec.h:168
rheolef::std type
rheolef::std Expr1
dot(x,y): see the expression page for the full documentation
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void assign_with_operator(ForwardIterator first, ForwardIterator last, InputIterator iter_rhs, OpAssign op_assign)
This file is part of Rheolef.
_RHEOLEF_vec_expr_v2_unary_operator(+, details::generic_unary_plus<>) _RHEOLEF_vec_expr_v2_unary_operator(-
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
std::iterator_traits< InputIterator >::value_type dis_accumulate(InputIterator first, Size n, const distributor::communicator_type &comm, sequential)
_RHEOLEF_vec_expr_v2_binary_operator(+, details::generic_plus<>) _RHEOLEF_vec_expr_v2_binary_operator(-
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
const_iterator(Op op, expr1_const_iterator iter1, expr2_const_iterator iter2)
Definition: vec_expr_v2.h:194
promote< typename Expr1::value_type, typename Expr2::value_type >::type value_type
Definition: vec_expr_v2.h:190
const_iterator begin() const
Definition: vec_expr_v2.h:203
Expr2::const_iterator expr2_const_iterator
Definition: vec_expr_v2.h:164
Expr1::const_iterator expr1_const_iterator
Definition: vec_expr_v2.h:163
vec_expr_v2_binary(const Op &op, const Expr1 &expr1, const Expr2 &expr2)
Definition: vec_expr_v2.h:169
const distributor & ownership() const
Definition: vec_expr_v2.h:181
promote< typename Expr1::value_type, typename Expr2::value_type >::type value_type
Definition: vec_expr_v2.h:161
float_traits< value_type >::type float_type
Definition: vec_expr_v2.h:165
const_iterator(Op op, expr_const_iterator expr_iter)
Definition: vec_expr_v2.h:111
const_iterator begin() const
Definition: vec_expr_v2.h:119
Expr::const_iterator expr_const_iterator
Definition: vec_expr_v2.h:82
const expr_const_iterator _expr_iter
Definition: vec_expr_v2.h:122
vec_expr_v2_unary(const Op &op, const Expr &expr)
Definition: vec_expr_v2.h:87
const distributor & ownership() const
Definition: vec_expr_v2.h:100
vec_expr_v2_unary(const BinaryOp &binop, const Constant &c, const Expr &expr)
Definition: vec_expr_v2.h:91
vec_expr_v2_unary(const BinaryOp &binop, const Expr &expr, const Constant &c)
Definition: vec_expr_v2.h:95
float_traits< value_type >::type float_type
Definition: vec_expr_v2.h:83
Expr1::memory_type M
Definition: vec_expr_v2.h:416
return dis_inner_product(expr1.begin(), expr2.begin(), expr1.size(), expr1.ownership().comm(), M())
#define _RHEOLEF_vec_expr_v2_op_assign(OP, FUNCTOR)
Definition: vec_expr_v2.h:352
#define _RHEOLEF_vec_expr_v2_op_assign_constant(OP, FUNCTOR)
Definition: vec_expr_v2.h:370