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
15namespace Dune
16{
17 template<typename FirstRow, typename... Args>
18 class MultiTypeBlockMatrix;
19
20 template<int I, int crow, int remain_row>
21 class MultiTypeBlockMatrix_Solver;
22}
23
24#include "gsetc.hh"
25
26namespace Dune {
27
41 template<typename FirstRow, typename... Args>
43 : public std::tuple<FirstRow, Args...>
44 {
45 public:
46
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 {
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 {
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) {
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
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
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
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition multitypeblockmatrix.hh:50
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:323
void mv(const X &x, Y &y) const
y = A x
Definition multitypeblockmatrix.hh:121
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:353
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:239
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:295
FirstRow::field_type field_type
Definition multitypeblockmatrix.hh:52
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition multitypeblockmatrix.hh:159
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:247
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:275
void operator=(const T &newval)
Definition multitypeblockmatrix.hh:107
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition multitypeblockmatrix.hh:220
void umv(const X &x, Y &y) const
y += A x
Definition multitypeblockmatrix.hh:131
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition multitypeblockmatrix.hh:210
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:268
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:365
static constexpr std::size_t N()
Return the number of matrix rows.
Definition multitypeblockmatrix.hh:55
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:361
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 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
void mmv(const X &x, Y &y) const
y -= A x
Definition multitypeblockmatrix.hh:145
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:302
static constexpr std::size_t M()
Return the number of matrix columns.
Definition multitypeblockmatrix.hh:61
Definition basearray.hh:19
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition bvector.hh:632
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:431
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:401
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:461
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:373
A Matrix class to support different block types.
Definition multitypeblockmatrix.hh:44
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition multitypeblockmatrix.hh:232
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition multitypeblockmatrix.hh:204