rheolef  6.3
vec_expr.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_VEC_EXPR_H
2 #define _RHEOLEF_VEC_EXPR_H
3 #include "rheolef/vec.h"
4 #include "rheolef/dis_inner_product.h"
5 #include "rheolef/dis_accumulate.h"
6 #include "rheolef/pretty_name.h"
7 
8 #include <boost/mpl/bool.hpp>
9 #include <boost/proto/core.hpp>
10 #include <boost/proto/debug.hpp>
11 #include <boost/proto/context.hpp>
12 #include <boost/proto/transform.hpp>
13 #include <boost/proto/operators.hpp>
14 #include <boost/utility/enable_if.hpp>
15 #include <boost/typeof/std/vector.hpp>
16 #include <boost/typeof/std/complex.hpp>
17 #include <boost/type_traits/remove_reference.hpp>
18 
19 namespace rheolef {
20 namespace mpl = boost::mpl;
21 namespace proto = boost::proto;
22 using proto::_;
23 
24 template <class Expr>
25 struct vec_expr;
26 
27 template<class Iterator>
28 struct vec_iterator_wrapper {
29  typedef Iterator iterator;
30  explicit vec_iterator_wrapper (Iterator iter1) : iter(iter1) {}
31  mutable Iterator iter;
32 };
33 struct vec_begin : proto::callable {
34  template<class Sig>
35  struct result;
36 
37  template<class This, class Container>
38  struct result<This(Container)> : proto::result_of::as_expr
39  <
40  vec_iterator_wrapper<
41  typename boost::remove_reference<Container>::type::const_iterator>
42  > {};
43 
44  template<class Container>
45  typename result<vec_begin(const Container&)>::type
46  operator ()(const Container& cont) const {
47  vec_iterator_wrapper<typename Container::const_iterator> iter (cont.begin());
48  return proto::as_expr (iter);
49  }
50 };
51 struct vec_grammar_begin
52  : proto::or_<
53  proto::when <proto::terminal<vec<_, _> >, vec_begin(proto::_value)>
54  , proto::when <proto::terminal<_> >
55  , proto::when <proto::nary_expr<_, proto::vararg<vec_grammar_begin> > >
56  >
57 {};
58 
59 struct vec_dereference_context {
60  template<class Expr, class EnableIf = void>
61  struct eval : proto::default_eval<Expr, const vec_dereference_context> {};
62 
63  template<class Expr>
64  struct eval<Expr,
65  typename boost::enable_if<
66  proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > >
67  >::type
68  > {
69  typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
70  typedef typename IteratorWrapper::iterator iterator;
71  typedef typename std::iterator_traits<iterator>::reference result_type;
72 
73  result_type operator() (Expr &expr, const vec_dereference_context&) const {
74  return *proto::value(expr).iter;
75  }
76  };
77 };
78 struct vec_increment_context {
79  template<class Expr, class EnableIf = void>
80  struct eval : proto::null_eval<Expr, const vec_increment_context> {};
81 
82  template<class Expr> struct eval<Expr,
83  typename boost::enable_if<
84  proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > >
85  >::type
86  > {
87  typedef void result_type;
88 
89  result_type operator() (Expr &expr, const vec_increment_context&) const {
90  ++proto::value(expr).iter;
91  }
92  };
93 };
94 struct vec_subscript_context {
95  typedef vec<Float>::size_type size_type;
96  vec_subscript_context (size_type i) : _i(i) {}
97  template<class Expr, class EnableIf = void>
98  struct eval : proto::default_eval<Expr, const vec_subscript_context> {};
99  template<class Expr>
100  struct eval<Expr, typename boost::enable_if<
101  proto::matches<Expr, proto::terminal<vec<_, _> > > >::type> {
102  typedef typename proto::result_of::value<Expr>::type::value_type result_type;
103  result_type operator() (Expr &expr, const vec_subscript_context& ctx) const {
104  return proto::value(expr)[ctx._i];
105  }
106  };
107  template<class Expr>
108  struct eval<Expr, typename boost::enable_if<
109  proto::matches<Expr, proto::terminal<vec_iterator_wrapper<_> > > >::type > {
110  typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
111  typedef typename IteratorWrapper::iterator iterator;
112  typedef typename std::iterator_traits<iterator>::reference result_type;
113  result_type operator() (Expr &expr, const vec_subscript_context& ctx) const {
114  return *proto::value(expr).iter[ctx._i];
115  }
116  };
117  mutable size_type _i;
118 };
119 struct vec_get_size_context {
120  typedef vec<Float>::size_type size_type;
121  vec_get_size_context () {}
122 
123  template<class Expr, class EnableIf = void>
124  struct eval : proto::null_eval<Expr, const vec_get_size_context> {};
125 
126  template<class Expr>
127  struct eval<Expr, typename boost::enable_if<
128  proto::matches<Expr, proto::terminal<vec<_,_> > > >::type> {
129  typedef void result_type;
130  result_type operator() (Expr &expr, const vec_get_size_context& ctx) const {
131  ctx._ownership = proto::value(expr).ownership();
132  }
133  };
134  size_type size() const { return _ownership.size(); };
135  distributor ownership() const { return _ownership; };
136  mutable distributor _ownership;
137 };
138 struct vec_check_size_context {
139  typedef vec<Float>::size_type size_type;
140  vec_check_size_context (size_type size) : _size(size) {}
141 
142  template<class Expr, class EnableIf = void>
143  struct eval : proto::null_eval<Expr, const vec_check_size_context> {};
144 
145  template<class Expr>
146  struct eval<Expr, typename boost::enable_if<
147  proto::matches<Expr, proto::terminal<vec<_,_> > > >::type> {
148  typedef void result_type;
149  result_type operator() (Expr &expr, const vec_check_size_context& ctx) const {
150  if (ctx._size != proto::value(expr).size()) {
151  error_macro ("incompatible sizes "<<ctx._size<< " and "
152  << proto::value(expr).size() << " in vec<T> expression "
153  <<typename_macro(Expr));
154  }
155  }
156  };
157  size_type _size;
158 };
159 struct vec_assign_operators : proto::switch_<struct vec_assign_operators_cases> {};
160 
161 struct vec_assign_operators_cases {
162  template<class Tag, int D = 0> struct case_ : proto::not_<_> {};
163 
164  template<int D> struct case_< proto::tag::plus_assign, D > : _ {};
165  template<int D> struct case_< proto::tag::minus_assign, D > : _ {};
166  template<int D> struct case_< proto::tag::multiplies_assign, D > : _ {};
167  template<int D> struct case_< proto::tag::divides_assign, D > : _ {};
168  template<int D> struct case_< proto::tag::modulus_assign, D > : _ {};
169  template<int D> struct case_< proto::tag::shift_left_assign, D > : _ {};
170  template<int D> struct case_< proto::tag::shift_right_assign, D > : _ {};
171  template<int D> struct case_< proto::tag::bitwise_and_assign, D > : _ {};
172  template<int D> struct case_< proto::tag::bitwise_or_assign, D > : _ {};
173  template<int D> struct case_< proto::tag::bitwise_xor_assign, D > : _ {};
174 };
175 struct vec_grammar
176  : proto::or_<
177  proto::terminal<_>
178  , proto::and_<
179  proto::nary_expr<_, proto::vararg<vec_grammar> >
180  , proto::not_<vec_assign_operators>
181  >
182  >
183 {};
184 
185 // Expressions in the vec_domain will be wrapped in vec_expr<>
186 struct vec_domain : proto::domain<proto::generator<vec_expr>, vec_grammar> {};
187 
188 // It mimics the array<T> and vec<T> interface
189 template<class Expr>
190 struct vec_expr : proto::extends<Expr, vec_expr<Expr>, vec_domain> {
191 
192 
193  typedef proto::extends<Expr, vec_expr<Expr>, vec_domain> base_type;
194 
195  typedef vec_get_size_context::size_type size_type;
196  typedef typename boost::result_of<vec_grammar_begin(const Expr&)>::type raw_iterator;
197 
198 
199  explicit vec_expr (const Expr& expr)
200  : base_type(expr)
201  {}
202 
203 
204  distributor ownership() const {
205  const vec_get_size_context get_size;
206  proto::eval (*this, get_size);
207  return get_size.ownership();
208  }
209  size_type size() const {
210  return ownership().size();
211  }
212  size_type dis_size() const {
213  return ownership().dis_size();
214  }
215  struct const_iterator_begin : raw_iterator {
216 
217  typedef std::input_iterator_tag iterator_category;
218  typedef typename proto::result_of::eval<raw_iterator,vec_dereference_context>::type value_type;
219  typedef value_type reference;
220  typedef value_type* pointer;
221  typedef std::ptrdiff_t difference_type;
222 
223 #if BOOST_VERSION < 104601
224  const_iterator_begin (const vec_expr<Expr>& expr)
225  : raw_iterator(vec_grammar_begin() (expr))
226  {}
227 #else // BOOST_VERSION < 104601
228  const_iterator_begin (const vec_expr<Expr>& expr)
229  : raw_iterator(vec_grammar_begin() (expr.proto_expr_))
230  // ICI: il faut recuperer base_type(expr) de type Expr au lieu vec_expr<Expr>
231  {}
232 #endif // BOOST_VERSION < 104601
233  const_iterator_begin& operator++ () {
234  const vec_increment_context inc = {};
235  proto::eval (*this, inc);
236  return *this;
237  }
238  const_iterator_begin operator++ (int) {
239  const_iterator_begin tmp = *this;
240  operator++();
241  return tmp;
242  }
243  reference dereference () const {
244  const vec_dereference_context deref = {};
245  typedef typename proto::result_of::eval<raw_iterator,vec_dereference_context>::type my_ref;
246  reference value = proto::eval(*this, deref);
247  return value;
248  }
249  reference operator* () const { return dereference(); }
250  reference operator* () { return dereference(); }
251  };
252  typedef const_iterator_begin const_iterator;
253  const_iterator begin() const { return const_iterator(*this); }
254 
255 
256  typedef typename proto::result_of::eval<const Expr, const vec_subscript_context>::type value_type;
257  typedef value_type scalar_type;
258  value_type operator[] (size_type i) const {
259  const vec_subscript_context get_subscript(i);
260  return proto::eval(*this, get_subscript);
261  }
262 private:
263  using proto::extends<Expr, vec_expr<Expr>, vec_domain>::operator[];
264 };
265 template<class T>
266 struct is_vec : mpl::false_ {};
267 
268 template<class T, class M>
269 struct is_vec<vec<T, M> > : mpl::true_ {};
270 
271 namespace vec_detail {
272  template<class ForwardIterator, class Expr, class Op>
273  void evaluate (ForwardIterator begin, ForwardIterator end, const Expr& expr, Op op) {
274  const vec_increment_context inc = {};
275  const vec_dereference_context deref = {};
276  typename boost::result_of<vec_grammar_begin(const Expr&)>::type expr_begin
277  = vec_grammar_begin() (expr);
278 #ifndef TRY_VEC_EXPR
279  for (; begin != end; ++begin) {
280  op (*begin, proto::eval(expr_begin, deref));
281  proto::eval (expr_begin, inc);
282  }
283 #endif // TRY_VEC_EXPR
284  }
285  struct assign_op {
286  template<class T, class U>
287  void operator() (T &t, U const &u) const { t = u; }
288  };
289 } // namespace vec_detail
290 
291 template<class T, class M>
292 template<class Expr>
293 inline
294 vec<T,M>&
295 vec<T,M>::operator= (const vec_expr<Expr>& expr) {
296  const vec_get_size_context get_size;
297  proto::eval (proto::as_expr<vec_domain>(expr), get_size);
298  size_type expr_size = get_size.size();
299  if (array<T,M>::dis_size() == 0) {
300  distributor expr_ownership = get_size.ownership();
301  resize (expr_ownership);
302  }
303  const vec_check_size_context check_size (array<T,M>::size());
304  proto::eval (proto::as_expr<vec_domain>(expr), check_size); // error if the sizes don't match
305  vec_detail::evaluate (array<T,M>::begin(), array<T,M>::end(), proto::as_expr<vec_domain>(expr), vec_detail::assign_op());
306  return *this;
307 }
308 template<class T, class M>
309 template<class Expr>
310 inline
311 vec<T,M>::vec (const vec_expr<Expr>& expr)
312  : array<T,M>()
313 {
314  operator= (expr);
315 }
316 template<class T, class M>
317 inline
318 T
319 norm2 (const vec<T,M>& x)
320 {
321  return dot(x,x);
322 }
323 template<class T, class M>
324 inline
325 T
326 norm (const vec<T,M>& x)
327 {
328  return sqrt(norm2(x));
329 }
330 
331 } // namespace rheolef
332 #endif // _RHEOLEF_VEC_EXPR_H
333