dune-istl  3.0-git
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 #include <tuple>
9 
10 #include <dune/common/hybridutilities.hh>
11 
12 #include "istlexception.hh"
13 
14 // forward declaration
15 namespace Dune
16 {
17  template<typename FirstRow, typename... Args>
18  class MultiTypeBlockMatrix;
19 
20  template<int I, int crow, int remain_row>
22 }
23 
24 #include "gsetc.hh"
25 
26 namespace Dune {
27 
41  template<typename FirstRow, typename... Args>
43  : public std::tuple<FirstRow, Args...>
44  {
45  public:
46 
50  typedef MultiTypeBlockMatrix<FirstRow, Args...> type;
51 
52  typedef typename FirstRow::field_type field_type;
53 
55  static constexpr std::size_t N()
56  {
57  return 1+sizeof...(Args);
58  }
59 
61  static constexpr std::size_t M()
62  {
63  return FirstRow::size();
64  }
65 
82  template< std::size_t index >
83  auto
84  operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) -> decltype(std::get<index>(*this))
85  {
86  DUNE_UNUSED_PARAMETER(indexVariable);
87  return std::get<index>(*this);
88  }
89 
95  template< std::size_t index >
96  auto
97  operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const -> decltype(std::get<index>(*this))
98  {
99  DUNE_UNUSED_PARAMETER(indexVariable);
100  return std::get<index>(*this);
101  }
102 
106  template<typename T>
107  void operator= (const T& newval) {
108  using namespace Dune::Hybrid;
109  auto size = index_constant<1+sizeof...(Args)>();
110  // Since Dune::Hybrid::size(MultiTypeBlockMatrix) is not implemented,
111  // we cannot use a plain forEach(*this, ...). This could be achieved,
112  // e.g., by implementing a static size() method.
113  forEach(integralRange(size), [&](auto&& i) {
114  (*this)[i] = newval;
115  });
116  }
117 
120  template<typename X, typename Y>
121  void mv (const X& x, Y& y) const {
122  static_assert(X::size() == M(), "length of x does not match row length");
123  static_assert(Y::size() == N(), "length of y does not match row count");
124  y = 0; //reset y (for mv uses umv)
125  umv(x,y);
126  }
127 
130  template<typename X, typename Y>
131  void umv (const X& x, Y& y) const {
132  static_assert(X::size() == M(), "length of x does not match row length");
133  static_assert(Y::size() == N(), "length of y does not match row count");
134  using namespace Dune::Hybrid;
135  forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
136  forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
137  (*this)[i][j].umv(x[j], y[i]);
138  });
139  });
140  }
141 
144  template<typename X, typename Y>
145  void mmv (const X& x, Y& y) const {
146  static_assert(X::size() == M(), "length of x does not match row length");
147  static_assert(Y::size() == N(), "length of y does not match row count");
148  using namespace Dune::Hybrid;
149  forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
150  forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
151  (*this)[i][j].mmv(x[j], y[i]);
152  });
153  });
154  }
155 
158  template<typename AlphaType, typename X, typename Y>
159  void usmv (const AlphaType& alpha, const X& x, Y& y) const {
160  static_assert(X::size() == M(), "length of x does not match row length");
161  static_assert(Y::size() == N(), "length of y does not match row count");
162  using namespace Dune::Hybrid;
163  forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
164  forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
165  (*this)[i][j].usmv(alpha, x[j], y[i]);
166  });
167  });
168  }
169 
170  };
171 
172 
173 
179  template<typename T1, typename... Args>
180  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,Args...>& m) {
181  auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
182  auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
183  using namespace Dune::Hybrid;
184  forEach(integralRange(N), [&](auto&& i) {
185  forEach(integralRange(M), [&](auto&& j) {
186  s << "\t(" << i << ", " << j << "): \n" << m[i][j];
187  });
188  });
189  s << std::endl;
190  return s;
191  }
192 
193  //make algmeta_itsteps known
194  template<int I, typename M>
195  struct algmeta_itsteps;
196 
203  template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
204  class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
205  public:
209  template <typename Trhs, typename TVector, typename TMatrix, typename K>
210  static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
211  std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
213  }
214 
215  };
216  template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
217  class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
218  public:
219  template <typename Trhs, typename TVector, typename TMatrix, typename K>
220  static void calc_rhs(const TMatrix&, TVector&, TVector&, Trhs&, const K&) {}
221  };
222 
223 
224 
231  template<int I, int crow, int remain_row>
233  public:
234 
238  template <typename TVector, typename TMatrix, typename K>
239  static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
240  TVector xold(x);
241  xold=x; //store old x values
243  x *= w;
244  x.axpy(1-w,xold); //improve x
245  }
246  template <typename TVector, typename TMatrix, typename K>
247  static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
248  auto rhs = std::get<crow> (b);
249 
250  MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
251  //solve on blocklevel I-1
252  using M =
253  typename std::remove_cv<
254  typename std::remove_reference<
255  decltype(std::get<crow>( std::get<crow>(A)))
256  >::type
257  >::type;
258  algmeta_itsteps<I-1,M>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
260  }
261 
262 
263 
267  template <typename TVector, typename TMatrix, typename K>
268  static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
269  TVector v;
270  v=x; //use latest x values in right side calculation
272 
273  }
274  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
275  static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
276  auto rhs = std::get<crow> (b);
277 
278  MultiTypeBlockMatrix_Solver_Col<I,crow,0,TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
279  //solve on blocklevel I-1
280  using M =
281  typename std::remove_cv<
282  typename std::remove_reference<
283  decltype(std::get<crow>( std::get<crow>(A)))
284  >::type
285  >::type;
286  algmeta_itsteps<I-1,M>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
287  std::get<crow>(x).axpy(w,std::get<crow>(v));
289  }
290 
294  template <typename TVector, typename TMatrix, typename K>
295  static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
296  TVector v;
297  v=x; //use latest x values in right side calculation
299 
300  }
301  template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
302  static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
303  auto rhs = std::get<crow> (b);
304 
305  MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
306  //solve on blocklevel I-1
307  using M =
308  typename std::remove_cv<
309  typename std::remove_reference<
310  decltype(std::get<crow>( std::get<crow>(A)))
311  >::type
312  >::type;
313  algmeta_itsteps<I-1,M>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
314  std::get<crow>(x).axpy(w,std::get<crow>(v));
316  }
317 
318 
322  template <typename TVector, typename TMatrix, typename K>
323  static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
324  TVector v(x);
325  v=0; //calc new x in v
327  x.axpy(w,v); //improve x
328  }
329  template <typename TVector, typename TMatrix, typename K>
330  static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
331  auto rhs = std::get<crow> (b);
332 
333  MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
334  //solve on blocklevel I-1
335  using M =
336  typename std::remove_cv<
337  typename std::remove_reference<
338  decltype(std::get<crow>( std::get<crow>(A)))
339  >::type
340  >::type;
341  algmeta_itsteps<I-1,M>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
343  }
344 
345 
346 
347 
348  };
349  template<int I, int crow> //recursion end for remain_row = 0
350  class MultiTypeBlockMatrix_Solver<I,crow,0> {
351  public:
352  template <typename TVector, typename TMatrix, typename K>
353  static void dbgs(const TMatrix&, TVector&, TVector&,
354  const TVector&, const K&) {}
355 
356  template <typename TVector, typename TMatrix, typename K>
357  static void bsorf(const TMatrix&, TVector&, TVector&,
358  const TVector&, const K&) {}
359 
360  template <typename TVector, typename TMatrix, typename K>
361  static void bsorb(const TMatrix&, TVector&, TVector&,
362  const TVector&, const K&) {}
363 
364  template <typename TVector, typename TMatrix, typename K>
365  static void dbjac(const TMatrix&, TVector&, TVector&,
366  const TVector&, const K&) {}
367  };
368 
369 } // end namespace
370 
371 #endif
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:131
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:220
A Matrix class to support different block types.
Definition: matrixutils.hh:121
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:302
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:330
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:357
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:239
Definition: gsetc.hh:370
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:361
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:365
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:247
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:210
auto operator[](const std::integral_constant< std::size_t, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:84
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:295
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:275
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:52
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:204
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:159
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:373
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:461
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:401
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:107
Definition: basearray.hh:19
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:145
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:268
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:431
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:121
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:632
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:323
static constexpr std::size_t M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:61
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:50
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:21
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:353