libdune-istl-dev 2.2.1-2 / usr / include / dune / istl / overlappingschwarz.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_OVERLAPPINGSCHWARZ_HH
#define DUNE_OVERLAPPINGSCHWARZ_HH
#include<cassert>
#include<algorithm>
#include<functional>
#include<vector>
#include<set>
#include<dune/common/dynmatrix.hh>
#include<dune/common/sllist.hh>
#include"preconditioners.hh"
#include"superlu.hh"
#include"bvector.hh"
#include"bcrsmatrix.hh"
#include"ilusubdomainsolver.hh"

namespace Dune
{

  /**
   * @addtogroup ISTL_Prec
   *
   * @{
   */
  /**
   * @file
   * @author Markus Blatt
   * @brief Contains one level overlapping Schwarz preconditioners
   */

  template<class M, class X, class TM, class TD, class TA>
  class SeqOverlappingSchwarz;
  
  /**
   * @brief Initializer for SuperLU Matrices representing the subdomains.
   */
  template<class I, class S, class D>
  class OverlappingSchwarzInitializer
  {
  public:
    /** @brief The vector type containing the subdomain to row index mapping. */
    typedef D subdomain_vector;

    typedef I InitializerList;
    typedef typename InitializerList::value_type AtomInitializer;
    typedef typename AtomInitializer::Matrix Matrix;
    typedef typename Matrix::const_iterator Iter;
    typedef typename Matrix::row_type::const_iterator CIter;
    
    typedef S IndexSet;
    typedef typename IndexSet::size_type size_type;
    
    OverlappingSchwarzInitializer(InitializerList& il,
                                  const IndexSet& indices,
                                  const subdomain_vector& domains);
    
    
    void addRowNnz(const Iter& row);
    
    void allocate();
        
    void countEntries(const Iter& row, const CIter& col) const;

    void calcColstart() const;
    
    void copyValue(const Iter& row, const CIter& col) const;
    
    void createMatrix() const;
  private:
    class IndexMap
    {
    public:
      typedef typename S::size_type size_type;
      typedef std::map<size_type,size_type> Map;      
      typedef typename Map::iterator iterator;
      typedef typename Map::const_iterator const_iterator;
      
      IndexMap();
      
      void insert(size_type grow);
      
      const_iterator find(size_type grow) const;
      
      iterator find(size_type grow);

      iterator begin();
      
      const_iterator begin()const;

      iterator end();
      
      const_iterator end() const;
            
    private:
      std::map<size_type,size_type> map_;
      size_type row;
    };
    
    
    typedef typename InitializerList::iterator InitIterator;
    typedef typename IndexSet::const_iterator IndexIteratur;
    InitializerList* initializers;
    const IndexSet *indices;
    mutable std::vector<IndexMap> indexMaps;
    const subdomain_vector& domains;
  };

  /** 
   * @brief Tag that the tells the schwarz method to be additive.
   */
  struct AdditiveSchwarzMode
  {};

  /**
   * @brief Tag that tells the Schwarz method to be multiplicative.
   */
  struct MultiplicativeSchwarzMode
  {};
  
  /**
   * @brief Tag that tells the Schwarz method to be multiplicative
   * and symmetric.
   */
  struct SymmetricMultiplicativeSchwarzMode
  {};  

  /**
   * @brief Exact subdomain solver using Dune::DynamicMatrix<T>::solve   
   * @tparam M The type of the matrix.
   */
  template<class M, class X, class Y>
  class DynamicMatrixSubdomainSolver;

  // Specialization for BCRSMatrix
  template<class K, int n, class Al, class X, class Y>
  class DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >
  {
    typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> M;
  public:
    //! \brief The matrix type the preconditioner is for.
    typedef typename Dune::remove_const<M>::type matrix_type;
    typedef K field_type;
    typedef typename Dune::remove_const<M>::type rilu_type;
    //! \brief The domain type of the preconditioner.
    typedef X domain_type;
    //! \brief The range type of the preconditioner.
    typedef Y range_type;
    
    /** 
     * @brief Apply the subdomain solver.
     * @copydoc ILUSubdomainSolver::apply
     */
    void apply (DynamicVector<field_type>& v, DynamicVector<field_type>& d)
    {
      assert(v.size() > 0);
      assert(v.size() == d.size());
      assert(A.rows() <= v.size());
      assert(A.cols() <= v.size());
      size_t sz = A.rows();
      v.resize(sz);
      d.resize(sz);
      A.solve(v,d);
      v.resize(v.capacity());
      d.resize(d.capacity());
    }

