dune-istl
2.2.0
|
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