Added a dense block matrix where each block can have a different type.

The code was developed by Theo Aouidet, Stefan Bayha, and Dominik Ull
during their software internship at the IPVS in Stuttgart.

[[Imported from SVN: r1259]]

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 \

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)
     {

dune/istl/multitypeblockmatrix.hh

+// -*- 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

dune/istl/multitypeblockvector.hh

+// -*- 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

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

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)
   {

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])