    /**
     * @brief Set the data of the local problem.
     * 
     * @param BCRS The global matrix.
     * @param rowset The global indices of the local problem.
     * @tparam S The type of the set with the indices.
     */
    template<class S>
    void setSubMatrix(const M& BCRS, S& rowset)
    {
      size_t sz = rowset.size();
      A.resize(sz*n,sz*n);
      typename DynamicMatrix<K>::RowIterator rIt = A.begin();
      typedef typename S::const_iterator SIter;
      size_t r = 0, c = 0;
      for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
          rowIdx!= rowEnd; ++rowIdx, r++)
      {
        c = 0;
        for(SIter colIdx = rowset.begin(), colEnd=rowset.end(); 
            colIdx!= colEnd; ++colIdx, c++)
        {
          if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
            continue;
          for (size_t i=0; i<n; i++)
          {
            for (size_t j=0; j<n; j++)
            {
              A[r*n+i][c*n+j] = BCRS[*rowIdx][*colIdx][i][j];
            }
          }
        }
      }
    }
  private:
    DynamicMatrix<K> A;
  };

  template<typename T>
  struct OverlappingAssigner
  {
  };
    
  // specialization for DynamicMatrix
  template<class K, int n, class Al, class X, class Y>
  class OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  {
  public:
    typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
    typedef K field_type;
    typedef Y range_type;
    typedef typename range_type::block_type block_type;
    typedef typename matrix_type::size_type size_type;
    
    /**
     * @brief Constructor.
     * @param maxlength The maximum entries over all subdomains.
     * @param mat_ The global matrix.
     * @param b_ the global right hand side.
     * @param x_ the global left hand side.
     */
    OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, const X& b_, Y& x_);

    /**
     * @brief Deallocates memory of the local vector.
     */
    inline
    void deallocate();
    
    /**
     * @brief Resets the local index to zero.
     */
    inline
    void resetIndexForNextDomain();
    
    /**
     * @brief Get the local left hand side.
     * @return The local left hand side.
     */
    inline
    DynamicVector<K> & lhs();
    
    /**
     * @brief Get the local right hand side.
     * @return The local right hand side.
     */
    inline
    DynamicVector<K> & rhs();
    
    /**
     * @brief relax the result.
     * @param relax The relaxation parameter.
     */
    inline
    void relaxResult(field_type relax);
    
    /**
     * @brief calculate one entry of the local defect.
     * @param domainIndex One index of the domain.
     */
    void operator()(const size_type& domainIndex);
    
    /**
     * @brief Assigns the block to the current local
     * index.
     * At the same time the local defect is calculated 
     * for the index and stored in the rhs.
     * Afterwards the is incremented for the next block.
     */
    inline
    void assignResult(block_type& res);
    
  private:
    /**
     * @brief The global matrix for the defect calculation.
     */
    const matrix_type* mat;
    /** @brief The local right hand side. */
    // we need a pointer, because we have to avoid deep copies
    DynamicVector<field_type> * rhs_;
    /** @brief The local left hand side. */
    // we need a pointer, because we have to avoid deep copies
    DynamicVector<field_type> * lhs_;
    /** @brief The global right hand side for the defect calculation. */
    const range_type* b;
    /** @brief The global right hand side for adding the local update to. */
    range_type* x;
    /** @brief The current local index. */
    std::size_t i;
    /** @brief The maximum entries over all subdomains. */
    std::size_t maxlength_;
  };

#if HAVE_SUPERLU
  template<typename T, typename A, int n, int m>
  struct OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
    {
      typedef BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
      typedef typename SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::range_type range_type;
      typedef typename range_type::field_type field_type;
      typedef typename range_type::block_type block_type;
      
      typedef typename matrix_type::size_type size_type;
      
      /**
       * @brief Constructor.
       * @param maxlength The maximum entries over all subdomains.
       * @param mat The global matrix.
       * @param b the global right hand side.
       * @param x the global left hand side.
       */
      OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat, 
                          const range_type& b, range_type& x);
      /**
       * @brief Deallocates memory of the local vector.
       * @warning memory is released by the destructor as this Functor
       * is copied and the copy needs to still have the data.
       */
      void deallocate();
      
      /*
       * @brief Resets the local index to zero.
       */
      void resetIndexForNextDomain();
      
      /**
       * @brief Get the local left hand side.
       * @return The local left hand side.
       */
      field_type* lhs();
            
      /**
       * @brief Get the local right hand side.
       * @return The local right hand side.
       */
      field_type* rhs();
      
      /**
       * @brief relax the result.
       * @param relax The relaxation parameter.
       */
      void relaxResult(field_type relax);

      /**
       * @brief calculate one entry of the local defect.
       * @param domain One index of the domain.
       */
      void operator()(const size_type& domain);

      /**
       * @brief Assigns the block to the current local
       * index.
       * At the same time the local defect is calculated 
       * for the index and stored in the rhs.
       * Afterwards the is incremented for the next block.
       */
      void assignResult(block_type& res);
      
    private:
      /**
       * @brief The global matrix for the defect calculation.
       */
      const matrix_type* mat;
      /** @brief The local right hand side. */
      field_type* rhs_;
      /** @brief The local left hand side. */
      field_type* lhs_;
      /** @brief The global right hand side for the defect calculation. */
      const range_type* b;
      /** @brief The global right hand side for adding the local update to. */
      range_type* x;
      /** @brief The current local index. */
      std::size_t i;
      /** @brief The maximum entries over all subdomains. */
      std::size_t maxlength_;
    };

