Dune Core Modules (2.5.0)

multitypeblockmatrix.hh
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
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 size()
62 {
63 return 1+sizeof...(Args);
64 }
65
67 static constexpr std::size_t M()
68 {
69 return FirstRow::size();
70 }
71
88 template< std::size_t index >
89 auto
90 operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) -> decltype(std::get<index>(*this))
91 {
92 DUNE_UNUSED_PARAMETER(indexVariable);
93 return std::get<index>(*this);
94 }
95
101 template< std::size_t index >
102 auto
103 operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const -> decltype(std::get<index>(*this))
104 {
105 DUNE_UNUSED_PARAMETER(indexVariable);
106 return std::get<index>(*this);
107 }
108
112 template<typename T>
113 void operator= (const T& newval) {
114 using namespace Dune::Hybrid;
115 auto size = index_constant<1+sizeof...(Args)>();
116 // Since Dune::Hybrid::size(MultiTypeBlockMatrix) is not implemented,
117 // we cannot use a plain forEach(*this, ...). This could be achieved,
118 // e.g., by implementing a static size() method.
119 forEach(integralRange(size), [&](auto&& i) {
120 (*this)[i] = newval;
121 });
122 }
123
126 template<typename X, typename Y>
127 void mv (const X& x, Y& y) const {
128 static_assert(X::size() == M(), "length of x does not match row length");
129 static_assert(Y::size() == N(), "length of y does not match row count");
130 y = 0; //reset y (for mv uses umv)
131 umv(x,y);
132 }
133
136 template<typename X, typename Y>
137 void umv (const X& x, Y& y) const {
138 static_assert(X::size() == M(), "length of x does not match row length");
139 static_assert(Y::size() == N(), "length of y does not match row count");
140 using namespace Dune::Hybrid;
141 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
142 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
143 (*this)[i][j].umv(x[j], y[i]);
144 });
145 });
146 }
147
150 template<typename X, typename Y>
151 void mmv (const X& x, Y& y) const {
152 static_assert(X::size() == M(), "length of x does not match row length");
153 static_assert(Y::size() == N(), "length of y does not match row count");
154 using namespace Dune::Hybrid;
155 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
156 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
157 (*this)[i][j].mmv(x[j], y[i]);
158 });
159 });
160 }
161
164 template<typename AlphaType, typename X, typename Y>
165 void usmv (const AlphaType& alpha, const X& x, Y& y) const {
166 static_assert(X::size() == M(), "length of x does not match row length");
167 static_assert(Y::size() == N(), "length of y does not match row count");
168 using namespace Dune::Hybrid;
169 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
170 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
171 (*this)[i][j].usmv(alpha, x[j], y[i]);
172 });
173 });
174 }
175
176 };
177
178
179
185 template<typename T1, typename... Args>
186 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,Args...>& m) {
189 using namespace Dune::Hybrid;
190 forEach(integralRange(N), [&](auto&& i) {
191 forEach(integralRange(M), [&](auto&& j) {
192 s << "\t(" << i << ", " << j << "): \n" << m[i][j];
193 });
194 });
195 s << std::endl;
196 return s;
197 }
198
199 //make algmeta_itsteps known
200 template<int I, typename M>
201 struct algmeta_itsteps;
202
209 template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
210 class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
211 public:
215 template <typename Trhs, typename TVector, typename TMatrix, typename K>
216 static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
217 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
219 }
220
221 };
222 template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
223 class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
224 public:
225 template <typename Trhs, typename TVector, typename TMatrix, typename K>
226 static void calc_rhs(const TMatrix&, TVector&, TVector&, Trhs&, const K&) {}
227 };
228
229
230
237 template<int I, int crow, int remain_row>
239 public:
240
244 template <typename TVector, typename TMatrix, typename K>
245 static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
246 TVector xold(x);
247 xold=x; //store old x values
249 x *= w;
250 x.axpy(1-w,xold); //improve x
251 }
252 template <typename TVector, typename TMatrix, typename K>
253 static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
254 auto rhs = std::get<crow> (b);
255
256 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
257 //solve on blocklevel I-1
258 using M =
259 typename std::remove_cv<
260 typename std::remove_reference<
261 decltype(std::get<crow>( std::get<crow>(A)))
262 >::type
263 >::type;
264 algmeta_itsteps<I-1,M>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
266 }
267
268
269
273 template <typename TVector, typename TMatrix, typename K>
274 static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
275 TVector v;
276 v=x; //use latest x values in right side calculation
278
279 }
280 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
281 static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
282 auto rhs = std::get<crow> (b);
283
284 MultiTypeBlockMatrix_Solver_Col<I,crow,0,TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
285 //solve on blocklevel I-1
286 using M =
287 typename std::remove_cv<
288 typename std::remove_reference<
289 decltype(std::get<crow>( std::get<crow>(A)))
290 >::type
291 >::type;
292 algmeta_itsteps<I-1,M>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
293 std::get<crow>(x).axpy(w,std::get<crow>(v));
295 }
296
300 template <typename TVector, typename TMatrix, typename K>
301 static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
302 TVector v;
303 v=x; //use latest x values in right side calculation
305
306 }
307 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
308 static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
309 auto rhs = std::get<crow> (b);
310
311 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
312 //solve on blocklevel I-1
313 using M =
314 typename std::remove_cv<
315 typename std::remove_reference<
316 decltype(std::get<crow>( std::get<crow>(A)))
317 >::type
318 >::type;
319 algmeta_itsteps<I-1,M>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
320 std::get<crow>(x).axpy(w,std::get<crow>(v));
322 }
323
324
328 template <typename TVector, typename TMatrix, typename K>
329 static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
330 TVector v(x);
331 v=0; //calc new x in v
333 x.axpy(w,v); //improve x
334 }
335 template <typename TVector, typename TMatrix, typename K>
336 static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
337 auto rhs = std::get<crow> (b);
338
339 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
340 //solve on blocklevel I-1
341 using M =
342 typename std::remove_cv<
343 typename std::remove_reference<
344 decltype(std::get<crow>( std::get<crow>(A)))
345 >::type
346 >::type;
347 algmeta_itsteps<I-1,M>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
349 }
350
351
352
353
354 };
355 template<int I, int crow> //recursion end for remain_row = 0
356 class MultiTypeBlockMatrix_Solver<I,crow,0> {
357 public:
358 template <typename TVector, typename TMatrix, typename K>
359 static void dbgs(const TMatrix&, TVector&, TVector&,
360 const TVector&, const K&) {}
361
362 template <typename TVector, typename TMatrix, typename K>
363 static void bsorf(const TMatrix&, TVector&, TVector&,
364 const TVector&, const K&) {}
365
366 template <typename TVector, typename TMatrix, typename K>
367 static void bsorb(const TMatrix&, TVector&, TVector&,
368 const TVector&, const K&) {}
369
370 template <typename TVector, typename TMatrix, typename K>
371 static void dbjac(const TMatrix&, TVector&, TVector&,
372 const TVector&, const K&) {}
373 };
374
375} // end namespace
376
377#endif
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:210
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:238
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:26
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:314
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:227
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:591
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:603
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:567
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:579
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:329
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:127
static constexpr std::size_t size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:61
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:245
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:301
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:165
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:113
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:137
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:216
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:274
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
auto operator[](const std::integral_constant< std::size_t, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:90
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:151
static constexpr std::size_t M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:67
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignment.hh:11
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)