novlpschwarz.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
#define DUNE_ISTL_NOVLPSCHWARZ_HH
 
#include <iostream>              // for input/output to shell
#include <fstream>               // for input/output to files
#include <vector>                // STL vector class
#include <sstream>
 
#include <cmath>                // Yes, we do some math here
 
#include <dune/common/timer.hh>
 
#include "io.hh"
#include "bvector.hh"
#include "vbvector.hh"
#include "bcrsmatrix.hh"
#include "io.hh"
#include "gsetc.hh"
#include "ilu.hh"
#include "operators.hh"
#include "solvers.hh"
#include "preconditioners.hh"
#include "scalarproducts.hh"
#include "owneroverlapcopy.hh"
 
namespace Dune {
 
  template<class M, class X, class Y, class C>
  class NonoverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
  {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef C communication_type;
 
    typedef typename C::PIS PIS;
    typedef typename C::RI RI;
    typedef typename RI::RemoteIndexList RIL;
    typedef typename RI::const_iterator RIIterator;
    typedef typename RIL::const_iterator RILIterator;
    typedef typename M::ConstColIterator ColIterator;
    typedef typename M::ConstRowIterator RowIterator;
    typedef std::multimap<int,int> MM;
    typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
    typedef typename RIMap::iterator RIMapit;
 
    enum {
      category=SolverCategory::nonoverlapping
    };
 
    NonoverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
      : _A_(A), communication(com), buildcomm(true)
    {}
 
    virtual void apply (const X& x, Y& y) const
    {
      y = 0;
      novlp_op_apply(x,y,1);
      communication.addOwnerCopyToOwnerCopy(y,y);
    }
 
    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
    {
      // only apply communication to alpha*A*x to make it consistent,
      // y already has to be consistent.
      Y y1(y);
      y = 0;
      novlp_op_apply(x,y,alpha);
      communication.addOwnerCopyToOwnerCopy(y,y);
      y += y1;
    }
 
    virtual const matrix_type& getmat () const
    {
      return _A_;
    }
 
    void novlp_op_apply (const X& x, Y& y, field_type alpha) const
    {
      //get index sets
      const PIS& pis=communication.indexSet();
      const RI& ri = communication.remoteIndices();
 
      // at the beginning make a multimap "bordercontribution".
      // process has i and j as border dofs but is not the owner
      // => only contribute to Ax if i,j is in bordercontribution
      if (buildcomm == true) {
 
        // set up mask vector
        if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
          mask.resize(x.size());
          for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
            mask[i] = 1;
          for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
            if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
              mask[i->local().local()] = 0;
            else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
              mask[i->local().local()] = 2;
        }
 
        for (MM::iterator iter = bordercontribution.begin();
             iter != bordercontribution.end(); ++iter)
          bordercontribution.erase(iter);
        std::map<int,int> owner; //key: local index i, value: process, that owns i
        RIMap rimap;
 
        // for each local index make multimap rimap:
        // key: local index i, data: pair of process that knows i and pointer to RI entry
        for (RowIterator i = _A_.begin(); i != _A_.end(); ++i)
          if (mask[i.index()] == 0)
            for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
              RIL& ril = *(remote->second.first);
              for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
                if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
                  if (rindex->localIndexPair().local().local() == i.index()) {
                    rimap.insert
                      (std::make_pair(i.index(),
                                      std::pair<int,RILIterator>(remote->first, rindex)));
                    if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
                      owner.insert(std::make_pair(i.index(),remote->first));
                  }
            }
 