#endif

  template<class M, class X, class Y>
  class OverlappingAssignerILUBase
  {
  public:
    typedef M matrix_type;
    
    typedef typename M::field_type field_type;
  
    typedef typename Y::block_type block_type;
      
    typedef typename matrix_type::size_type size_type;
    /**
     * @brief Constructor.
     * @param maxlength The maximum entries over all subdomains.
     * @param mat The global matrix.
     * @param b the global right hand side.
     * @param x the global left hand side.
     */
    OverlappingAssignerILUBase(std::size_t maxlength, const M& mat, 
                        const Y& b, X& x);
    /**
     * @brief Deallocates memory of the local vector.
     * @warning memory is released by the destructor as this Functor
     * is copied and the copy needs to still have the data.
     */
    void deallocate();
    
    /**
     * @brief Resets the local index to zero.
     */
    void resetIndexForNextDomain();
    
    /**
     * @brief Get the local left hand side.
     * @return The local left hand side.
     */
    X& lhs();
    
    /**
     * @brief Get the local right hand side.
     * @return The local right hand side.
     */
    Y& rhs();
    
    /**
     * @brief relax the result.
     * @param relax The relaxation parameter.
     */
    void relaxResult(field_type relax);
    
    /**
     * @brief calculate one entry of the local defect.
     * @param domain One index of the domain.
     */
    void operator()(const size_type& domain);
    
    /**
       * @brief Assigns the block to the current local
       * index.
       * At the same time the local defect is calculated 
       * for the index and stored in the rhs.
       * Afterwards the is incremented for the next block.
       */
    void assignResult(block_type& res);

  private:
    /**
     * @brief The global matrix for the defect calculation.
     */
    const M* mat;
    /** @brief The local left hand side. */
    X* lhs_;
    /** @brief The local right hand side. */
    Y* rhs_;
    /** @brief The global right hand side for the defect calculation. */
    const Y* b;
    /** @brief The global left hand side for adding the local update to. */
    X* x;
    /** @brief The maximum entries over all subdomains. */
    size_type i;
  };
  
  // specialization for ILU0
  template<class M, class X, class Y>
  class OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >
    : public OverlappingAssignerILUBase<M,X,Y>
  {
  public:
    /**
     * @brief Constructor.
     * @param maxlength The maximum entries over all subdomains.
     * @param mat The global matrix.
     * @param b the global right hand side.
     * @param x the global left hand side.
     */
    OverlappingAssigner(std::size_t maxlength, const M& mat, 
                        const Y& b, X& x)
      : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
    {}
  };

  // specialization for ILUN
  template<class M, class X, class Y>
  class OverlappingAssigner<ILUNSubdomainSolver<M,X,Y> >
    : public OverlappingAssignerILUBase<M,X,Y>
  {
  public:
    /**
     * @brief Constructor.
     * @param maxlength The maximum entries over all subdomains.
     * @param mat The global matrix.
     * @param b the global right hand side.
     * @param x the global left hand side.
     */
    OverlappingAssigner(std::size_t maxlength, const M& mat, 
                        const Y& b, X& x)
      : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
    {}
  };

    template<typename S, typename T>
    struct AdditiveAdder
    {
    };
    
    template<typename S, typename T, typename A, int n>
    struct AdditiveAdder<S, BlockVector<FieldVector<T,n>,A> >
    {
      typedef typename A::size_type size_type;
      AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
                    OverlappingAssigner<S>& assigner, const T& relax_);
      void operator()(const size_type& domain);
      void axpy();
      
    private:
      BlockVector<FieldVector<T,n>,A>* v;
      BlockVector<FieldVector<T,n>,A>* x;
      OverlappingAssigner<S>* assigner;
      T relax;
    };

    template<typename S,typename T>
    struct MultiplicativeAdder
    {
    };
    
    template<typename S, typename T, typename A, int n>
    struct MultiplicativeAdder<S, BlockVector<FieldVector<T,n>,A> >
    {
      typedef typename A::size_type size_type;
      MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x, 
                          OverlappingAssigner<S>& assigner_, const T& relax_);
      void operator()(const size_type& domain);
      void axpy();
      
    private:
      BlockVector<FieldVector<T,n>,A>* x;
      OverlappingAssigner<S>* assigner;
      T relax;
    };

    /**
     * @brief template meta program for choosing  how to add the correction.
     *
     * There are specialization for the additive, the multiplicative, and the symmetric multiplicative mode.
     *
     * \tparam T The Schwarz mode (either AdditiveSchwarzMode or MuliplicativeSchwarzMode or 
     * SymmetricMultiplicativeSchwarzMode)
     * \tparam X The vector field type
     */
    template<typename T, class X, class S>
    struct AdderSelector
    {};
    
  template<class X, class S>
    struct AdderSelector<AdditiveSchwarzMode,X,S>
    {
      typedef AdditiveAdder<S,X> Adder;
    };
    
  template<class X, class S>
    struct AdderSelector<MultiplicativeSchwarzMode,X,S>
    {
      typedef MultiplicativeAdder<S,X> Adder;
    };    

  template<class X, class S>
    struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
    {
      typedef MultiplicativeAdder<S,X> Adder;
    };    

    /**
     * @brief Helper template meta program for application of overlapping schwarz.
     *
     * The is needed because when using the multiplicative schwarz version one
     * might still want to make multigrid symmetric, i.e. forward sweep when pre- 
     * and backward sweep when post-smoothing.
     *
     * @tparam T1 type of the vector with the subdomain solvers.
     * @tparam T2 type of the vector with the subdomain vector fields.
     * @tparam forward If true apply in a forward sweep.
     */
    template<typename T1, typename T2, bool forward>
    struct IteratorDirectionSelector
    {
      typedef T1 solver_vector;
      typedef typename solver_vector::iterator solver_iterator;
      typedef T2 subdomain_vector;
      typedef typename subdomain_vector::const_iterator domain_iterator;
    
      static solver_iterator begin(solver_vector& sv)
      {
        return sv.begin();
      }
      
      static solver_iterator end(solver_vector& sv)
      {
        return sv.end();
      }
      static domain_iterator begin(const subdomain_vector& sv)
      {
        return sv.begin();
      }
      
      static domain_iterator end(const subdomain_vector& sv)
      {
        return sv.end();
      }
    };

    template<typename T1, typename T2>
    struct IteratorDirectionSelector<T1,T2,false>
    {
      typedef T1 solver_vector;
      typedef typename solver_vector::reverse_iterator solver_iterator;
      typedef T2 subdomain_vector;
      typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
    
      static solver_iterator begin(solver_vector& sv)
      {
        return sv.rbegin();
      }
      
      static solver_iterator end(solver_vector& sv)
      {
        return sv.rend();
      }
      static domain_iterator begin(const subdomain_vector& sv)
      {
        return sv.rbegin();
      }
      
      static domain_iterator end(const subdomain_vector& sv)
      {
        return sv.rend();
      }
    };
    
    /**
     * @brief Helper template meta program for application of overlapping schwarz.
     *
     * The is needed because when using the multiplicative schwarz version one
     * might still want to make multigrid symmetric, i.e. forward sweep when pre- 
     * and backward sweep when post-smoothing.
     * @tparam T The smoother to apply.
     */
    template<class T>
    struct SeqOverlappingSchwarzApplier
    {
      typedef T smoother;
      typedef typename smoother::range_type range_type;
      
      static void apply(smoother& sm, range_type& v, const range_type& b)
      {
        sm.template apply<true>(v, b);
      }
    };
    
    template<class M, class X, class TD, class TA>
    struct SeqOverlappingSchwarzApplier<SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode,TD,TA> >
    {
      typedef SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode,TD,TA> smoother;
      typedef typename smoother::range_type range_type;
      
      static void apply(smoother& sm, range_type& v, const range_type& b)
      {
        sm.template apply<true>(v, b);
        sm.template apply<false>(v, b);
      }
    };

  template<class T>
  struct SeqOverlappingSchwarzAssembler
  {};
  

  template<class K, int n, class Al, class X, class Y>
  struct SeqOverlappingSchwarzAssembler< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  {
    typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
    template<class RowToDomain, class Solvers, class SubDomains>
    static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
                                             Solvers& solvers, const SubDomains& domains,
                                             bool onTheFly);
  };

