diff --git a/dune/istl/Makefile.am b/dune/istl/Makefile.am
index 1e466dd7ec26a4fa28b20f1a669b9c6248cf065d..259915881a38446ada98958f40091d6fe9a93000 100644
--- a/dune/istl/Makefile.am
+++ b/dune/istl/Makefile.am
@@ -25,6 +25,8 @@ istl_HEADERS = allocator.hh \
matrixmatrix.hh \
matrixredistribute.hh \
matrixutils.hh \
+ multitypeblockmatrix.hh \
+ multitypeblockvector.hh \
operators.hh \
overlappingschwarz.hh \
owneroverlapcopy.hh \
diff --git a/dune/istl/gsetc.hh b/dune/istl/gsetc.hh
index 205bb95aaf0f66072581c45a94edd17e1e954534..18055349477ce2c3e6f5f9ec6e16f21679c5999e 100644
--- a/dune/istl/gsetc.hh
+++ b/dune/istl/gsetc.hh
@@ -8,6 +8,8 @@
#include <iostream>
#include <iomanip>
#include <string>
+#include "cvector.hh"
+#include "cmatrix.hh"
#include "istlexception.hh"
@@ -369,6 +371,25 @@ namespace Dune {
// template meta program for iterative solver steps
template<int I>
struct algmeta_itsteps {
+
+#ifdef HAVE_BOOST_FUSION
+
+ template<typename T11, typename T12, typename T13, typename T14,
+ typename T15, typename T16, typename T17, typename T18,
+ typename T19, typename T21, typename T22, typename T23,
+ typename T24, typename T25, typename T26, typename T27,
+ typename T28, typename T29, class K>
+ static void dbgs (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
+ MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
+ const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
+ const K& w)
+ {
+ const int rowcount =
+ boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
+ Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount>::dbgs(A, x, b, w);
+ }
+#endif
+
template<class M, class X, class Y, class K>
static void dbgs (const M& A, X& x, const Y& b, const K& w)
{
@@ -396,6 +417,24 @@ namespace Dune {
x *= w;
x.axpy(1-w,xold);
}
+
+#ifdef HAVE_BOOST_FUSION
+ template<typename T11, typename T12, typename T13, typename T14,
+ typename T15, typename T16, typename T17, typename T18,
+ typename T19, typename T21, typename T22, typename T23,
+ typename T24, typename T25, typename T26, typename T27,
+ typename T28, typename T29, class K>
+ static void bsorf (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
+ MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
+ const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
+ const K& w)
+ {
+ const int rowcount =
+ boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
+ Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount>::bsorf(A, x, b, w);
+ }
+#endif
+
template<class M, class X, class Y, class K>
static void bsorf (const M& A, X& x, const Y& b, const K& w)
{
@@ -425,6 +464,25 @@ namespace Dune {
x[i.index()].axpy(w,v);
}
}
+
+#ifdef HAVE_BOOST_FUSION
+
+ template<typename T11, typename T12, typename T13, typename T14,
+ typename T15, typename T16, typename T17, typename T18,
+ typename T19, typename T21, typename T22, typename T23,
+ typename T24, typename T25, typename T26, typename T27,
+ typename T28, typename T29, class K>
+ static void bsorb (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
+ MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
+ const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
+ const K& w)
+ {
+ const int rowcount =
+ mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
+ Dune::MultiTypeBlockMatrix_Solver<I,rowcount-1,rowcount>::bsorb(A, x, b, w);
+ }
+#endif
+
template<class M, class X, class Y, class K>
static void bsorb (const M& A, X& x, const Y& b, const K& w)
{
@@ -454,6 +512,24 @@ namespace Dune {
x[i.index()].axpy(w,v);
}
}
+
+#ifdef HAVE_BOOST_FUSION
+ template<typename T11, typename T12, typename T13, typename T14,
+ typename T15, typename T16, typename T17, typename T18,
+ typename T19, typename T21, typename T22, typename T23,
+ typename T24, typename T25, typename T26, typename T27,
+ typename T28, typename T29, class K>
+ static void dbjac (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
+ MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
+ const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
+ const K& w)
+ {
+ const int rowcount =
+ boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value
+ Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount >::dbjac(A, x, b, w);
+ }
+#endif
+
template<class M, class X, class Y, class K>
static void dbjac (const M& A, X& x, const Y& b, const K& w)
{
diff --git a/dune/istl/multitypeblockmatrix.hh b/dune/istl/multitypeblockmatrix.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7bce230ca3e1fdb34dd37e140c23246ff79c3dcd
--- /dev/null
+++ b/dune/istl/multitypeblockmatrix.hh
@@ -0,0 +1,455 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MultiTypeMATRIX_HH
+#define DUNE_MultiTypeMATRIX_HH
+
+#ifdef HAVE_BOOST_FUSION
+
+#include <cmath>
+#include <iostream>
+
+#include "istlexception.hh"
+#include "allocator.hh"
+
+#include "gsetc.hh"
+
+#include <boost/fusion/sequence.hpp>
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/iterator.hpp>
+#include <boost/typeof/typeof.hpp>
+#include <boost/fusion/algorithm.hpp>
+
+#include "cvector.hh"
+
+namespace mpl=boost::mpl;
+namespace fusion=boost::fusion;
+
+namespace Dune {
+
+ /**
+ @addtogroup ISTL_SPMV
+ @{
+ */
+
+
+
+
+ /**
+ @brief prints out a matrix (type MultiTypeBlockMatrix)
+
+ The parameters "crow" and "ccol" are the indices of
+ the element to be printed out. Via recursive calling
+ all other elements are printed, too. This is internally
+ called when a MultiTypeBlockMatrix object is passed to an output stream.
+ */
+
+ template<int crow, int remain_rows, int ccol, int remain_cols,
+ typename TMatrix>
+ class MultiTypeBlockMatrix_Print {
+ public:
+
+ /**
+ * print out a matrix block and all following
+ */
+ static void print(const TMatrix& m) {
+ std::cout << "\t(" << crow << ", " << ccol << "): \n" << fusion::at_c<ccol>( fusion::at_c<crow>(m));
+ MultiTypeBlockMatrix_Print<crow,remain_rows,ccol+1,remain_cols-1,TMatrix>::print(m); //next column
+ }
+ };
+ template<int crow, int remain_rows, int ccol, typename TMatrix> //specialization for remain_cols=0
+ class MultiTypeBlockMatrix_Print<crow,remain_rows,ccol,0,TMatrix> {
+ public: static void print(const TMatrix& m) {
+ static const int xlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
+ MultiTypeBlockMatrix_Print<crow+1,remain_rows-1,0,xlen,TMatrix>::print(m); //next row
+ }
+ };
+
+ template<int crow, int ccol, int remain_cols, typename TMatrix> //recursion end: specialization for remain_rows=0
+ class MultiTypeBlockMatrix_Print<crow,0,ccol,remain_cols,TMatrix> {
+ public:
+ static void print(const TMatrix& m)
+ {std::cout << std::endl;}
+ };
+
+
+
+ //make MultiTypeBlockVector_Ident known (for MultiTypeBlockMatrix_Ident)
+ template<int count, typename T1, typename T2>
+ class MultiTypeBlockVector_Ident;
+
+
+ /**
+ @brief set a MultiTypeBlockMatrix to some specific value
+
+ This class is used by the MultiTypeBlockMatrix class' internal = operator.
+ Whenever a vector is assigned to a value, each block element
+ has to provide the = operator for the right side. Example:
+ \code
+ typedef MultiTypeBlockVector<int,int,int> CVect;
+ MultiTypeBlockMatrix<CVect,CVect> CMat;
+ CMat = 3; //sets all 3*2 integer elements to 3
+ \endcode
+ */
+ template<int rowcount, typename T1, typename T2>
+ class MultiTypeBlockMatrix_Ident {
+ public:
+
+ /**
+ * equalize two matrix' element
+ * note: uses MultiTypeBlockVector_Ident to equalize each row (which is of MultiTypeBlockVector type)
+ */
+ static void equalize(T1& a, const T2& b) {
+ MultiTypeBlockVector_Ident< mpl::size< typename mpl::at_c<T1,rowcount-1>::type >::value ,T1,T2>::equalize(a,b); //rows are cvectors
+ MultiTypeBlockMatrix_Ident<rowcount-1,T1,T2>::equalize(a,b); //iterate over rows
+ }
+ };
+
+ //recursion end for rowcount=0
+ template<typename T1, typename T2>
+ class MultiTypeBlockMatrix_Ident<0,T1,T2> {
+ public:
+ static void equalize (T1& a, const T2& b)
+ {}
+ };
+
+ /**
+ @brief Matrix Vector Multiplication
+
+ This class implements matrix vector multiplication for MultiTypeBlockMatrix/MultiTypeBlockVector types
+ */
+ template<int crow, int remain_rows, int ccol, int remain_cols,
+ typename TVecY, typename TMatrix, typename TVecX>
+ class MultiTypeBlockMatrix_VectMul {
+ public:
+
+ /**
+ * y += A x
+ */
+ static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
+ fusion::at_c<ccol>( fusion::at_c<crow>(A) ).umv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
+ MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::umv(y, A, x);
+ }
+
+ /**
+ * y -= A x
+ */
+ static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
+ fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
+ MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::mmv(y, A, x);
+ }
+
+ template<typename AlphaType>
+ static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
+ fusion::at_c<ccol>( fusion::at_c<crow>(A) ).usmv(alpha, fusion::at_c<ccol>(x), fusion::at_c<crow>(y) );
+ MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
+ }
+
+
+ };
+
+ //specialization for remain_cols = 0
+ template<int crow, int remain_rows,int ccol, typename TVecY,
+ typename TMatrix, typename TVecX>
+ class MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol,0,TVecY,TMatrix,TVecX> { //start iteration over next row
+
+ public:
+ /**
+ * do y -= A x in next row
+ */
+ static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
+ static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
+ MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::umv(y, A, x);
+ }
+
+ /**
+ * do y -= A x in next row
+ */
+ static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
+ static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
+ MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::mmv(y, A, x);
+ }
+
+ template <typename AlphaType>
+ static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
+ static const int rowlen = mpl::size< typename mpl::at_c<TMatrix,crow>::type >::value;
+ MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,rowlen,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
+ }
+ };
+
+ //specialization for remain_rows = 0
+ template<int crow, int ccol, int remain_cols, typename TVecY,
+ typename TMatrix, typename TVecX>
+ class MultiTypeBlockMatrix_VectMul<crow,0,ccol,remain_cols,TVecY,TMatrix,TVecX> {
+ //end recursion
+ public:
+ static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {}
+ static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {}
+
+ template<typename AlphaType>
+ static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {}
+ };
+
+
+
+
+
+
+ /**
+ @brief A Matrix class to support different block types
+
+ This matrix class combines MultiTypeBlockVector elements as rows.
+ */
+ template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
+ typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_, typename T8=fusion::void_, typename T9=fusion::void_>
+ class MultiTypeBlockMatrix : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
+
+ public:
+
+ /**
+ * own class' type
+ */
+ typedef MultiTypeBlockMatrix<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
+
+ typedef typename mpl::at_c<T1,0>::type field_type;
+
+ /**
+ * assignment operator
+ */
+ template<typename T>
+ void operator= (const T& newval) {MultiTypeBlockMatrix_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
+
+ /**
+ * y = A x
+ */
+ template<typename X, typename Y>
+ void mv (const X& x, Y& y) const {
+ BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
+ BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
+
+ y = 0; //reset y (for mv uses umv)
+ MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
+ }
+
+ /**
+ * y += A x
+ */
+ template<typename X, typename Y>
+ void umv (const X& x, Y& y) const {
+ BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
+ BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
+
+ MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
+ }
+
+ /**
+ * y -= A x
+ */
+ template<typename X, typename Y>
+ void mmv (const X& x, Y& y) const {
+ BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
+ BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
+
+ MultiTypeBlockMatrix_VectMul<0,mpl::size<type>::value,0,mpl::size<T1>::value,Y,type,X>::mmv(y, *this, x); //iterate over all matrix elements
+ }
+
+ //! y += alpha A x
+ template<typename AlphaType, typename X, typename Y>
+ void usmv (const AlphaType& alpha, const X& x, Y& y) const {
+ BOOST_STATIC_ASSERT(mpl::size<X>::value == mpl::size<T1>::value); //make sure x's length matches row length
+ BOOST_STATIC_ASSERT(mpl::size<Y>::value == mpl::size<type>::value); //make sure y's length matches row count
+
+ 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
+
+ }
+
+
+
+ };
+
+
+
+ /**
+ @brief << operator for a MultiTypeBlockMatrix
+
+ operator<< for printing out a MultiTypeBlockMatrix
+ */
+ template<typename T1, typename T2, typename T3, typename T4, typename T5,
+ typename T6, typename T7, typename T8, typename T9>
+ std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9>& m) {
+ static const int i = mpl::size<MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value; //row count
+ 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
+ MultiTypeBlockMatrix_Print<0,i,0,j,MultiTypeBlockMatrix<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(m);
+ return s;
+ }
+
+
+
+
+
+ //make algmeta_itsteps known
+ template<int I>
+ struct algmeta_itsteps;
+
+
+
+
+
+
+ /**
+ @brief part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
+
+ For the given row (index "crow") each element is used to
+ calculate the equation's right side.
+ */
+ template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
+ class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
+ public:
+ /**
+ * iterating over one row in MultiTypeBlockMatrix to calculate right side b-A[i][j]*x[j]
+ */
+ template <typename Trhs, typename TVector, typename TMatrix, typename K>
+ static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
+ fusion::at_c<ccol>( fusion::at_c<crow>(A) ).mmv( fusion::at_c<ccol>(x), b );
+ MultiTypeBlockMatrix_Solver_Col<I, crow, ccol+1, remain_col-1>::calc_rhs(A,x,v,b,w); //next column element
+ }
+
+ };
+ template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
+ class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
+ public:
+ template <typename Trhs, typename TVector, typename TMatrix, typename K>
+ static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {}
+ };
+
+
+
+ /**
+ @brief solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
+
+ The methods of this class are called by the solver specializations
+ for MultiTypeBlockVector & MultiTypeBlockMatrix types (dbgs, bsorf, bsorb, dbjac).
+ */
+ template<int I, int crow, int remain_row>
+ class MultiTypeBlockMatrix_Solver {
+ public:
+
+ /**
+ * Gauss-Seidel solver for MultiTypeBlockMatrix & MultiTypeBlockVector types
+ */
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
+ TVector xold(x);
+ xold=x; //store old x values
+ MultiTypeBlockMatrix_Solver<I,crow,remain_row>::dbgs(A,x,x,b,w);
+ x *= w;
+ x.axpy(1-w,xold); //improve x
+ }
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
+ typename mpl::at_c<TVector,crow>::type rhs;
+ rhs = fusion::at_c<crow> (b);
+
+ 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
+ //solve on blocklevel I-1
+ algmeta_itsteps<I-1>::dbgs(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(x),rhs,w);
+ MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::dbgs(A,x,v,b,w); //next row
+ }
+
+
+
+ /**
+ * bsorf for MultiTypeBlockMatrix & MultiTypeBlockVector types
+ */
+ template <typename TVector, typename TMatrix, typename K>
+ static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
+ TVector v;
+ v=x; //use latest x values in right side calculation
+ MultiTypeBlockMatrix_Solver<I,crow,remain_row>::bsorf(A,x,v,b,w);
+
+ }
+ template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
+ static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
+ typename mpl::at_c<TVector,crow>::type rhs;
+ rhs = fusion::at_c<crow> (b);
+
+ 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
+ //solve on blocklevel I-1
+ algmeta_itsteps<I-1>::bsorf(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
+ fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
+ MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::bsorf(A,x,v,b,w); //next row
+ }
+
+ /**
+ * bsorb for MultiTypeBlockMatrix & MultiTypeBlockVector types
+ */
+ template <typename TVector, typename TMatrix, typename K>
+ static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
+ TVector v;
+ v=x; //use latest x values in right side calculation
+ MultiTypeBlockMatrix_Solver<I,crow,remain_row>::bsorb(A,x,v,b,w);
+
+ }
+ template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
+ static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
+ typename mpl::at_c<TVector,crow>::type rhs;
+ rhs = fusion::at_c<crow> (b);
+
+ 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
+ //solve on blocklevel I-1
+ algmeta_itsteps<I-1>::bsorb(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
+ fusion::at_c<crow>(x).axpy(w,fusion::at_c<crow>(v));
+ MultiTypeBlockMatrix_Solver<I,crow-1,remain_row-1>::bsorb(A,x,v,b,w); //next row
+ }
+
+
+ /**
+ * Jacobi solver for MultiTypeBlockMatrix & MultiTypeBlockVector types
+ */
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
+ TVector v(x);
+ v=0; //calc new x in v
+ MultiTypeBlockMatrix_Solver<I,crow,remain_row>::dbjac(A,x,v,b,w);
+ x.axpy(w,v); //improve x
+ }
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
+ typename mpl::at_c<TVector,crow>::type rhs;
+ rhs = fusion::at_c<crow> (b);
+
+ 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
+ //solve on blocklevel I-1
+ algmeta_itsteps<I-1>::dbjac(fusion::at_c<crow>( fusion::at_c<crow>(A)), fusion::at_c<crow>(v),rhs,w);
+ MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::dbjac(A,x,v,b,w); //next row
+ }
+
+
+
+
+ };
+ template<int I, int crow> //recursion end for remain_row = 0
+ class MultiTypeBlockMatrix_Solver<I,crow,0> {
+ public:
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbgs(const TMatrix& A, TVector& x, TVector& v,
+ const TVector& b, const K& w) {}
+
+ template <typename TVector, typename TMatrix, typename K>
+ static void bsorf(const TMatrix& A, TVector& x, TVector& v,
+ const TVector& b, const K& w) {}
+
+ template <typename TVector, typename TMatrix, typename K>
+ static void bsorb(const TMatrix& A, TVector& x, TVector& v,
+ const TVector& b, const K& w) {}
+
+ template <typename TVector, typename TMatrix, typename K>
+ static void dbjac(const TMatrix& A, TVector& x, TVector& v,
+ const TVector& b, const K& w) {}
+ };
+
+
+} // end namespace
+
+#endif
+
+#endif
diff --git a/dune/istl/multitypeblockvector.hh b/dune/istl/multitypeblockvector.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2481dec36335c9988f2244b723e68ff7918b6b9a
--- /dev/null
+++ b/dune/istl/multitypeblockvector.hh
@@ -0,0 +1,306 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTITYPEVECTOR_HH
+#define DUNE_MULTITYPEVECTOR_HH
+
+#ifdef HAVE_BOOST_FUSION
+
+#include <cmath>
+#include <iostream>
+
+#include "istlexception.hh"
+#include "allocator.hh"
+
+#include "gsetc.hh"
+
+#include <boost/fusion/sequence.hpp>
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/iterator.hpp>
+#include <boost/typeof/typeof.hpp>
+#include <boost/fusion/algorithm.hpp>
+
+
+
+namespace mpl=boost::mpl;
+namespace fusion=boost::fusion;
+
+namespace Dune {
+
+ /**
+ @addtogroup ISTL_SPMV
+ @{
+ */
+
+
+
+
+ /**
+ @brief prints out a vector (type MultiTypeBlockVector)
+
+ The parameter "current_element" is the index of
+ the element to be printed out. Via recursive calling
+ all other elements (number is "remaining_elements")
+ are printed, too. This is internally called when
+ a MultiTypeBlockVector object is passed to an output stream.
+ Example:
+ \code
+ MultiTypeBlockVector<int,int> CVect;
+ std::cout << CVect;
+ \endcode
+ */
+
+ template<int current_element, int remaining_elements, typename TVec>
+ class MultiTypeBlockVector_Print {
+ public:
+ /**
+ * print out the current vector element and all following
+ */
+ static void print(const TVec& v) {
+ std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
+ MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v); //next element
+ }
+ };
+ template<int current_element, typename TVec> //recursion end (remaining_elements=0)
+ class MultiTypeBlockVector_Print<current_element,0,TVec> {
+ public:
+ static void print(const TVec& v) {std::cout << "\n";}
+ };
+
+
+
+ /**
+ @brief set a MultiTypeBlockVector to some specific value
+
+ This class is used by the MultiTypeBlockVector class' internal = operator.
+ Whenever a vector is assigned to a value, each vector element
+ has to provide the = operator for the right side. Example:
+ \code
+ MultiTypeBlockVector<int,int,int> CVect;
+ CVect = 3; //sets all integer elements to 3
+ \endcode
+ */
+ template<int count, typename T1, typename T2>
+ class MultiTypeBlockVector_Ident {
+ public:
+
+ /**
+ * equalize two vectors' element (index is template parameter count)
+ * note: each MultiTypeBlockVector element has to provide the = operator with type T2
+ */
+ static void equalize(T1& a, const T2& b) {
+ fusion::at_c<count-1>(a) = b; //equalize current elements
+ MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b); //next elements
+ }
+ };
+ template<typename T1, typename T2> //recursion end (count=0)
+ class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
+
+
+
+
+
+
+ /**
+ @brief add/sub second vector to/from the first (v1 += v2)
+
+ This class implements vector addition/subtraction for any MultiTypeBlockVector-Class type.
+ */
+ template<int count, typename T>
+ class MultiTypeBlockVector_Add {
+ public:
+
+ /**
+ * add vector to vector
+ */
+ static void add (T& a, const T& b) { //add vector elements
+ fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
+ MultiTypeBlockVector_Add<count-1,T>::add(a,b);
+ }
+
+ /**
+ * sub vector from vector
+ */
+ static void sub (T& a, const T& b) { //sub vector elements
+ fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
+ MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
+ }
+ };
+ template<typename T> //recursion end; specialization for count=0
+ class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
+
+
+
+ /**
+ @brief AXPY operation on vectors
+
+ calculates x += a * y
+ */
+ template<int count, typename TVec, typename Ta>
+ class MultiTypeBlockVector_AXPY {
+ public:
+
+ /**
+ * calculates x += a * y
+ */
+ static void axpy(TVec& x, const Ta& a, const TVec& y) {
+ fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
+ MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
+ }
+ };
+ template<typename TVec, typename Ta> //specialization for count=0
+ class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
+
+
+ /**
+ @brief Scalar * Vector Multiplication
+
+ calculates v *= a for each element of the given vector
+ */
+ template<int count, typename TVec, typename Ta>
+ class MultiTypeBlockVector_Mulscal {
+ public:
+
+ /**
+ * calculates x *= a
+ */
+ static void mul(TVec& x, const Ta& a) {
+ fusion::at_c<(count-1)>(x) *= a;
+ MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
+ }
+ };
+ template<typename TVec, typename Ta> //specialization for count=0
+ class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
+
+
+
+ /**
+ @brief Vector scalar multiplication
+
+ multiplies the the current elements of x and y and recursively
+ and sums it all up.
+ */
+ template<int count, typename TVec>
+ class MultiTypeBlockVector_Mul {
+ public:
+ static typename TVec::field_type mul(const TVec& x, const TVec& y) {
+ return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y);
+ }
+ };
+ template<typename TVec>
+ class MultiTypeBlockVector_Mul<0,TVec> {
+ public: static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
+ };
+
+
+
+
+
+ /**
+ @brief calulate the 2-norm out of vector elements
+
+ Each element of the vector has to provide the method "two_norm2()"
+ in order to calculate the whole vector's 2-norm.
+ */
+ template<int count, typename T>
+ class MultiTypeBlockVector_Norm {
+ public:
+
+ /**
+ * sum up all elements' 2-norms
+ */
+ static double result (const T& a) { //result = sum of all elements' 2-norms
+ return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
+ }
+ };
+ template<typename T> //recursion end: no more vector elements to add...
+ class MultiTypeBlockVector_Norm<0,T> {public: static double result (const T& a) {return 0.0;} };
+
+
+
+
+
+ /**
+ @brief A Vector class to support different block types
+
+ This vector class combines elements of different types known at compile-time.
+ */
+ template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
+ typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_, typename T8=fusion::void_, typename T9=fusion::void_>
+ class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
+
+ public:
+
+ /**
+ * own class' type
+ */
+ typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
+
+ typedef typename T1::field_type field_type;
+
+ /**
+ * number of elements
+ */
+ const int count() {return mpl::size<type>::value;}
+
+ /**
+ * assignment operator
+ */
+ template<typename T>
+ void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
+
+ /**
+ * operator for MultiTypeBlockVector += MultiTypeBlockVector operations
+ */
+ void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
+
+ /**
+ * operator for MultiTypeBlockVector -= MultiTypeBlockVector operations
+ */
+ void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
+
+ void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
+ void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
+ void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
+
+ field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
+
+ /**
+ * two-norm^2
+ */
+ double two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
+
+ /**
+ * the real two-norm
+ */
+ double two_norm() const {return sqrt(this->two_norm2());}
+
+ /**
+ * axpy operation on this vector (*this += a * y)
+ */
+ template<typename Ta>
+ void axpy (const Ta& a, const type& y) {
+ MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
+ }
+
+ };
+
+
+
+ /**
+ @brief << operator for a MultiTypeBlockVector
+
+ operator<< for printing out a MultiTypeBlockVector
+ */
+ template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
+ std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
+ MultiTypeBlockVector_Print<0,mpl::size<MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value,MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(v);
+ return s;
+ }
+
+
+
+} // end namespace
+
+#endif
+
+#endif
diff --git a/dune/istl/tutorial/Makefile.am b/dune/istl/tutorial/Makefile.am
index 2cb421440402f4706fa26047e98b7ac9203bc896..a7ecf900217cb36b7d0d83090a434a4a19d81de4 100644
--- a/dune/istl/tutorial/Makefile.am
+++ b/dune/istl/tutorial/Makefile.am
@@ -4,5 +4,6 @@ dist_noinst_DATA = example.cc
noinst_PROGRAMS = example
example_SOURCES = example.cc
+example_CPPFLAGS = $(BOOST_CPPFALGS)
include $(top_srcdir)/am/global-rules
diff --git a/dune/istl/tutorial/example.cc b/dune/istl/tutorial/example.cc
index 217c85dd6c93d61ff25500b29f2f53e4266ec429..60ad1a3cb75912965c8452e35607649fde5c19c5 100644
--- a/dune/istl/tutorial/example.cc
+++ b/dune/istl/tutorial/example.cc
@@ -26,6 +26,10 @@
#include <dune/istl/preconditioners.hh>
#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+
+
// a simple stop watch
class Timer
{
@@ -503,6 +507,100 @@ void test_Interface ()
}
+#ifdef HAVE_BOOST_FUSION
+
+void test_MultiTypeBlockVector_MultiTypeBlockMatrix() { //Jacobi Solver Test MultiTypeBlockMatrix_Solver::dbjac on MultiTypeBlockMatrix<BCRSMatrix>
+
+ std::cout << "\n\n\nJacobi Solver Test on MultiTypeBlockMatrix<BCRSMatrix>\n";
+
+ typedef Dune::FieldMatrix<double,1,1> LittleBlock; //matrix block type
+ typedef Dune::BCRSMatrix<LittleBlock> BCRSMat; //matrix type
+ typedef Dune::BlockVector<Dune::FieldVector<double,1> > VecType;
+
+ const int X1=3; //index bounds of all four matrices
+ const int X2=2;
+ const int Y1=3;
+ const int Y2=2;
+ BCRSMat A11 = BCRSMat(X1,Y1,X1*Y1,BCRSMat::random); //A11 is 3x3
+ BCRSMat A12 = BCRSMat(X1,Y2,X1*Y2,BCRSMat::random); //A12 is 2x3
+ BCRSMat A21 = BCRSMat(X2,Y1,X2*Y1,BCRSMat::random); //A11 is 3x2
+ BCRSMat A22 = BCRSMat(X2,Y2,X2*Y2,BCRSMat::random); //A12 is 2x2
+
+ typedef Dune::MultiTypeBlockVector<Dune::BlockVector<Dune::FieldVector<double,1> >,Dune::BlockVector<Dune::FieldVector<double,1> > > TestVector;
+ TestVector x, b;
+
+ fusion::at_c<0>(x).resize(Y1);
+ fusion::at_c<1>(x).resize(Y2);
+ fusion::at_c<0>(b).resize(X1);
+ fusion::at_c<1>(b).resize(X2);
+
+ x = 1; b = 1;
+
+ //set row sizes
+ for (int i=0; i<Y1; i++) {A11.setrowsize(i,X1); A12.setrowsize(i,X2);} A11.endrowsizes(); A12.endrowsizes();
+ for (int i=0; i<Y2; i++) {A21.setrowsize(i,X1); A22.setrowsize(i,X2);} A21.endrowsizes(); A22.endrowsizes();
+
+ //set indices
+ for (int i=0; i<X1+X2; i++) {
+ for (int j=0; j<Y1+Y2; j++) {
+ if (i<X1 && j<Y1) {A11.addindex(i,j);}
+ if (i<X1 && j>=Y1) {A12.addindex(i,j-Y1);}
+ if (i>=X1 && j<Y1) {A21.addindex(i-X1,j);}
+ if (i>=X1 && j>=Y1) {A22.addindex(i-X1,j-Y1);}
+ }
+ }
+ A11.endindices(); A12.endindices(); A21.endindices(); A22.endindices();
+ A11 = 0; A12 = 0; A21 = 0; A22 = 0;
+
+ //fill in values (row-wise) in A11 and A22
+ for (int i=0; i<Y1; i++) {
+ if (i>0) {A11[i][i-1]=-1;}
+ A11[i][i]=2; //diag
+ if (i<Y1-1) {A11[i][i+1]=-1;}
+ if (i<Y2) { //also in A22
+ if (i>0) {A22[i][i-1]=-1;}
+ A22[i][i]=2;
+ if (i<Y2-1) {A22[i][i+1]=-1;}
+ }
+ }
+ A12[2][0] = -1; A21[0][2] = -1; //fill A12 & A21
+
+
+ typedef Dune::MultiTypeBlockVector<BCRSMat,BCRSMat> BCRS_Row;
+ typedef Dune::MultiTypeBlockMatrix<BCRS_Row,BCRS_Row> CM_BCRS;
+ CM_BCRS A;
+ fusion::at_c<0>(fusion::at_c<0>(A)) = A11;
+ fusion::at_c<1>(fusion::at_c<0>(A)) = A12;
+ fusion::at_c<0>(fusion::at_c<1>(A)) = A21;
+ fusion::at_c<1>(fusion::at_c<1>(A)) = A22;
+
+ printmatrix(std::cout,A11,"matrix A11","row",9,1);
+ printmatrix(std::cout,A12,"matrix A12","row",9,1);
+ printmatrix(std::cout,A21,"matrix A21","row",9,1);
+ printmatrix(std::cout,A22,"matrix A22","row",9,1);
+
+ x = 1;
+ b = 1;
+
+ Dune::MatrixAdapter<CM_BCRS,TestVector,TestVector> op(A); // make linear operator from A
+ Dune::SeqJac<CM_BCRS,TestVector,TestVector,2> jac(A,1,1); // Jacobi preconditioner
+ Dune::SeqGS<CM_BCRS,TestVector,TestVector,2> gs(A,1,1); // GS preconditioner
+ Dune::SeqSOR<CM_BCRS,TestVector,TestVector,2> sor(A,1,1.9520932); // SSOR preconditioner
+ Dune::SeqSSOR<CM_BCRS,TestVector,TestVector,2> ssor(A,1,1.0); // SSOR preconditioner
+
+ Dune::LoopSolver<TestVector> loop(op,gs,1E-4,18000,2); // an inverse operator
+ Dune::InverseOperatorResult r;
+
+ loop.apply(x,b,r);
+
+ printvector(std::cout,fusion::at_c<0>(x),"solution x1","entry",11,9,1);
+ printvector(std::cout,fusion::at_c<1>(x),"solution x2","entry",11,9,1);
+
+}
+#endif
+
+
+
int main (int argc , char ** argv)
{
@@ -515,6 +613,9 @@ int main (int argc , char ** argv)
test_IO();
test_Iter();
test_Interface();
+#ifdef HAVE_BOOST_FUSION
+ test_MultiTypeBlockVector_MultiTypeBlockMatrix();
+#endif
}
catch (Dune::ISTLError& error)
{
diff --git a/m4/dune_istl.m4 b/m4/dune_istl.m4
index 3789841df2e8999b606f6a123bf5a1d8376208f9..191523531c924af540627b5616509ef6c265471f 100644
--- a/m4/dune_istl.m4
+++ b/m4/dune_istl.m4
@@ -10,6 +10,8 @@ AC_DEFUN([DUNE_ISTL_CHECKS],
AC_REQUIRE([__AC_FC_NAME_MANGLING])
AC_REQUIRE([AC_PROG_F77])
AC_REQUIRE([ACX_BLAS])
+ AX_BOOST_BASE(1.35, [ DUNE_BOOST_FUSION ] , [] )
+
# add summary entries for tests not maintained by dune
DUNE_ADD_SUMMARY_ENTRY([METIS],[$with_metis])
DUNE_ADD_SUMMARY_ENTRY([BLAS],[$acx_blas_ok])