3 #ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
4 #define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
12 #ifdef HAVE_BOOST_FUSION
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>
20 namespace mpl=boost::mpl;
21 namespace fusion=boost::fusion;
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;
31 template<
int I,
int crow,
int remain_row>
32 class MultiTypeBlockMatrix_Solver;
56 template<
int crow,
int remain_rows,
int ccol,
int remain_cols,
58 class MultiTypeBlockMatrix_Print {
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);
69 template<
int crow,
int remain_rows,
int ccol,
typename TMatrix>
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);
77 template<
int crow,
int ccol,
int remain_cols,
typename TMatrix>
78 class MultiTypeBlockMatrix_Print<crow,0,ccol,remain_cols,TMatrix> {
80 static void print(
const TMatrix& m)
81 {std::cout << std::endl;}
87 template<
int count,
typename T1,
typename T2>
88 class MultiTypeBlockVector_Ident;
103 template<
int rowcount,
typename T1,
typename T2>
104 class MultiTypeBlockMatrix_Ident {
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);
113 MultiTypeBlockMatrix_Ident<rowcount-1,T1,T2>::equalize(a,b);
118 template<
typename T1,
typename T2>
119 class MultiTypeBlockMatrix_Ident<0,T1,T2> {
121 static void equalize (T1&,
const T2&)
130 template<
int crow,
int remain_rows,
int ccol,
int remain_cols,
131 typename TVecY,
typename TMatrix,
typename TVecX>
132 class MultiTypeBlockMatrix_VectMul {
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);
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);
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);
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> {
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);
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);
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);
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> {
195 static void umv(TVecY&,
const TMatrix&,
const TVecX&) {}
196 static void mmv(TVecY&,
const TMatrix&,
const TVecX&) {}
198 template<
typename AlphaType>
199 static void usmv(
const AlphaType&, TVecY&,
const TMatrix&,
const TVecX&) {}
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> {
224 typedef MultiTypeBlockMatrix<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
226 typedef typename mpl::at_c<T1,0>::type field_type;
232 void operator= (
const T& newval) {MultiTypeBlockMatrix_Ident<mpl::size<type>::value,type,T>::equalize(*
this, newval); }
237 template<
typename X,
typename Y>
238 void mv (
const X& x, Y& y)
const {
239 BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value);
240 BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value);
243 MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *
this, x);
249 template<
typename X,
typename Y>
250 void umv (
const X& x, Y& y)
const {
251 BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value);
252 BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value);
254 MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *
this, x);
260 template<
typename X,
typename Y>
261 void mmv (
const X& x, Y& y)
const {
262 BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value);
263 BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value);
265 MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::mmv(y, *
this, x);
269 template<
typename AlphaType,
typename X,
typename Y>
270 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
271 BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value);
272 BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value);
274 MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::usmv(alpha,y, *
this, x);
289 template<
typename T1,
typename T2,
typename T3,
typename T4,
typename T5,
290 typename T6,
typename T7,
typename T8,
typename T9>
291 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>& m) {
292 static const int i = mpl::size<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value;
293 static const int j = mpl::size< typename mpl::at_c<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>,0>::type >::value;
294 MultiTypeBlockMatrix_Print<0,i,0,j,MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(m);
304 struct algmeta_itsteps;
317 template<
int I,
int crow,
int ccol,
int remain_col>
318 class MultiTypeBlockMatrix_Solver_Col {
323 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
324 static void calc_rhs(
const TMatrix&
A, TVector& x, TVector& v, Trhs& b,
const K& w) {
325 fusion::at_c<ccol>( fusion::at_c<crow>(
A) ).mmv( fusion::at_c<ccol>(x), b );
326 MultiTypeBlockMatrix_Solver_Col<I, crow, ccol+1, remain_col-1>::calc_rhs(
A,x,v,b,w);
330 template<
int I,
int crow,
int ccol>
331 class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
333 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
334 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
345 template<
int I,
int crow,
int remain_row>
346 class MultiTypeBlockMatrix_Solver {
352 template <
typename TVector,
typename TMatrix,
typename K>
353 static void dbgs(
const TMatrix&
A, TVector& x,
const TVector& b,
const K& w) {
360 template <
typename TVector,
typename TMatrix,
typename K>
361 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
362 typename mpl::at_c<TVector,crow>::type rhs;
363 rhs = fusion::at_c<crow> (b);
365 MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w);
376 template <
typename TVector,
typename TMatrix,
typename K>
377 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
383 template <
typename TVector,
typename TMatrix,
typename K>
384 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
385 typename mpl::at_c<TVector,crow>::type rhs;
386 rhs = fusion::at_c<crow> (b);
388 MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w);
391 fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
398 template <
typename TVector,
typename TMatrix,
typename K>
399 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
405 template <
typename TVector,
typename TMatrix,
typename K>
406 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
407 typename mpl::at_c<TVector,crow>::type rhs;
408 rhs = fusion::at_c<crow> (b);
410 MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w);
413 fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
421 template <
typename TVector,
typename TMatrix,
typename K>
422 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
428 template <
typename TVector,
typename TMatrix,
typename K>
429 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
430 typename mpl::at_c<TVector,crow>::type rhs;
431 rhs = fusion::at_c<crow> (b);
433 MultiTypeBlockMatrix_Solver_Col<I,crow,0, mpl::size<typename mpl::at_c<TMatrix,crow>::type>::value>::calc_rhs(A,x,v,rhs,w);
443 template<
int I,
int crow>
444 class MultiTypeBlockMatrix_Solver<I,crow,0> {
446 template <
typename TVector,
typename TMatrix,
typename K>
447 static void dbgs(
const TMatrix&, TVector&, TVector&,
448 const TVector&,
const K&) {}
450 template <
typename TVector,
typename TMatrix,
typename K>
451 static void bsorf(
const TMatrix&, TVector&, TVector&,
452 const TVector&,
const K&) {}
454 template <
typename TVector,
typename TMatrix,
typename K>
455 static void bsorb(
const TMatrix&, TVector&, TVector&,
456 const TVector&,
const K&) {}
458 template <
typename TVector,
typename TMatrix,
typename K>
459 static void dbjac(
const TMatrix&, TVector&, TVector&,
460 const TVector&,
const K&) {}
465 #endif // HAVE_BOOST_FUSION
466 #endif // HAVE_DUNE_BOOST
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:612
Definition: basearray.hh:19
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 dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:600
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:636
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:495
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.