#if HAVE_SUPERLU
  template<class T>
  struct SeqOverlappingSchwarzAssembler<SuperLU<T> >
  {
    typedef T matrix_type;
    template<class RowToDomain, class Solvers, class SubDomains>
    static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
                                             Solvers& solvers, const SubDomains& domains,
                                             bool onTheFly);
  };
#endif

  template<class M,class X, class Y>
  struct SeqOverlappingSchwarzAssemblerILUBase
  {
    typedef M matrix_type;
    template<class RowToDomain, class Solvers, class SubDomains>
    static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
                                             Solvers& solvers, const SubDomains& domains,
                                             bool onTheFly);
  };

  template<class M,class X, class Y>
  struct SeqOverlappingSchwarzAssembler<ILU0SubdomainSolver<M,X,Y> >
    : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
  {};

  template<class M,class X, class Y>
  struct SeqOverlappingSchwarzAssembler<ILUNSubdomainSolver<M,X,Y> >
    : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
  {};
  
  /**
   * @brief Sequential overlapping Schwarz preconditioner
   *
   * @tparam M The matrix type.
   * @tparam X The range and domain type.
   * @tparam TM The Schwarz mode. Currently supported modes are AdditiveSchwarzMode,
   * MultiplicativeSchwarzMode, and SymmetricMultiplicativeSchwarzMode. (Default values is AdditiveSchwarzMode)
   * @tparam TD The type of the local subdomain solver to be used.
   * @tparam TA The type of the allocator to use.
   */
  template<class M, class X, class TM=AdditiveSchwarzMode, 
           class TD=ILU0SubdomainSolver<M,X,X>, class TA=std::allocator<X> >
  class SeqOverlappingSchwarz
    : public Preconditioner<X,X>
  {
  public:
    /**
     * @brief The type of the matrix to precondition.
     */
    typedef M matrix_type;

    /**
     * @brief The domain type of the preconditioner
     */
    typedef X domain_type;

    /**
     * @brief The range type of the preconditioner.
     */
    typedef X range_type;

    /**
     * @brief The mode (additive or multiplicative) of the Schwarz
     * method.
     *
     * Either AdditiveSchwarzMode or MultiplicativeSchwarzMode
     */
    typedef TM Mode;

    /**
     * @brief The field type of the preconditioner.
     */
    typedef typename X::field_type field_type;

    /** @brief The return type of the size method. */
    typedef typename matrix_type::size_type size_type;

    /** @brief The allocator to use. */
    typedef TA allocator;
        
    /** @brief The type for the subdomain to row index mapping. */
    typedef std::set<size_type, std::less<size_type>, 
		     typename TA::template rebind<size_type>::other> 
      subdomain_type;

    /** @brief The vector type containing the subdomain to row index mapping. */
    typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other> subdomain_vector;

    /** @brief The type for the row to subdomain mapping. */
    typedef SLList<size_type, typename TA::template rebind<size_type>::other> subdomain_list;
      
    /** @brief The vector type containing the row index to subdomain mapping. */
    typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other > rowtodomain_vector;

    /** @brief The type for the subdomain solver in use. */
    typedef TD slu;
    
    /** @brief The vector type containing subdomain solvers. */
    typedef std::vector<slu, typename TA::template rebind<slu>::other> slu_vector;

    enum{
      //! \brief The category the precondtioner is part of.
      category = SolverCategory::sequential
        };
    
    /**
     * @brief Construct the overlapping Schwarz method.
     * @param mat The matrix to precondition.
     * @param subDomains Array of sets of rowindices belonging to an overlapping
     * subdomain
     * @param relaxationFactor relaxation factor
     * @param onTheFly_ If true the decomposition of the exact local solvers is 
     * computed on the fly for each subdomain and 
     * iteration step. If false all decompositions are computed in pre and 
     * only forward and backward substitution takes place 
     * in the iteration steps.
     * @warning Each rowindex should be part of at least one subdomain!
     */
    SeqOverlappingSchwarz(const matrix_type& mat, const subdomain_vector& subDomains,
                          field_type relaxationFactor=1, bool onTheFly_=true);

    /**
     * Construct the overlapping Schwarz method
     * @param mat The matrix to precondition.
     * @param rowToDomain The mapping of the rows onto the domains.
     * @param relaxationFactor relaxation factor
     * @param onTheFly_ If true the decomposition of the exact local solvers is 
     * computed on the fly for each subdomain and 
     * iteration step. If false all decompositions are computed in pre and 
     * only forward and backward substitution takes place 
     * in the iteration steps.
     */
    SeqOverlappingSchwarz(const matrix_type& mat, const rowtodomain_vector& rowToDomain,
                          field_type relaxationFactor=1, bool onTheFly_=true);
                          
    /*! 
      \brief Prepare the preconditioner.
      
      \copydoc Preconditioner::pre(X&,Y&)
    */
    virtual void pre (X& x, X& b) {}

    /*! 
      \brief Apply the precondtioner
      
      \copydoc Preconditioner::apply(X&,const Y&)
    */
    virtual void apply (X& v, const X& d);

    template<bool forward>
    void apply(X& v, const X& d);
    
    /*! 
      \brief Clean up.
      
      \copydoc Preconditioner::post(X&)
    */
    virtual void post (X& x) {
      Dune::dverb<<" avg nnz over subdomain is "<<nnz<<std::endl;
    }
  
  private:
    const M& mat;
    slu_vector solvers;
    subdomain_vector subDomains;
    field_type relax;
    
    typename M::size_type maxlength;
    std::size_t nnz;

    bool onTheFly;
  };
  

    
  template<class I, class S, class D>
  OverlappingSchwarzInitializer<I,S,D>::OverlappingSchwarzInitializer(InitializerList& il,
                                                                    const IndexSet& idx,
                                                                    const subdomain_vector& domains_)
    : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
  {
  }
  
  
  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(const Iter& row)
  {
    typedef typename IndexSet::value_type::const_iterator iterator;
    for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
      (*initializers)[*domain].addRowNnz(row, domains[*domain]);
      indexMaps[*domain].insert(row.index());
    }
  }
  
  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::allocate()
  {
    std::for_each(initializers->begin(), initializers->end(),
                  std::mem_fun_ref(&AtomInitializer::allocateMatrixStorage));
    std::for_each(initializers->begin(), initializers->end(), 
                  std::mem_fun_ref(&AtomInitializer::allocateMarker));
  }
  
  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::countEntries(const Iter& row, const CIter& col) const
  {
    typedef typename IndexSet::value_type::const_iterator iterator;
    for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
      typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
      if(v!= indexMaps[*domain].end()){
        (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
      }
    }
  }

  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::calcColstart() const
  {
    std::for_each(initializers->begin(), initializers->end(),
                  std::mem_fun_ref(&AtomInitializer::calcColstart));
  }

  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::copyValue(const Iter& row, const CIter& col) const
  {
    typedef typename IndexSet::value_type::const_iterator iterator;
    for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain){
      typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
      if(v!= indexMaps[*domain].end()){
        assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
        (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second, 
                                           v->second);
      }
    }
  }
    
  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::createMatrix() const
  {
    indexMaps.clear();
    indexMaps.swap(std::vector<IndexMap>(indexMaps));
    std::for_each(initializers->begin(), initializers->end(),
                  std::mem_fun_ref(&AtomInitializer::createMatrix));
  }

  template<class I, class S, class D>
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
    : row(0)
  {}

  template<class I, class S, class D>
  void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
  {
    assert(map_.find(grow)==map_.end());
    map_.insert(std::make_pair(grow, row++));
  }

  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow) const
  {
    return map_.find(grow);
  }
  
  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
  {
    return map_.find(grow);
  }

  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::end() const
  {
    return map_.end();
  }
  
  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
  {
    return map_.end();
  }

  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin() const
  {
    return map_.begin();
  }
  
  template<class I, class S, class D>
  typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
  OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
  {
    return map_.begin();
  }

  template<class M, class X, class TM, class TD, class TA>
  SeqOverlappingSchwarz<M,X,TM,TD,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const rowtodomain_vector& rowToDomain,
                                                          field_type relaxationFactor, bool fly)
    : mat(mat_), relax(relaxationFactor), onTheFly(fly)
  {
    typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
    typedef typename subdomain_list::const_iterator DomainIterator;
#ifdef DUNE_ISTL_WITH_CHECKING
    assert(rowToDomain.size()==mat.N());
    assert(rowToDomain.size()==mat.M());
    
    for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
      assert(iter->size()>0);
    
#endif
    // calculate the number of domains
    size_type domains=0;
    for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
      for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
        domains=std::max(domains, *d);
    ++domains;
    
    solvers.resize(domains);
    subDomains.resize(domains);
    
    // initialize subdomains to row mapping from row to subdomain mapping
    size_type row=0;
    for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
      for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
        subDomains[*d].insert(row);

#ifdef DUNE_ISTL_WITH_CHECKING
    size_type i=0;
    typedef typename subdomain_vector::const_iterator iterator;
    for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter){
      typedef typename subdomain_type::const_iterator entry_iterator;
      Dune::dvverb<<"domain "<<i++<<":";
      for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry){
        Dune::dvverb<<" "<<*entry;
      }
      Dune::dvverb<<std::endl;
    }
