Dune Core Modules (2.4.1)

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
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
20namespace mpl=boost::mpl;
21namespace fusion=boost::fusion;
22
23// forward decl
24namespace 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
37namespace 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
std::ostream & operator<<(std::ostream &s, const array< T, N > &e)
Output operator for array.
Definition: array.hh:26
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:624
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:636
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:600
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:612
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignment.hh:10
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional 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 13, 23:29, 2024)