        int iowner = 0;
        for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
          if (mask[i.index()] == 0) {
            std::map<int,int>::iterator it = owner.find(i.index());
            iowner = it->second;
            std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
            for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
              if (mask[j.index()] == 0) {
                bool flag = true;
                for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
                  std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
                  for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
                    if (foundj->second.first == foundi->second.first)
                      if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
                          || foundj->second.first == iowner
                          || foundj->second.first  < communication.communicator().rank()) {
                        flag = false;
                        continue;
                      }
                  if (flag == false)
                    continue;
                }
                // don´t contribute to Ax if
                // 1. the owner of j has i as interior/border dof
                // 2. iowner has j as interior/border dof
                // 3. there is another process with smaller rank that has i and j
                // as interor/border dofs
                // if the owner of j does not have i as interior/border dof,
                // it will not be taken into account
                if (flag==true)
                  bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
              }
            }
          }
        }
        buildcomm = false;
      }
 
      //compute alpha*A*x nonoverlapping case
      for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
        if (mask[i.index()] == 0) {
          //dof doesn't belong to process but is border (not ghost)
          for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
            if (mask[j.index()] == 1) //j is owner => then sum entries
              (*j).usmv(alpha,x[j.index()],y[i.index()]);
            else if (mask[j.index()] == 0) {
              std::pair<MM::iterator, MM::iterator> itp =
                bordercontribution.equal_range(i.index());
              for (MM::iterator it = itp.first; it != itp.second; ++it)
                if ((*it).second == (int)j.index())
                  (*j).usmv(alpha,x[j.index()],y[i.index()]);
            }
          }
        }
        else if (mask[i.index()] == 1) {
          for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j)
            if (mask[j.index()] != 2)
              (*j).usmv(alpha,x[j.index()],y[i.index()]);
        }
      }
    }
 
  private:
    const matrix_type& _A_;
    const communication_type& communication;
    mutable bool buildcomm;
    mutable std::vector<double> mask;
    mutable std::multimap<int,int>  bordercontribution;
  };
 
  template<class X, class C>
  class NonoverlappingSchwarzScalarProduct : public ScalarProduct<X>
  {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
    typedef C communication_type;
 
    enum {category=SolverCategory::nonoverlapping};
 
    NonoverlappingSchwarzScalarProduct (const communication_type& com)
      : communication(com)
    {}
 
    virtual field_type dot (const X& x, const X& y)
    {
      field_type result;
      communication.dot(x,y,result);
      return result;
    }
 
    virtual real_type norm (const X& x)
    {
      return communication.norm(x);
    }
 
    void make_consistent (X& x) const
    {
      communication.copyOwnerToAll(x,x);
    }
 
  private:
    const communication_type& communication;
  };
 
  template<class X, class C>
  struct ScalarProductChooser<X,C,SolverCategory::nonoverlapping>
  {
    typedef NonoverlappingSchwarzScalarProduct<X,C> ScalarProduct;
    typedef C communication_type;
 
    enum {
      solverCategory=SolverCategory::nonoverlapping
    };
 
    static ScalarProduct* construct(const communication_type& comm)
    {
      return new ScalarProduct(comm);
    }
  };
 
  namespace Amg
  {
    template<class T> class ConstructionTraits;
  }
 
  template<class C, class P>
  class NonoverlappingBlockPreconditioner
    : public Dune::Preconditioner<typename P::domain_type,typename P::range_type> {
    friend class Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
  public:
    typedef typename P::domain_type domain_type;
    typedef typename P::range_type range_type;
    typedef C communication_type;
 
    // define the category
    enum {
      category=SolverCategory::nonoverlapping
    };
 
    NonoverlappingBlockPreconditioner (P& prec, const communication_type& c)
      : preconditioner(prec), communication(c)
    {}
 
    virtual void pre (domain_type& x, range_type& b)
    {
      preconditioner.pre(x,b);
    }
 
    virtual void apply (domain_type& v, const range_type& d)
    {
      // block preconditioner equivalent to WrappedPreconditioner from
      // pdelab/backend/ovlpistsolverbackend.hh,
      // but not to BlockPreconditioner from schwarz.hh
      preconditioner.apply(v,d);
      communication.addOwnerCopyToOwnerCopy(v,v);
    }
 
    virtual void post (domain_type& x)
    {
      preconditioner.post(x);
    }
 
  private:
    P& preconditioner;
 
    const communication_type& communication;
  };
 
} // end namespace
 
#endif