#endif
    maxlength = SeqOverlappingSchwarzAssembler<slu>
      ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
  }
  
  template<class M, class X, class TM, class TD, class TA>
  SeqOverlappingSchwarz<M,X,TM,TD,TA>::SeqOverlappingSchwarz(const matrix_type& mat_,
                                                          const subdomain_vector& sd,
                                                          field_type relaxationFactor,
                                                          bool fly)
    :  mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
       onTheFly(fly)
  {
    typedef typename subdomain_vector::const_iterator DomainIterator;

#ifdef DUNE_ISTL_WITH_CHECKING
    size_type i=0;
    
    for(DomainIterator d=sd.begin(); d != sd.end(); ++d,++i){
      //std::cout<<i<<": "<<d->size()<<std::endl;
      assert(d->size()>0);
      typedef typename DomainIterator::value_type::const_iterator entry_iterator;
      Dune::dvverb<<"domain "<<i<<":";
      for(entry_iterator entry = d->begin(); entry != d->end(); ++entry){
        Dune::dvverb<<" "<<*entry;
      }
      Dune::dvverb<<std::endl;
    }
    
#endif
    
    // Create a row to subdomain mapping
    rowtodomain_vector rowToDomain(mat.N());

    size_type domainId=0;
    
    for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId){
      typedef typename subdomain_type::const_iterator iterator;
      for(iterator row=domain->begin(); row != domain->end(); ++row)
        rowToDomain[*row].push_back(domainId);
    }

    maxlength = SeqOverlappingSchwarzAssembler<slu>
      ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
  }

  /** 
     template helper struct to determine the size of a domain for the
     SeqOverlappingSchwarz solver
     
     only implemented for BCRSMatrix<FieldMatrix<T,n,m>
  */ 
  template<class M>
  struct SeqOverlappingSchwarzDomainSize {};
  
  template<typename T, typename A, int n, int m>
  struct SeqOverlappingSchwarzDomainSize<BCRSMatrix<FieldMatrix<T,n,m>,A > >
  {
    template<class Domain>
    static int size(const Domain & d)
    {
      assert(n==m);
      return m*d.size();
    }
  };
  
  template<class K, int n, class Al, class X, class Y>
  template<class RowToDomain, class Solvers, class SubDomains>
  std::size_t
  SeqOverlappingSchwarzAssembler< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >::
  assembleLocalProblems(const RowToDomain& rowToDomain, 
    const matrix_type& mat,
    Solvers& solvers,
    const SubDomains& subDomains,
    bool onTheFly)
  {
    typedef typename SubDomains::const_iterator DomainIterator;
    std::size_t maxlength = 0;
    
    assert(onTheFly);
    
    for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
      maxlength=std::max(maxlength, domain->size());
    maxlength*=n;
    
    return maxlength;
  }

