dune-istl  2.2.0
multitypeblockmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_MultiTypeMATRIX_HH
2 #define DUNE_MultiTypeMATRIX_HH
3 
4 #include<cmath>
5 #include<iostream>
6 
7 #include "istlexception.hh"
8 
9 #if HAVE_BOOST
10 #ifdef HAVE_BOOST_FUSION
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 {
24  template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
25  typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
26  typename T8=fusion::void_, typename T9=fusion::void_>
27  class MultiTypeBlockMatrix;
28 
29  template<int I, int crow, int remain_row>
30  class MultiTypeBlockMatrix_Solver;
31 }
32 
33 #include "gsetc.hh"
34 
35 namespace Dune {
36 
54  template<int crow, int remain_rows, int ccol, int remain_cols,
55  typename TMatrix>
56  class MultiTypeBlockMatrix_Print {
57  public:
58 
62  static void print(const TMatrix& m) {
63  std::cout << "\t(" << crow << ", " << ccol << "): \n" << fusion::at_c<ccol>( fusion::at_c<crow>(m));
64  MultiTypeBlockMatrix_Print<crow,remain_rows,ccol+1,remain_cols-1,TMatrix>::print(m); //next column
65  }
66  };
67  template<int crow, int remain_rows, int ccol, typename TMatrix> //specialization for remain_cols=0
68  class MultiTypeBlockMatrix_Print<crow,remain_rows,ccol,0,TMatrix> {
69  public: static void print(const TMatrix& m) {
70  static const int xlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
71  MultiTypeBlockMatrix_Print<crow+1,remain_rows-1,0,xlen,TMatrix>::print(m); //next row
72  }
73  };
74 
75  template<int crow, int ccol, int remain_cols, typename TMatrix> //recursion end: specialization for remain_rows=0
76  class MultiTypeBlockMatrix_Print<crow,0,ccol,remain_cols,TMatrix> {
77  public:
78  static void print(const TMatrix& m)
79  {std::cout << std::endl;}
80  };
81 
82 
83 
84  //make MultiTypeBlockVector_Ident known (for MultiTypeBlockMatrix_Ident)
85  template<int count, typename T1, typename T2>
86  class MultiTypeBlockVector_Ident;
87 
88 
101  template<int rowcount, typename T1, typename T2>
102  class MultiTypeBlockMatrix_Ident {
103  public:
104 
109  static void equalize(T1& a, const T2& b) {
110  MultiTypeBlockVector_Ident< mpl::size< typename mpl::at_c<T1,rowcount-1>::type >::value ,T1,T2>::equalize(a,b); //rows are cvectors
111  MultiTypeBlockMatrix_Ident<rowcount-1,T1,T2>::equalize(a,b); //iterate over rows
112  }
113  };
114 
115  //recursion end for rowcount=0
116  template<typename T1, typename T2>
117  class MultiTypeBlockMatrix_Ident<0,T1,T2> {
118  public:
119  static void equalize (T1& a, const T2& b)
120  {}
121  };
122 
128  template<int crow, int remain_rows, int ccol, int remain_cols,
129  typename TVecY, typename TMatrix, typename TVecX>
130  class MultiTypeBlockMatrix_VectMul {
131  public:
132 
136  static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
137  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).umv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
138  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::umv(y, A, x);
139  }
140 
144  static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
145  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
146  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::mmv(y, A, x);
147  }
148 
149  template<typename AlphaType>
150  static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
151  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).usmv(alpha, fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
152  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
153  }
154 
155 
156  };
157 
158  //specialization for remain_cols = 0
159  template<int crow, int remain_rows,int ccol, typename TVecY,
160  typename TMatrix, typename TVecX>
161  class MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol,0,TVecY,TMatrix,TVecX> { //start iteration over next row
162 
163  public:
167  static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
168  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
169  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::umv(y, A, x);
170  }
171 
175  static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
176  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
177  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::mmv(y, A, x);
178  }
179 
180  template <typename AlphaType>
181  static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
182  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
183  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
184  }
185  };
186 
187  //specialization for remain_rows = 0
188  template<int crow, int ccol, int remain_cols, typename TVecY,
189  typename TMatrix, typename TVecX>
190  class MultiTypeBlockMatrix_VectMul<crow,0,ccol,remain_cols,TVecY,TMatrix,TVecX> {
191  //end recursion
192  public:
193  static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {}
194  static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {}
195 
196  template<typename AlphaType>
197  static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {}
198  };
199 
200 
201 
202 
203 
204 
213  template<typename T1, typename T2, typename T3, typename T4,
214  typename T5, typename T6, typename T7, typename T8, typename T9>
215  class MultiTypeBlockMatrix : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
216 
217  public:
218 
222  typedef MultiTypeBlockMatrix<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
223 
224  typedef typename mpl::at_c<T1,0>::type field_type;
225 
229  template<typename T>
230  void operator= (const T& newval) {MultiTypeBlockMatrix_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
231 
235  template<typename X, typename Y>
236  void mv (const X& x, Y& y) const {
237  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
238  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
239 
240  y = 0; //reset y (for mv uses umv)
241  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
242  }
243 
247  template<typename X, typename Y>
248  void umv (const X& x, Y& y) const {
249  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
250  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
251 
252  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
253  }
254 
258  template<typename X, typename Y>
259  void mmv (const X& x, Y& y) const {
260  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
261  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
262 
263  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::mmv(y, *this, x); //iterate over all matrix elements
264  }
265 
267  template<typename AlphaType, typename X, typename Y>
268  void usmv (const AlphaType& alpha, const X& x, Y& y) const {
269  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
270  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
271 
272  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::usmv(alpha,y, *this, x); //iterate over all matrix elements
273 
274  }
275 
276 
277 
278  };
279 
280 
281 
287  template<typename T1, typename T2, typename T3, typename T4, typename T5,
288  typename T6, typename T7, typename T8, typename T9>
289  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>& m) {
290  static const int i = mpl::size<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value; //row count
291  static const int j = mpl::size< typename mpl::at_c<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>,0>::type >::value; //col count of first row
292  MultiTypeBlockMatrix_Print<0,i,0,j,MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(m);
293  return s;
294  }
295 
296 
297 
298 
299 
300  //make algmeta_itsteps known
301  template<int I>
302  struct algmeta_itsteps;
303 
304 
305 
306 
307 
308 
315  template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
316  class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
317  public:
321  template <typename Trhs, typename TVector, typename TMatrix, typename K>
322  static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
323  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), b );
324  MultiTypeBlockMatrix_Solver_Col<I, crow, ccol+1, remain_col-1>::calc_rhs(A,x,v,b,w); //next column element
325  }
326 
327  };
328  template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
329  class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
330  public:
331  template <typename Trhs, typename TVector, typename TMatrix, typename K>
332  static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {}
333  };
334 
335 
336 
343  template<int I, int crow, int remain_row>
344  class MultiTypeBlockMatrix_Solver {
345  public:
346 
350  template <typename TVector, typename TMatrix, typename K>
351  static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
352  TVector xold(x);
353  xold=x; //store old x values
355  x *= w;
356  x.axpy(1-w,xold); //improve x
357  }
358  template <typename TVector, typename TMatrix, typename K>
359  static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
360  typename mpl::at_c<TVector,crow>::type rhs;
361  rhs = fusion::at_c<crow> (b);
362 
363  MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
364  //solve on blocklevel I-1
365  algmeta_itsteps<I-1>::dbgs(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(x),rhs,w);
367  }
368 
369 
370 
374  template <typename TVector, typename TMatrix, typename K>
375  static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
376  TVector v;
377  v=x; //use latest x values in right side calculation
379 
380  }
381  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
382  static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
383  typename mpl::at_c<TVector,crow>::type rhs;
384  rhs = fusion::at_c<crow> (b);
385 
386  MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
387  //solve on blocklevel I-1
388  algmeta_itsteps<I-1>::bsorf(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
389  fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
391  }
392 
396  template <typename TVector, typename TMatrix, typename K>
397  static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
398  TVector v;
399  v=x; //use latest x values in right side calculation
401 
402  }
403  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
404  static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
405  typename mpl::at_c<TVector,crow>::type rhs;
406  rhs = fusion::at_c<crow> (b);
407 
408  MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
409  //solve on blocklevel I-1
410  algmeta_itsteps<I-1>::bsorb(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
411  fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
413  }
414 
415 
419  template <typename TVector, typename TMatrix, typename K>
420  static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
421  TVector v(x);
422  v=0; //calc new x in v
424  x.axpy(w,v); //improve x
425  }
426  template <typename TVector, typename TMatrix, typename K>
427  static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
428  typename mpl::at_c<TVector,crow>::type rhs;
429  rhs = fusion::at_c<crow> (b);
430 
431  MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
432  //solve on blocklevel I-1
433  algmeta_itsteps<I-1>::dbjac(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
435  }
436 
437 
438 
439 
440  };
441  template<int I, int crow> //recursion end for remain_row = 0
442  class MultiTypeBlockMatrix_Solver<I,crow,0> {
443  public:
444  template <typename TVector, typename TMatrix, typename K>
445  static void dbgs(const TMatrix& A, TVector& x, TVector& v,
446  const TVector& b, const K& w) {}
447 
448  template <typename TVector, typename TMatrix, typename K>
449  static void bsorf(const TMatrix& A, TVector& x, TVector& v,
450  const TVector& b, const K& w) {}
451 
452  template <typename TVector, typename TMatrix, typename K>
453  static void bsorb(const TMatrix& A, TVector& x, TVector& v,
454  const TVector& b, const K& w) {}
455 
456  template <typename TVector, typename TMatrix, typename K>
457  static void dbjac(const TMatrix& A, TVector& x, TVector& v,
458  const TVector& b, const K& w) {}
459  };
460 
461 } // end namespace
462 
463 #endif // HAVE_BOOST_FUSION
464 #endif // HAVE_BOOST
465 #endif
466