overlappingschwarz.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
#define DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
#include <cassert>
#include <algorithm>
#include <functional>
#include <vector>
#include <set>
#include <dune/common/dynmatrix.hh>
#include <dune/common/sllist.hh>
#include <dune/common/unused.hh>
#include "preconditioners.hh"
#include "superlu.hh"
#include "umfpack.hh"
#include "bvector.hh"
#include "bcrsmatrix.hh"
#include "ilusubdomainsolver.hh"
#include <dune/istl/solvertype.hh>
 
namespace Dune
{
 
  template<class M, class X, class TM, class TD, class TA>
  class SeqOverlappingSchwarz;
 
  template<class I, class S, class D>
  class OverlappingSchwarzInitializer
  {
  public:
    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;
  };
 
  struct AdditiveSchwarzMode
  {};
 
  struct MultiplicativeSchwarzMode
  {};
 
  struct SymmetricMultiplicativeSchwarzMode
  {};
 
  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:
    typedef typename std::remove_const<M>::type matrix_type;
    typedef K field_type;
    typedef typename std::remove_const<M>::type rilu_type;
    typedef X domain_type;
    typedef Y range_type;
 
    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());
    }
 
    template<class S>
    void setSubMatrix(const M& BCRS, S& rowset)
    {
      size_t sz = rowset.size();
      A.resize(sz*n,sz*n);
      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, bool tag>
  class OverlappingAssignerHelper
  {};
 
  template<typename T>
  using OverlappingAssigner = OverlappingAssignerHelper<T, Dune::StoresColumnCompressed<T>::value>;
 
  // specialization for DynamicMatrix
  template<class K, int n, class Al, class X, class Y>
  class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  {
  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;
 
    OverlappingAssignerHelper(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, const X& b_, Y& x_);
 
    inline
    void deallocate();
 
    inline
    void resetIndexForNextDomain();
 
    inline
    DynamicVector<K> & lhs();
 
    inline
    DynamicVector<K> & rhs();
 
    inline
    void relaxResult(field_type relax);
 
    void operator()(const size_type& domainIndex);
 
    inline
    void assignResult(block_type& res);
 
  private:
    const matrix_type* mat;
    // we need a pointer, because we have to avoid deep copies
    DynamicVector<field_type> * rhs_;
    // we need a pointer, because we have to avoid deep copies
    DynamicVector<field_type> * lhs_;
    const range_type* b;
    range_type* x;
    std::size_t i;
    std::size_t maxlength_;
  };
 
#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
  template<template<class> class S, int n, int m, typename T, typename A>
  struct OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>, A> >, true>
  {
    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> matrix_type;
    typedef typename S<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;
 
    OverlappingAssignerHelper(std::size_t maxlength, const matrix_type& mat,
                        const range_type& b, range_type& x);
    void deallocate();
 
    /*
     * @brief Resets the local index to zero.
     */
    void resetIndexForNextDomain();
 
    field_type* lhs();
 
    field_type* rhs();
 
    void relaxResult(field_type relax);
 
    void operator()(const size_type& domain);
 
    void assignResult(block_type& res);
 
  private:
    const matrix_type* mat;
    field_type* rhs_;
    field_type* lhs_;
    const range_type* b;
    range_type* x;
    std::size_t i;
    std::size_t maxlength_;
  };
 
#endif // HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
 
  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;
    OverlappingAssignerILUBase(std::size_t maxlength, const M& mat,
                               const Y& b, X& x);
    void deallocate();
 
    void resetIndexForNextDomain();
 
    X& lhs();
 
    Y& rhs();
 
    void relaxResult(field_type relax);
 
    void operator()(const size_type& domain);
 
    void assignResult(block_type& res);
 
  private:
    const M* mat;
    X* lhs_;
    Y* rhs_;
    const Y* b;
    X* x;
    size_type i;
  };
 
  // specialization for ILU0
  template<class M, class X, class Y>
  class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
    : public OverlappingAssignerILUBase<M,X,Y>
  {
  public:
    OverlappingAssignerHelper(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 OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
    : public OverlappingAssignerILUBase<M,X,Y>
  {
  public:
    OverlappingAssignerHelper(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;
  };
 
  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;
  };
 
  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();
    }
  };
 
  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, bool tag>
  struct SeqOverlappingSchwarzAssemblerHelper
  {};
 
  template<class T>
  using SeqOverlappingSchwarzAssembler = SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>;
 
  template<class K, int n, class Al, class X, class Y>
  struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  {
    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);
  };
 
  template<template<class> class S, typename T, typename A, int m, int n>
  struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,true>
  {
    typedef BCRSMatrix<FieldMatrix<T,m,n>,A> 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 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 SeqOverlappingSchwarzAssemblerHelper<ILU0SubdomainSolver<M,X,Y>,false>
    : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
  {};
 
  template<class M,class X, class Y>
  struct SeqOverlappingSchwarzAssemblerHelper<ILUNSubdomainSolver<M,X,Y>,false>
    : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
  {};
 
  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:
    typedef M matrix_type;
 
    typedef X domain_type;
 
    typedef X range_type;
 
    typedef TM Mode;
 
    typedef typename X::field_type field_type;
 
    typedef typename matrix_type::size_type size_type;
 
    typedef TA allocator;
 
    typedef std::set<size_type, std::less<size_type>,
        typename TA::template rebind<size_type>::other>
    subdomain_type;
 
    typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other> subdomain_vector;
 
    typedef SLList<size_type, typename TA::template rebind<size_type>::other> subdomain_list;
 
    typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other > rowtodomain_vector;
 
    typedef TD slu;
 
    typedef std::vector<slu, typename TA::template rebind<slu>::other> slu_vector;
 
    enum {
      category = SolverCategory::sequential
    };
 
    SeqOverlappingSchwarz(const matrix_type& mat, const subdomain_vector& subDomains,
                          field_type relaxationFactor=1, bool onTheFly_=true);
 
    SeqOverlappingSchwarz(const matrix_type& mat, const rowtodomain_vector& rowToDomain,
                          field_type relaxationFactor=1, bool onTheFly_=true);
 
    virtual void pre (X& x, X& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(b);
    }
 
    virtual void apply (X& v, const X& d);
 
    virtual void post (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
    template<bool forward>
    void apply(X& v, const X& d);
 
  private:
    const M& mat;
    slu_vector solvers;
    subdomain_vector subDomains;
    field_type relax;
 
    typename M::size_type maxlength;
 
    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
  {
    std::vector<IndexMap>().swap(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<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
  SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>::
  assembleLocalProblems(const RowToDomain& rowToDomain,
                        const matrix_type& mat,
                        Solvers& solvers,
                        const SubDomains& subDomains,
                        bool onTheFly)
  {
    DUNE_UNUSED_PARAMETER(onTheFly);
    DUNE_UNUSED_PARAMETER(rowToDomain);
    DUNE_UNUSED_PARAMETER(mat);
    DUNE_UNUSED_PARAMETER(solvers);
    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 || HAVE_SUITESPARSE_UMFPACK
  template<template<class> class S, typename T, typename A, int m, int n>
  template<class RowToDomain, class Solvers, class SubDomains>
  std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,true>::assembleLocalProblems(const RowToDomain& rowToDomain,
                                                                                 const matrix_type& mat,
                                                                                 Solvers& solvers,
                                                                                 const SubDomains& subDomains,
                                                                                 bool onTheFly)
  {
    typedef typename S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::MatrixInitializer MatrixInitializer;
    typedef typename std::vector<MatrixInitializer>::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<MatrixInitializer> initializers(subDomains.size());
 
      SolverIterator solver=solvers.begin();
      for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
          ++initializer, ++solver, ++domain) {
        solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
        solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
        //solver->setVerbosity(true);
        *initializer=MatrixInitializer(solver->getInternalMatrix());
      }
 
      // Set up the supermatrices according to the subdomains
      typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
          RowToDomain, SubDomains> Initializer;
 
      Initializer initializer(initializers, rowToDomain, subDomains);
      copyToColCompMatrix(initializer, mat);
 
      // Calculate the LU decompositions
      std::for_each(solvers.begin(), solvers.end(), std::mem_fun_ref(&S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::decompose));
      for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
      {
        assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
        maxlength = std::max(maxlength, solverIt->getInternalMatrix().N());
      }
    }
    return maxlength;
  }
 
#endif // HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
 
  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)
  {
    DUNE_UNUSED_PARAMETER(rowToDomain);
    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 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);
 
    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());
      }else{
        solver->apply(assigner.lhs(), assigner.rhs());
        ++solver;
      }
 
      //Add relaxed correction to from SuperLU to v
      std::for_each(domain->begin(), domain->end(), adder);
      assigner.resetIndexForNextDomain();
 
    }
 
    adder.axpy();
    assigner.deallocate();
  }
 
  template<class K, int n, class Al, class X, class Y>
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::OverlappingAssignerHelper(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
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::deallocate()
  {
    delete rhs_;
    delete lhs_;
  }
 
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::resetIndexForNextDomain()
  {
    i=0;
  }
 
  template<class K, int n, class Al, class X, class Y>
  DynamicVector<K> &
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::lhs()
  {
    return *lhs_;
  }
 
  template<class K, int n, class Al, class X, class Y>
  DynamicVector<K> &
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::rhs()
  {
    return *rhs_;
  }
 
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::relaxResult(field_type relax)
  {
    lhs() *= relax;
  }
 
  template<class K, int n, class Al, class X, class Y>
  void
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::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
  OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
  ::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 || HAVE_SUITESPARSE_UMFPACK
 
  template<template<class> class S, int n, int m, typename T, typename A>
  OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>
  ::OverlappingAssignerHelper(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<template<class> class S, int n, int m, typename T, typename A>
  void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::deallocate()
  {
    delete[] rhs_;
    delete[] lhs_;
  }
 
  template<template<class> class S, int n, int m, typename T, typename A>
  void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::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<template<class> class S, int n, int m, typename T, typename A>
  void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::relaxResult(field_type relax)
  {
    for(size_type j=i+n; i<j; ++i) {
      assert(i<maxlength_);
      lhs_[i]*=relax;
    }
    i-=n;
  }
 
  template<template<class> class S, int n, int m, typename T, typename A>
  void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::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<template<class> class S, int n, int m, typename T, typename A>
  void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::resetIndexForNextDomain()
  {
    i=0;
  }
 
  template<template<class> class S, int n, int m, typename T, typename A>
  typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::field_type*
  OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::lhs()
  {
    return lhs_;
  }
 
  template<template<class> class S, int n, int m, typename T, typename A>
  typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::field_type*
  OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::rhs()
  {
    return rhs_;
  }
 
#endif // HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
 
  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_)
  {
    DUNE_UNUSED_PARAMETER(v_);
  }
 
 
  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