#if HAVE_SUPERLU
  template<class T>
  template<class RowToDomain, class Solvers, class SubDomains>
  std::size_t SeqOverlappingSchwarzAssembler<SuperLU<T> >::assembleLocalProblems(const RowToDomain& rowToDomain, 
                                                                                const matrix_type& mat,
                                                                                Solvers& solvers,
                                                                                const SubDomains& subDomains,
                                                                                bool onTheFly)
  {
    typedef typename std::vector<SuperMatrixInitializer<matrix_type> >::iterator InitializerIterator;
    typedef typename SubDomains::const_iterator DomainIterator;
    typedef typename Solvers::iterator SolverIterator;
    std::size_t maxlength = 0;
    
    if(onTheFly){
      for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
	maxlength=std::max(maxlength, domain->size());
      maxlength*=mat[0].begin()->N();
    }else{
      // initialize the initializers
      DomainIterator domain=subDomains.begin();
    
      // Create the initializers list.
      std::vector<SuperMatrixInitializer<matrix_type> > initializers(subDomains.size());    
      
      SolverIterator solver=solvers.begin();
      for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
	  ++initializer, ++solver, ++domain){
	solver->mat.N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
	solver->mat.M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
	//solver->setVerbosity(true);
	*initializer=SuperMatrixInitializer<matrix_type>(solver->mat);
      }

      // Set up the supermatrices according to the subdomains
      typedef OverlappingSchwarzInitializer<std::vector<SuperMatrixInitializer<matrix_type> >,
	RowToDomain, SubDomains> Initializer;

      Initializer initializer(initializers, rowToDomain, subDomains);
      copyToSuperMatrix(initializer, mat);
      if(solvers.size()==1)
	assert(solvers[0].mat==mat);
    
      /*    for(SolverIterator solver=solvers.begin(); solver!=solvers.end(); ++solver)
	    dPrint_CompCol_Matrix("superlu", &static_cast<SuperMatrix&>(solver->mat)); */
      
      // Calculate the LU decompositions
      std::for_each(solvers.begin(), solvers.end(), std::mem_fun_ref(&SuperLU<matrix_type>::decompose));
      for(SolverIterator solver=solvers.begin(); solver!=solvers.end(); ++solver){
	assert(solver->mat.N()==solver->mat.M());
	maxlength=std::max(maxlength, solver->mat.N());
	//writeCompColMatrixToMatlab(static_cast<SuperLUMatrix<M>&>(solver->mat), std::cout);
      }
    }
    return maxlength;
  }

