dune-istl  2.2.0
multitypeblockvector.hh
Go to the documentation of this file.
1 #ifndef DUNE_MULTITYPEVECTOR_HH
2 #define DUNE_MULTITYPEVECTOR_HH
3 
4 #if HAVE_BOOST
5 #ifdef HAVE_BOOST_FUSION
6 
7 #include<cmath>
8 #include<iostream>
9 
10 #include "istlexception.hh"
11 
12 #include <boost/fusion/sequence.hpp>
13 #include <boost/fusion/container.hpp>
14 #include <boost/fusion/iterator.hpp>
15 #include <boost/typeof/typeof.hpp>
16 #include <boost/fusion/algorithm.hpp>
17 
18 namespace mpl=boost::mpl;
19 namespace fusion=boost::fusion;
20 
21 // forward decl
22 namespace Dune {
23  template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
24  typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
25  typename T8=fusion::void_, typename T9=fusion::void_>
26  class MultiTypeBlockVector;
27 }
28 
29 #include "gsetc.hh"
30 
31 namespace Dune {
32 
56  template<int current_element, int remaining_elements, typename TVec>
57  class MultiTypeBlockVector_Print {
58  public:
62  static void print(const TVec& v) {
63  std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
64  MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v); //next element
65  }
66  };
67  template<int current_element, typename TVec> //recursion end (remaining_elements=0)
68  class MultiTypeBlockVector_Print<current_element,0,TVec> {
69  public:
70  static void print(const TVec& v) {std::cout << "\n";}
71  };
72 
73 
74 
86  template<int count, typename T1, typename T2>
87  class MultiTypeBlockVector_Ident {
88  public:
89 
94  static void equalize(T1& a, const T2& b) {
95  fusion::at_c<count-1>(a) = b; //equalize current elements
96  MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b); //next elements
97  }
98  };
99  template<typename T1, typename T2> //recursion end (count=0)
100  class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
101 
102 
103 
104 
105 
106 
112  template<int count, typename T>
113  class MultiTypeBlockVector_Add {
114  public:
115 
119  static void add (T& a, const T& b) { //add vector elements
120  fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
121  MultiTypeBlockVector_Add<count-1,T>::add(a,b);
122  }
123 
127  static void sub (T& a, const T& b) { //sub vector elements
128  fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
129  MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
130  }
131  };
132  template<typename T> //recursion end; specialization for count=0
133  class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
134 
135 
136 
142  template<int count, typename TVec, typename Ta>
143  class MultiTypeBlockVector_AXPY {
144  public:
145 
149  static void axpy(TVec& x, const Ta& a, const TVec& y) {
150  fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
151  MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
152  }
153  };
154  template<typename TVec, typename Ta> //specialization for count=0
155  class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
156 
157 
163  template<int count, typename TVec, typename Ta>
164  class MultiTypeBlockVector_Mulscal {
165  public:
166 
170  static void mul(TVec& x, const Ta& a) {
171  fusion::at_c<(count-1)>(x) *= a;
172  MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
173  }
174  };
175  template<typename TVec, typename Ta> //specialization for count=0
176  class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
177 
178 
179 
186  template<int count, typename TVec>
187  class MultiTypeBlockVector_Mul {
188  public:
189  static typename TVec::field_type mul(const TVec& x, const TVec& y) {
190  return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y);
191  }
192  };
193  template<typename TVec>
194  class MultiTypeBlockVector_Mul<0,TVec> {
195  public: static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
196  };
197 
198 
199 
200 
201 
208  template<int count, typename T>
209  class MultiTypeBlockVector_Norm {
210  public:
211 
215  static double result (const T& a) { //result = sum of all elements' 2-norms
216  return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
217  }
218  };
219 
220  template<typename T> //recursion end: no more vector elements to add...
221  class MultiTypeBlockVector_Norm<0,T> {
222  public:
223  static double result (const T& a) {return 0.0;}
224  };
225 
234  template<typename T1, typename T2, typename T3, typename T4,
235  typename T5, typename T6, typename T7, typename T8, typename T9>
236  class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
237 
238  public:
239 
243  typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
244 
245  typedef typename T1::field_type field_type;
246 
250  const int count() {return mpl::size<type>::value;}
251 
255  template<typename T>
256  void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
257 
261  void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
262 
266  void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
267 
268  void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
269  void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
270  void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
271 
272  field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
273 
277  double two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
278 
282  double two_norm() const {return sqrt(this->two_norm2());}
283 
287  template<typename Ta>
288  void axpy (const Ta& a, const type& y) {
289  MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
290  }
291 
292  };
293 
294 
295 
301  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
302  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
303  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);
304  return s;
305  }
306 
307 
308 
309 } // end namespace
310 
311 #endif // end HAVE_BOOST_FUSION
312 #endif // end HAVE_BOOST
313 
314 #endif