dune-istl  2.4.1
multitypeblockmatrix.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_MULTITYPEBLOCKMATRIX_HH
4 #define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
5 
6 #include <cmath>
7 #include <iostream>
8 
9 #include "istlexception.hh"
10 
11 #if HAVE_DUNE_BOOST
12 #ifdef HAVE_BOOST_FUSION
13 
14 #include <boost/fusion/sequence.hpp>
15 #include <boost/fusion/container.hpp>
16 #include <boost/fusion/iterator.hpp>
17 #include <boost/typeof/typeof.hpp>
18 #include <boost/fusion/algorithm.hpp>
19 
20 namespace mpl=boost::mpl;
21 namespace fusion=boost::fusion;
22 
23 // forward decl
24 namespace Dune
25 {
26  template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
27  typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
28  typename T8=fusion::void_, typename T9=fusion::void_>
29  class MultiTypeBlockMatrix;
30 
31  template<int I, int crow, int remain_row>
32  class MultiTypeBlockMatrix_Solver;
33 }
34 
35 #include "gsetc.hh"
36 
37 namespace Dune {
38 
56  template<int crow, int remain_rows, int ccol, int remain_cols,
57  typename TMatrix>
58  class MultiTypeBlockMatrix_Print {
59  public:
60 
64  static void print(const TMatrix& m) {
65  std::cout << "\t(" << crow << ", " << ccol << "): \n" << fusion::at_c<ccol>( fusion::at_c<crow>(m));
66  MultiTypeBlockMatrix_Print<crow,remain_rows,ccol+1,remain_cols-1,TMatrix>::print(m); //next column
67  }
68  };
69  template<int crow, int remain_rows, int ccol, typename TMatrix> //specialization for remain_cols=0
70  class MultiTypeBlockMatrix_Print<crow,remain_rows,ccol,0,TMatrix> {
71  public: static void print(const TMatrix& m) {
72  static const int xlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
73  MultiTypeBlockMatrix_Print<crow+1,remain_rows-1,0,xlen,TMatrix>::print(m); //next row
74  }
75  };
76 
77  template<int crow, int ccol, int remain_cols, typename TMatrix> //recursion end: specialization for remain_rows=0
78  class MultiTypeBlockMatrix_Print<crow,0,ccol,remain_cols,TMatrix> {
79  public:
80  static void print(const TMatrix& m)
81  {std::cout << std::endl;}
82  };
83 
84 
85 
86  //make MultiTypeBlockVector_Ident known (for MultiTypeBlockMatrix_Ident)
87  template<int count, typename T1, typename T2>
88  class MultiTypeBlockVector_Ident;
89 
90 
103  template<int rowcount, typename T1, typename T2>
104  class MultiTypeBlockMatrix_Ident {
105  public:
106 
111  static void equalize(T1& a, const T2& b) {
112  MultiTypeBlockVector_Ident< mpl::size< typename mpl::at_c<T1,rowcount-1>::type >::value ,T1,T2>::equalize(a,b); //rows are cvectors
113  MultiTypeBlockMatrix_Ident<rowcount-1,T1,T2>::equalize(a,b); //iterate over rows
114  }
115  };
116 
117  //recursion end for rowcount=0
118  template<typename T1, typename T2>
119  class MultiTypeBlockMatrix_Ident<0,T1,T2> {
120  public:
121  static void equalize (T1&, const T2&)
122  {}
123  };
124 
130  template<int crow, int remain_rows, int ccol, int remain_cols,
131  typename TVecY, typename TMatrix, typename TVecX>
132  class MultiTypeBlockMatrix_VectMul {
133  public:
134 
138  static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
139  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).umv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
140  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::umv(y, A, x);
141  }
142 
146  static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
147  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
148  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::mmv(y, A, x);
149  }
150 
151  template<typename AlphaType>
152  static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
153  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).usmv(alpha, fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
154  MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
155  }
156 
157 
158  };
159 
160  //specialization for remain_cols = 0
161  template<int crow, int remain_rows,int ccol, typename TVecY,
162  typename TMatrix, typename TVecX>
163  class MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol,0,TVecY,TMatrix,TVecX> { //start iteration over next row
164 
165  public:
169  static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
170  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
171  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::umv(y, A, x);
172  }
173 
177  static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
178  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
179  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::mmv(y, A, x);
180  }
181 
182  template <typename AlphaType>
183  static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
184  static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
185  MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
186  }
187  };
188 
189  //specialization for remain_rows = 0
190  template<int crow, int ccol, int remain_cols, typename TVecY,
191  typename TMatrix, typename TVecX>
192  class MultiTypeBlockMatrix_VectMul<crow,0,ccol,remain_cols,TVecY,TMatrix,TVecX> {
193  //end recursion
194  public:
195  static void umv(TVecY&, const TMatrix&, const TVecX&) {}
196  static void mmv(TVecY&, const TMatrix&, const TVecX&) {}
197 
198  template<typename AlphaType>
199  static void usmv(const AlphaType&, TVecY&, const TMatrix&, const TVecX&) {}
200  };
201 
202 
203 
204 
205 
206 
215  template<typename T1, typename T2, typename T3, typename T4,
216  typename T5, typename T6, typename T7, typename T8, typename T9>
217  class MultiTypeBlockMatrix : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
218 
219  public:
220 
224  typedef MultiTypeBlockMatrix<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
225 
226  typedef typename mpl::at_c<T1,0>::type field_type;
227 
244  template< std::size_t index >
245  auto
246  operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) -> decltype(fusion::at_c<index>(*this))
247  {
248  DUNE_UNUSED_PARAMETER(indexVariable);
249  return fusion::at_c<index>(*this);
250  }
251 
257  template< std::size_t index >
258  auto
259  operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const -> decltype(fusion::at_c<index>(*this))
260  {
261  DUNE_UNUSED_PARAMETER(indexVariable);
262  return fusion::at_c<index>(*this);
263  }
264 
265 
269  template<typename T>
270  void operator= (const T& newval) {MultiTypeBlockMatrix_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
271 
275  template<typename X, typename Y>
276  void mv (const X& x, Y& y) const {
277  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
278  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
279 
280  y = 0; //reset y (for mv uses umv)
281  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
282  }
283 
287  template<typename X, typename Y>
288  void umv (const X& x, Y& y) const {
289  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
290  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
291 
292  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
293  }
294 
298  template<typename X, typename Y>
299  void mmv (const X& x, Y& y) const {
300  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
301  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
302 
303  MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::mmv(y, *this, x); //iterate over all matrix elements
304  }
305 
307  template<typename AlphaType, typename X, typename Y>
308  void usmv (const AlphaType& alpha, const X& x, Y& y) const {
309  BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
310  BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
311 
312  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
313 
314  }
315 
316 
317 
318  };
319 
320 
321 
327  template<typename T1, typename T2, typename T3, typename T4, typename T5,
328  typename T6, typename T7, typename T8, typename T9>
329  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>& m) {
330  static const int i = mpl::size<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value; //row count
331  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
332  MultiTypeBlockMatrix_Print<0,i,0,j,MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(m);
333  return s;
334  }
335 
336 
337 
338 
339 
340  //make algmeta_itsteps known
341  template<int I>
342  struct algmeta_itsteps;
343 
344 
345 
346 
347 
348 
355  template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
356  class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
357  public:
361  template <typename Trhs, typename TVector, typename TMatrix, typename K>
362  static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
363  fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), b );
364  MultiTypeBlockMatrix_Solver_Col<I, crow, ccol+1, remain_col-1>::calc_rhs(A,x,v,b,w); //next column element
365  }
366 
367  };
368  template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
369  class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
370  public:
371  template <typename Trhs, typename TVector, typename TMatrix, typename K>
372  static void calc_rhs(const TMatrix&, TVector&, TVector&, Trhs&, const K&) {}
373  };
374 
375 
376 
383  template<int I, int crow, int remain_row>
384  class MultiTypeBlockMatrix_Solver {
385  public:
386 
390  template <typename TVector, typename TMatrix, typename K>
391  static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
392  TVector xold(x);
393  xold=x; //store old x values
395  x *= w;
396  x.axpy(1-w,xold); //improve x
397  }
398  template <typename TVector, typename TMatrix, typename K>
399  static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
400  typename mpl::at_c<TVector,crow>::type rhs;
401  rhs = fusion::at_c<crow> (b);
402 
403  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
404  //solve on blocklevel I-1
405  algmeta_itsteps<I-1>::dbgs(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(x),rhs,w);
407  }
408 
409 
410 
414  template <typename TVector, typename TMatrix, typename K>
415  static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
416  TVector v;
417  v=x; //use latest x values in right side calculation
419 
420  }
421  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
422  static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
423  typename mpl::at_c<TVector,crow>::type rhs;
424  rhs = fusion::at_c<crow> (b);
425 
426  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
427  //solve on blocklevel I-1
428  algmeta_itsteps<I-1>::bsorf(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
429  fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
431  }
432 
436  template <typename TVector, typename TMatrix, typename K>
437  static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
438  TVector v;
439  v=x; //use latest x values in right side calculation
441 
442  }
443  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
444  static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
445  typename mpl::at_c<TVector,crow>::type rhs;
446  rhs = fusion::at_c<crow> (b);
447 
448  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
449  //solve on blocklevel I-1
450  algmeta_itsteps<I-1>::bsorb(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
451  fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
453  }
454 
455 
459  template <typename TVector, typename TMatrix, typename K>
460  static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
461  TVector v(x);
462  v=0; //calc new x in v
464  x.axpy(w,v); //improve x
465  }
466  template <typename TVector, typename TMatrix, typename K>
467  static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
468  typename mpl::at_c<TVector,crow>::type rhs;
469  rhs = fusion::at_c<crow> (b);
470 
471  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
472  //solve on blocklevel I-1
473  algmeta_itsteps<I-1>::dbjac(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
475  }
476 
477 
478 
479 
480  };
481  template<int I, int crow> //recursion end for remain_row = 0
482  class MultiTypeBlockMatrix_Solver<I,crow,0> {
483  public:
484  template <typename TVector, typename TMatrix, typename K>
485  static void dbgs(const TMatrix&, TVector&, TVector&,
486  const TVector&, const K&) {}
487 
488  template <typename TVector, typename TMatrix, typename K>
489  static void bsorf(const TMatrix&, TVector&, TVector&,
490  const TVector&, const K&) {}
491 
492  template <typename TVector, typename TMatrix, typename K>
493  static void bsorb(const TMatrix&, TVector&, TVector&,
494  const TVector&, const K&) {}
495 
496  template <typename TVector, typename TMatrix, typename K>
497  static void dbjac(const TMatrix&, TVector&, TVector&,
498  const TVector&, const K&) {}
499  };
500 
501 } // end namespace
502 
503 #endif // HAVE_BOOST_FUSION
504 #endif // HAVE_DUNE_BOOST
505 #endif
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:612
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:545
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:624
Matrix & A
Definition: matrixmatrix.hh:216
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:397
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:445
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:636
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Definition: basearray.hh:19
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:495
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:600