#endif

  template<class M,class X,class Y>
  template<class RowToDomain, class Solvers, class SubDomains>
  std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(const RowToDomain& rowToDomain, 
                                                                                                 const matrix_type& mat,
                                                                                                 Solvers& solvers,
                                                                                                 const SubDomains& subDomains,
                                                                                                 bool onTheFly)
  {
    typedef typename SubDomains::const_iterator DomainIterator;
    typedef typename Solvers::iterator SolverIterator;
    std::size_t maxlength = 0;
    
    if(onTheFly){
      for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
	maxlength=std::max(maxlength, domain->size());
    }else{
      // initialize the solvers of the local prolems.
      SolverIterator solver=solvers.begin();
      for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
          ++domain, ++solver){
        solver->setSubMatrix(mat, *domain);
	maxlength=std::max(maxlength, domain->size());
      }
    }
    
    return maxlength;
      
  }

  
  template<class M, class X, class TM, class TD, class TA>
  void SeqOverlappingSchwarz<M,X,TM,TD,TA>::apply(X& x, const X& b)
  {
    SeqOverlappingSchwarzApplier<SeqOverlappingSchwarz>::apply(*this, x, b);
  }
  
  template<class M, class X, class TM, class TD, class TA>
  template<bool forward>
  void SeqOverlappingSchwarz<M,X,TM,TD,TA>::apply(X& x, const X& b)
  {
    typedef typename X::block_type block;
    typedef slu_vector solver_vector;
    typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
    typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
      domain_iterator;
    
    OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
        
    domain_iterator domain=IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::begin(subDomains);
    iterator solver = IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::begin(solvers);
    X v(x); // temporary for the update
    v=0;
    
    typedef typename AdderSelector<TM,X,TD >::Adder Adder;    
    Adder adder(v, x, assigner, relax);
    
    nnz=0;
    std::size_t no=0;
    for(;domain != IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::end(subDomains); ++domain){
      //Copy rhs to C-array for SuperLU
      std::for_each(domain->begin(), domain->end(), assigner);
      assigner.resetIndexForNextDomain();
      if(onTheFly){
	// Create the subdomain solver
	slu sdsolver;
	sdsolver.setSubMatrix(mat, *domain);
	// Apply
	sdsolver.apply(assigner.lhs(), assigner.rhs());
        //nnz+=sdsolver.nnz();
      }else{
        solver->apply(assigner.lhs(), assigner.rhs());
        //nnz+=solver->nnz();
        ++solver;
      }
      ++no;
      //Add relaxed correction to from SuperLU to v
      std::for_each(domain->begin(), domain->end(), adder);
      assigner.resetIndexForNextDomain();
      
    }
    nnz/=no;
    
    adder.axpy();
    assigner.deallocate();
  }
  
  template<class K, int n, class Al, class X, class Y>
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, 
    const X& b_, Y& x_) :
    mat(&mat_), 
    rhs_( new DynamicVector<field_type>(maxlength, 42) ),
    lhs_( new DynamicVector<field_type>(maxlength, -42) ),
    b(&b_),
    x(&x_),
    i(0),
    maxlength_(maxlength)
  {}

  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::deallocate()
  {
    delete rhs_;
    delete lhs_;
  }
    
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::resetIndexForNextDomain()
  {
    i=0;
  }
  
  template<class K, int n, class Al, class X, class Y>
  DynamicVector<K> &
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::lhs()
  {
    return *lhs_;
  }
    
  template<class K, int n, class Al, class X, class Y>
  DynamicVector<K> &   
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::rhs()
  {
    return *rhs_;
  }
    
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::relaxResult(field_type relax)
  {
    lhs() *= relax;
  }
  
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::operator()(const size_type& domainIndex)
  {
    lhs() = 0.0;
#if 0
    //assign right hand side of current domainindex block
    for(size_type j=0; j<n; ++j, ++i){
      assert(i<maxlength_);
      rhs()[i]=(*b)[domainIndex][j];
    }
    
    // loop over all Matrix row entries and calculate defect.
    typedef typename matrix_type::ConstColIterator col_iterator;
    
    // calculate defect for current row index block
    for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
      block_type tmp(0.0);
      (*col).mv((*x)[col.index()], tmp);
      i-=n;
      for(size_type j=0; j<n; ++j, ++i){
        assert(i<maxlength_);
        rhs()[i]-=tmp[j];
      }
    }
#else
    //assign right hand side of current domainindex block
    for(size_type j=0; j<n; ++j, ++i){
      assert(i<maxlength_);
      rhs()[i]=(*b)[domainIndex][j];
      
      // loop over all Matrix row entries and calculate defect.
      typedef typename matrix_type::ConstColIterator col_iterator;
      
      // calculate defect for current row index block
      for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
        for(size_type k=0; k<n; ++k){
          rhs()[i]-=(*col)[j][k] * (*x)[col.index()][k];
        }
      }
    }
