Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
50namespace rheolef {
51
52// ===================================================================
53// 1. linear/affine algebra operators
54// ===================================================================
55// -------------------------------------------------------------------
56// 1.1. vec_expr, a type concept for expression types
57// -------------------------------------------------------------------
58namespace details {
59
60// Define a trait type for detecting vec expression valid arguments
61template<class T> struct is_vec : std::false_type {};
62template<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
65template<class T> struct is_vec_expr_v2_arg : std::false_type {};
66template<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// -------------------------------------------
72namespace details {
73
74template <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
106 typedef std::forward_iterator_tag iterator_category;
107 typedef typename Expr::value_type value_type;
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 };
120protected:
121 const Op _op;
124};
125template<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) \
130template <class Expr> \
131inline \
132typename \
133std::enable_if< \
134 details::is_vec_expr_v2_arg<Expr>::value, \
135 details::vec_expr_v2_unary< \
136 FUNCTOR, \
137 Expr \
138 > \
139>::type \
140operator 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// -------------------------------------------
153namespace details {
154
155template <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
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;
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 };
204protected:
205 const Op _op;
209};
210template<class Op, class Expr1, class Expr2> struct is_vec_expr_v2_arg <vec_expr_v2_binary<Op,Expr1,Expr2> > : std::true_type {};
211
212template <class Op, class Expr1, class Expr2, class Sfinae = void>
213struct vec_expr_v2_binary_traits { /* catch-all case */ };
214
215// vec_expr +- vec_expr
216template <class Op, class Expr1, class Expr2>
217struct 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{
222 typedef vec_expr_v2_binary <Op,Expr1,Expr2> type;
223};
224// constant +-* vec_expr
225template <class Op, class Expr1, class Expr2>
226struct 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{
231 typedef generic_binder1st <Op, Expr1> fun_t;
232 typedef vec_expr_v2_unary <fun_t,Expr2> type;
233};
234// vec_expr +-* constant
235template <class Op, class Expr1, class Expr2>
236struct 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{
241 typedef generic_binder2nd <Op, Expr2> fun_t;
242 typedef vec_expr_v2_unary <fun_t,Expr1> type;
243};
244
245} // namespace details
246
247// x+y; x+c ; c+x
248#define _RHEOLEF_vec_expr_v2_binary_operator(OP, FUNCTOR) \
249template <class Expr1, class Expr2> \
250inline \
251typename \
252std::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 \
265operator 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
275template <class Expr1, class Expr2>
276inline
277typename
278std::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
289operator* (const Expr1& expr1, const Expr2& expr2)
290{
291 typedef details::generic_multiplies<> fun_t;
292 typedef typename details::vec_expr_v2_binary_traits <fun_t, Expr1, Expr2>::type expr_t;
293 return expr_t (fun_t(), expr1, expr2);
294}
295
296// x/c
297template <class Expr1, class Expr2>
298inline
299typename
300std::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
309operator/ (const Expr1& expr1, const Expr2& expr2)
310{
311 typedef details::generic_divides<> fun_t;
312 typedef typename details::vec_expr_v2_binary_traits <fun_t, Expr1, Expr2>::type expr_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;
322template<class T, class M>
323template<class Expr, class Sfinae>
324inline
325vec<T,M>::vec (const Expr& expr)
326 : disarray<T,M>()
327{
328 operator= (expr);
329}
330// x = expr;
331template< class T, class M>
332template <class Expr, class Sfinae>
333inline
335vec<T,M>::operator= (const Expr& expr)
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) \
353template<class T, class M, class Expr> \
354inline \
355typename std::enable_if< \
356 details::is_vec_expr_v2_arg<Expr>::value, \
357 vec<T,M>&>::type \
358operator 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) \
371template<class T, class M, class Expr> \
372inline \
373typename std::enable_if< \
374 details::is_rheolef_arithmetic<Expr>::value, \
375 vec<T,M>&>::type \
376operator 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// ---------------------------------------------------------------------------
391template <class Expr1, class Expr2>
392inline
393typename
394std::enable_if<
397 typename promote<
398 typename Expr1::float_type,
399 typename Expr2::float_type>::type
400>::type
401dot (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}
407template <class Expr1, class Expr2>
408inline
409typename
410std::enable_if<
411 details::is_rheolef_arithmetic<Expr1>::value &&
413 typename Expr2::float_type
414>::type
415dot (const Expr1& expr1, const Expr2& expr2)
417 typedef typename Expr2::memory_type M;
418 return expr1*dis_accumulate (expr2.begin(), expr2.size(), expr2.ownership().comm(), M());
419}
421template <class Expr1, class Expr2>
422inline
423typename
424std::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
429dot (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
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type size(size_type iproc) const
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
Expr1::float_type T
Definition field_expr.h:230
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.
std::iterator_traits< InputIterator >::value_type dis_accumulate(InputIterator first, Size n, const distributor::communicator_type &comm, sequential)
promote< typenamestd::iterator_traits< InputIterator1 >::value_type, typenamestd::iterator_traits< InputIterator2 >::value_type >::type dis_inner_product(InputIterator1 first1, InputIterator2 first2, Size n, const distributor::communicator_type &comm, sequential)
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
STL namespace.
const_iterator(Op op, expr1_const_iterator iter1, expr2_const_iterator iter2)
promote< typenameExpr1::value_type, typenameExpr2::value_type >::type value_type
Expr2::const_iterator expr2_const_iterator
Expr1::const_iterator expr1_const_iterator
vec_expr_v2_binary(const Op &op, const Expr1 &expr1, const Expr2 &expr2)
promote< typenameExpr1::value_type, typenameExpr2::value_type >::type value_type
const distributor & ownership() const
float_traits< value_type >::type float_type
const_iterator(Op op, expr_const_iterator expr_iter)
Expr::const_iterator expr_const_iterator
Definition vec_expr_v2.h:82
const expr_const_iterator _expr_iter
vec_expr_v2_unary(const Op &op, const Expr &expr)
Definition vec_expr_v2.h:87
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
const distributor & ownership() const
float_traits< value_type >::type float_type
Definition vec_expr_v2.h:83
Expr1::memory_type M
#define _RHEOLEF_vec_expr_v2_op_assign(OP, FUNCTOR)
#define _RHEOLEF_vec_expr_v2_op_assign_constant(OP, FUNCTOR)
#define _RHEOLEF_vec_expr_v2_unary_operator(OP, FUNCTOR)
#define _RHEOLEF_vec_expr_v2_binary_operator(OP, FUNCTOR)