dune-istl  2.4.1-rc2
multitypeblockvector.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
4 #define DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
5 
6 #if HAVE_DUNE_BOOST
7 #ifdef HAVE_BOOST_FUSION
8 
9 #include <cmath>
10 #include <iostream>
11 
12 #include <dune/common/dotproduct.hh>
13 #include <dune/common/ftraits.hh>
14 
15 #include "istlexception.hh"
16 
17 #include <boost/fusion/sequence.hpp>
18 #include <boost/fusion/container.hpp>
19 #include <boost/fusion/iterator.hpp>
20 #include <boost/typeof/typeof.hpp>
21 #include <boost/fusion/algorithm.hpp>
22 
23 namespace mpl=boost::mpl;
24 namespace fusion=boost::fusion;
25 
26 // forward decl
27 namespace Dune {
28  template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
29  typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
30  typename T8=fusion::void_, typename T9=fusion::void_>
31  class MultiTypeBlockVector;
32 }
33 
34 #include "gsetc.hh"
35 
36 namespace Dune {
37 
61  template<int current_element, int remaining_elements, typename TVec>
62  class MultiTypeBlockVector_Print {
63  public:
67  static void print(const TVec& v) {
68  std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
69  MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v); //next element
70  }
71  };
72  template<int current_element, typename TVec> //recursion end (remaining_elements=0)
73  class MultiTypeBlockVector_Print<current_element,0,TVec> {
74  public:
75  static void print(const TVec& v) {std::cout << "\n";}
76  };
77 
78 
79 
91  template<int count, typename T1, typename T2>
92  class MultiTypeBlockVector_Ident {
93  public:
94 
99  static void equalize(T1& a, const T2& b) {
100  fusion::at_c<count-1>(a) = b; //equalize current elements
101  MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b); //next elements
102  }
103  };
104  template<typename T1, typename T2> //recursion end (count=0)
105  class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
106 
107 
108 
109 
110 
111 
117  template<int count, typename T>
118  class MultiTypeBlockVector_Add {
119  public:
120 
124  static void add (T& a, const T& b) { //add vector elements
125  fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
126  MultiTypeBlockVector_Add<count-1,T>::add(a,b);
127  }
128 
132  static void sub (T& a, const T& b) { //sub vector elements
133  fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
134  MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
135  }
136  };
137  template<typename T> //recursion end; specialization for count=0
138  class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
139 
140 
141 
147  template<int count, typename TVec, typename Ta>
148  class MultiTypeBlockVector_AXPY {
149  public:
150 
154  static void axpy(TVec& x, const Ta& a, const TVec& y) {
155  fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
156  MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
157  }
158  };
159  template<typename TVec, typename Ta> //specialization for count=0
160  class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
161 
162 
168  template<int count, typename TVec, typename Ta>
169  class MultiTypeBlockVector_Mulscal {
170  public:
171 
175  static void mul(TVec& x, const Ta& a) {
176  fusion::at_c<(count-1)>(x) *= a;
177  MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
178  }
179  };
180  template<typename TVec, typename Ta> //specialization for count=0
181  class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
182 
183 
184 
193  template<int count, typename TVec>
194  class MultiTypeBlockVector_Mul {
195  public:
196  static typename TVec::field_type mul(const TVec& x, const TVec& y) { return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y); }
197  static typename TVec::field_type dot(const TVec& x, const TVec& y) { return (Dune::dot(fusion::at_c<count-1>(x),fusion::at_c<count-1>(y))) + MultiTypeBlockVector_Mul<count-1,TVec>::dot(x,y); }
198  };
199  template<typename TVec>
200  class MultiTypeBlockVector_Mul<0,TVec> {
201  public:
202  static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
203  static typename TVec::field_type dot(const TVec& x, const TVec& y) {return 0;}
204  };
205 
206 
207 
208 
209 
216  template<int count, typename T>
217  class MultiTypeBlockVector_Norm {
218  public:
219  typedef typename T::field_type field_type;
220  typedef typename FieldTraits<field_type>::real_type real_type;
221 
225  static real_type result (const T& a) { //result = sum of all elements' 2-norms
226  return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
227  }
228  };
229 
230  template<typename T> //recursion end: no more vector elements to add...
231  class MultiTypeBlockVector_Norm<0,T> {
232  public:
233  typedef typename T::field_type field_type;
234  typedef typename FieldTraits<field_type>::real_type real_type;
235  static real_type result (const T& a) {return 0.0;}
236  };
237 
246  template<typename T1, typename T2, typename T3, typename T4,
247  typename T5, typename T6, typename T7, typename T8, typename T9>
248  class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
249 
250  public:
251 
255  typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
256 
257  typedef typename T1::field_type field_type;
258 
262  int count() {return mpl::size<type>::value;}
263 
282  template< std::size_t index >
283  auto
284  operator[] ( const std::integral_constant< std::size_t, index > indexVariable )
285  -> decltype(fusion::at_c<index>(*this))
286  {
287  DUNE_UNUSED_PARAMETER(indexVariable);
288  return fusion::at_c<index>(*this);
289  }
290 
296  template< std::size_t index >
297  auto
298  operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const
299  -> decltype(fusion::at_c<index>(*this))
300  {
301  DUNE_UNUSED_PARAMETER(indexVariable);
302  return fusion::at_c<index>(*this);
303  }
304 
308  template<typename T>
309  void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
310 
314  void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
315 
319  void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
320 
321  void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
322  void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
323  void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
324 
325  field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
326  field_type dot (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::dot(*this,newv);}
327 
331  typename FieldTraits<field_type>::real_type two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
332 
336  typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
337 
341  template<typename Ta>
342  void axpy (const Ta& a, const type& y) {
343  MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
344  }
345 
346  };
347 
348 
349 
355  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
356  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
357  MultiTypeBlockVector_Print<0,mpl::size<MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value,MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(v);
358  return s;
359  }
360 
361 
362 
363 } // end namespace
364 
365 #endif // end HAVE_BOOST_FUSION
366 #endif // end HAVE_DUNE_BOOST
367 
368 #endif
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Definition: basearray.hh:19