#endif 
  }

  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
  ::assignResult(block_type& res)
  {
    // assign the result of the local solve to the global vector
    for(size_type j=0; j<n; ++j, ++i){
      assert(i<maxlength_);
      res[j]+=lhs()[i];
    }
  }

#if HAVE_SUPERLU

  template<typename T, typename A, int n, int m>
  OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
  ::OverlappingAssigner(std::size_t maxlength, 
                        const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat_, 
                        const range_type& b_,
                        range_type& x_)
    : mat(&mat_), 
      b(&b_), 
      x(&x_), i(0), maxlength_(maxlength)
  {
    rhs_ = new field_type[maxlength];
    lhs_ = new field_type[maxlength];
    
  }
  
  template<typename T, typename A, int n, int m>
  void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::deallocate()
  {
    delete[] rhs_;
    delete[] lhs_;
  }
  
  template<typename T, typename A, int n, int m>
  void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::operator()(const size_type& domainIndex)
  {
    //assign right hand side of current domainindex block
    // rhs is an array of doubles!
    // rhs[starti] = b[domainindex]
    for(size_type j=0; j<n; ++j, ++i){
      assert(i<maxlength_);
      rhs_[i]=(*b)[domainIndex][j];
    }
    
    
    // loop over all Matrix row entries and calculate defect.
    typedef typename matrix_type::ConstColIterator col_iterator;

    // calculate defect for current row index block
    for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
        block_type tmp;
        (*col).mv((*x)[col.index()], tmp);
        i-=n;
        for(size_type j=0; j<n; ++j, ++i){
          assert(i<maxlength_);
          rhs_[i]-=tmp[j];
        }
        
    }
    
  }
  template<typename T, typename A, int n, int m>
  void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::relaxResult(field_type relax)
  {
    for(size_type j=i+n; i<j; ++i){
      assert(i<maxlength_);
      lhs_[i]*=relax;
    }
    i-=n;
  }
  
  template<typename T, typename A, int n, int m>
  void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::assignResult(block_type& res)
  {
    // assign the result of the local solve to the global vector
    for(size_type j=0; j<n; ++j, ++i){
      assert(i<maxlength_);
      res[j]+=lhs_[i];
    }
  }
  
  template<typename T, typename A, int n, int m>
  void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::resetIndexForNextDomain()
  {
    i=0;
  }
  
  template<typename T, typename A, int n, int m>
  typename OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::field_type* 
  OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::lhs()
  {
    return lhs_;
  }
  
  template<typename T, typename A, int n, int m>
  typename OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::field_type* 
  OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::rhs()
  {
    return rhs_;
  }

#endif

  template<class M, class X, class Y>
  OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength, 
                                                                        const M& mat_,
                                                                        const Y& b_,
                                                                        X& x_)
    : mat(&mat_), 
      b(&b_), 
      x(&x_), i(0)
  {
    rhs_= new Y(maxlength);
    lhs_ = new X(maxlength);
  }
  
  template<class M, class X, class Y>
  void OverlappingAssignerILUBase<M,X,Y>::deallocate()
  {
    delete rhs_;
    delete lhs_;
  }
  
  template<class M, class X, class Y>
  void OverlappingAssignerILUBase<M,X,Y>::operator()(const size_type& domainIndex)
  {
    (*rhs_)[i]=(*b)[domainIndex];
    
     // loop over all Matrix row entries and calculate defect.
    typedef typename matrix_type::ConstColIterator col_iterator;

    // calculate defect for current row index block
    for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
      (*col).mmv((*x)[col.index()], (*rhs_)[i]);
    }
    // Goto next local index
    ++i;
  }
  
  template<class M, class X, class Y>
  void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
  {
    (*lhs_)[i]*=relax;
  }
  
  template<class M, class X, class Y>
  void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
  {
    res+=(*lhs_)[i++];
  }

  template<class M, class X, class Y>
  X& OverlappingAssignerILUBase<M,X,Y>::lhs()
  {
    return *lhs_;
  }
  
  template<class M, class X, class Y>
  Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
  {
    return *rhs_;
  }

  template<class M, class X, class Y>
  void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
  {
    i=0;
  }
  
  template<typename S, typename T, typename A, int n>
  AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v_, 
                                                                   BlockVector<FieldVector<T,n>,A>& x_, 
                                                                   OverlappingAssigner<S>& assigner_, 
                                                                   const T& relax_)
    : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
  {}

  template<typename S, typename T, typename A, int n>
  void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
  {
    // add the result of the local solve to the current update
    assigner->assignResult((*v)[domainIndex]);
  }

 
  template<typename S, typename T, typename A, int n>
  void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
  {
    // relax the update and add it to the current guess.
    x->axpy(relax,*v);
  }

 
  template<typename S, typename T, typename A, int n>
  MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >
  ::MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v_, 
                        BlockVector<FieldVector<T,n>,A>& x_, 
                        OverlappingAssigner<S>& assigner_, const T& relax_)
    : x(&x_), assigner(&assigner_), relax(relax_)
  {}

  
  template<typename S,typename T, typename A, int n>
  void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
  {
    // add the result of the local solve to the current guess
    assigner->relaxResult(relax);
    assigner->assignResult((*x)[domainIndex]);
  }

  
  template<typename S,typename T, typename A, int n>
  void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
  {
    // nothing to do, as the corrections already relaxed and added in operator()
  }
  
  
  /** @} */
}

#endif