dune-istl  2.2.0
multitypeblockvector.hh
Go to the documentation of this file.
00001 #ifndef DUNE_MULTITYPEVECTOR_HH
00002 #define DUNE_MULTITYPEVECTOR_HH
00003 
00004 #if HAVE_BOOST
00005 #ifdef HAVE_BOOST_FUSION
00006 
00007 #include<cmath>
00008 #include<iostream>
00009 
00010 #include "istlexception.hh"
00011 
00012 #include <boost/fusion/sequence.hpp>
00013 #include <boost/fusion/container.hpp>
00014 #include <boost/fusion/iterator.hpp>
00015 #include <boost/typeof/typeof.hpp>
00016 #include <boost/fusion/algorithm.hpp>
00017 
00018 namespace mpl=boost::mpl;
00019 namespace fusion=boost::fusion;
00020 
00021 // forward decl
00022 namespace Dune {
00023   template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
00024            typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
00025            typename T8=fusion::void_, typename T9=fusion::void_>
00026   class MultiTypeBlockVector;
00027 }
00028 
00029 #include "gsetc.hh"
00030 
00031 namespace Dune {
00032    
00056   template<int current_element, int remaining_elements, typename TVec>
00057   class MultiTypeBlockVector_Print {
00058   public:
00062     static void print(const TVec& v) {
00063       std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
00064       MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v);   //next element
00065     }
00066   };
00067   template<int current_element, typename TVec>            //recursion end (remaining_elements=0)
00068   class MultiTypeBlockVector_Print<current_element,0,TVec> {
00069   public:
00070     static void print(const TVec& v) {std::cout << "\n";}
00071   };
00072 
00073 
00074 
00086   template<int count, typename T1, typename T2>
00087   class MultiTypeBlockVector_Ident {
00088   public:
00089 
00094     static void equalize(T1& a, const T2& b) {
00095       fusion::at_c<count-1>(a) = b;           //equalize current elements
00096       MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b);    //next elements
00097     }
00098   };
00099   template<typename T1, typename T2>                      //recursion end (count=0)
00100   class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
00101 
00102 
00103 
00104 
00105 
00106 
00112   template<int count, typename T>
00113   class MultiTypeBlockVector_Add {
00114   public:
00115         
00119     static void add (T& a, const T& b) {    //add vector elements
00120       fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
00121       MultiTypeBlockVector_Add<count-1,T>::add(a,b);
00122     }
00123 
00127     static void sub (T& a, const T& b) {    //sub vector elements
00128       fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
00129       MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
00130     }
00131   };
00132   template<typename T>                                    //recursion end; specialization for count=0
00133   class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
00134 
00135 
00136 
00142   template<int count, typename TVec, typename Ta>
00143   class MultiTypeBlockVector_AXPY {
00144   public:
00145 
00149     static void axpy(TVec& x, const Ta& a, const TVec& y) {
00150       fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
00151       MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
00152     }
00153   };
00154   template<typename TVec, typename Ta>                    //specialization for count=0
00155   class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
00156 
00157 
00163   template<int count, typename TVec, typename Ta>
00164   class MultiTypeBlockVector_Mulscal {
00165   public:
00166 
00170     static void mul(TVec& x, const Ta& a) {
00171       fusion::at_c<(count-1)>(x) *= a;
00172       MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
00173     }
00174   };
00175   template<typename TVec, typename Ta>                    //specialization for count=0
00176   class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
00177 
00178 
00179 
00186   template<int count, typename TVec>
00187   class MultiTypeBlockVector_Mul {
00188   public:
00189     static typename TVec::field_type mul(const TVec& x, const TVec& y) {
00190       return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y);
00191     }
00192   };
00193   template<typename TVec>
00194   class MultiTypeBlockVector_Mul<0,TVec> {
00195   public: static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
00196   };
00197 
00198 
00199 
00200 
00201 
00208   template<int count, typename T>
00209   class MultiTypeBlockVector_Norm {
00210   public:
00211 
00215     static double result (const T& a) {             //result = sum of all elements' 2-norms
00216       return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
00217     }
00218   };
00219 
00220   template<typename T>                                    //recursion end: no more vector elements to add...
00221   class MultiTypeBlockVector_Norm<0,T> {
00222   public: 
00223     static double result (const T& a) {return 0.0;} 
00224   };
00225 
00234   template<typename T1, typename T2, typename T3, typename T4,
00235            typename T5, typename T6, typename T7, typename T8, typename T9>
00236   class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
00237 
00238   public:
00239 
00243     typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
00244 
00245     typedef typename T1::field_type field_type;
00246 
00250     const int count() {return mpl::size<type>::value;}
00251 
00255     template<typename T>
00256     void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
00257 
00261     void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
00262 
00266     void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
00267 
00268     void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
00269     void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
00270     void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
00271 
00272     field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
00273 
00277     double two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
00278 
00282     double two_norm() const {return sqrt(this->two_norm2());}
00283 
00287     template<typename Ta>
00288     void axpy (const Ta& a, const type& y) {
00289       MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
00290     }
00291 
00292   };
00293 
00294 
00295 
00301   template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
00302   std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
00303     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);
00304     return s;
00305   }
00306 
00307 
00308 
00309 } // end namespace
00310 
00311 #endif // end HAVE_BOOST_FUSION
00312 #endif // end HAVE_BOOST
00313 
00314 #endif