hierarchy.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_AMGHIERARCHY_HH
#define DUNE_AMGHIERARCHY_HH
 
#include <list>
#include <memory>
#include <limits>
#include <algorithm>
#include <tuple>
#include "aggregates.hh"
#include "graph.hh"
#include "galerkin.hh"
#include "renumberer.hh"
#include "graphcreator.hh"
#include <dune/common/stdstreams.hh>
#include <dune/common/unused.hh>
#include <dune/common/timer.hh>
#include <dune/common/bigunsignedint.hh>
#include <dune/istl/bvector.hh>
#include <dune/common/parallel/indexset.hh>
#include <dune/istl/matrixutils.hh>
#include <dune/istl/matrixredistribute.hh>
#include <dune/istl/paamg/dependency.hh>
#include <dune/istl/paamg/graph.hh>
#include <dune/istl/paamg/indicescoarsener.hh>
#include <dune/istl/paamg/globalaggregates.hh>
#include <dune/istl/paamg/construction.hh>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/transfer.hh>
 
namespace Dune
{
  namespace Amg
  {
    enum {
      MAX_PROCESSES = 72000
    };
 
    template<typename T, typename A=std::allocator<T> >
    class Hierarchy
    {
    public:
      typedef T MemberType;
 
      template<typename T1, typename T2>
      class LevelIterator;
 
    private:
      struct Element
      {
        friend class LevelIterator<Hierarchy<T,A>, T>;
        friend class LevelIterator<const Hierarchy<T,A>, const T>;
 
        Element* coarser_;
 
        Element* finer_;
 
        MemberType* element_;
 
        MemberType* redistributed_;
      };
    public:
      //       enum{
      //        /**
      //         * @brief If true only the method addCoarser will be usable
      //         * otherwise only the method addFiner will be usable.
      //         */
      //        coarsen = b
      //          };
 
      typedef typename A::template rebind<Element>::other Allocator;
 
      typedef typename ConstructionTraits<T>::Arguments Arguments;
 
      Hierarchy(MemberType& first);
 
      Hierarchy(MemberType* first);
 
      Hierarchy();
 
      Hierarchy(const Hierarchy& other);
      void addCoarser(Arguments& args);
 
      void addRedistributedOnCoarsest(Arguments& args);
 
      void addFiner(Arguments& args);
 
      template<class C, class T1>
      class LevelIterator
        : public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
      {
        friend class LevelIterator<typename std::remove_const<C>::type,
            typename std::remove_const<T1>::type >;
        friend class LevelIterator<const typename std::remove_const<C>::type,
            const typename std::remove_const<T1>::type >;
 
      public:
        LevelIterator()
          : element_(0)
        {}
 
        LevelIterator(Element* element)
          : element_(element)
        {}
 
        LevelIterator(const LevelIterator<typename std::remove_const<C>::type,
                          typename std::remove_const<T1>::type>& other)
          : element_(other.element_)
        {}
 
        LevelIterator(const LevelIterator<const typename std::remove_const<C>::type,
                          const typename std::remove_const<T1>::type>& other)
          : element_(other.element_)
        {}
 
        bool equals(const LevelIterator<typename std::remove_const<C>::type,
                        typename std::remove_const<T1>::type>& other) const
        {
          return element_ == other.element_;
        }
 
        bool equals(const LevelIterator<const typename std::remove_const<C>::type,
                        const typename std::remove_const<T1>::type>& other) const
        {
          return element_ == other.element_;
        }
 
        T1& dereference() const
        {
          return *(element_->element_);
        }
 
        void increment()
        {
          element_ = element_->coarser_;
        }
 
        void decrement()
        {
          element_ = element_->finer_;
        }
 
        bool isRedistributed() const
        {
          return element_->redistributed_;
        }
 
        T1& getRedistributed() const
        {
          assert(element_->redistributed_);
          return *element_->redistributed_;
        }
        void addRedistributed(T1* t)
        {
          element_->redistributed_ = t;
        }
 
        void deleteRedistributed()
        {
          element_->redistributed_ = nullptr;
        }
 
      private:
        Element* element_;
      };
 
      typedef LevelIterator<Hierarchy<T,A>,T> Iterator;
 
      typedef LevelIterator<const Hierarchy<T,A>, const T> ConstIterator;
 
      Iterator finest();
 
      Iterator coarsest();
 
 
      ConstIterator finest() const;
 
      ConstIterator coarsest() const;
 
      std::size_t levels() const;
 
      ~Hierarchy();
 
    private:
      Element* finest_;
      Element* coarsest_;
      Element* nonAllocated_;
      Allocator allocator_;
      int levels_;
    };
 
    template<class M, class PI, class A=std::allocator<M> >
    class MatrixHierarchy
    {
    public:
      typedef M MatrixOperator;
 
      typedef typename MatrixOperator::matrix_type Matrix;
 
      typedef PI ParallelInformation;
 
      typedef A Allocator;
 
      typedef Dune::Amg::AggregatesMap<typename MatrixGraph<Matrix>::VertexDescriptor> AggregatesMap;
 
      typedef Dune::Amg::Hierarchy<MatrixOperator,Allocator> ParallelMatrixHierarchy;
 
      typedef Dune::Amg::Hierarchy<ParallelInformation,Allocator> ParallelInformationHierarchy;
 
      typedef typename Allocator::template rebind<AggregatesMap*>::other AAllocator;
 
      typedef std::list<AggregatesMap*,AAllocator> AggregatesMapList;
 
      typedef RedistributeInformation<ParallelInformation> RedistributeInfoType;
 
      typedef typename Allocator::template rebind<RedistributeInfoType>::other RILAllocator;
 
      typedef std::list<RedistributeInfoType,RILAllocator> RedistributeInfoList;
 
      MatrixHierarchy(const MatrixOperator& fineMatrix,
                      const ParallelInformation& pinfo=ParallelInformation());
 
 
      ~MatrixHierarchy();
 
      template<typename O, typename T>
      void build(const T& criterion);
 
      template<class F>
      void recalculateGalerkin(const F& copyFlags);
 
      template<class V, class TA>
      void coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const;
 
      template<class S, class TA>
      void coarsenSmoother(Hierarchy<S,TA>& smoothers,
                           const typename SmootherTraits<S>::Arguments& args) const;
 
      std::size_t levels() const;
 
      std::size_t maxlevels() const;
 
      bool hasCoarsest() const;
 
      bool isBuilt() const;
 
      const ParallelMatrixHierarchy& matrices() const;
 
      const ParallelInformationHierarchy& parallelInformation() const;
 
      const AggregatesMapList& aggregatesMaps() const;
 
      const RedistributeInfoList& redistributeInformation() const;
 
 
      typename MatrixOperator::field_type getProlongationDampingFactor() const
      {
        return prolongDamp_;
      }
 
      void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const;
 
    private:
      typedef typename ConstructionTraits<MatrixOperator>::Arguments MatrixArgs;
      typedef typename ConstructionTraits<ParallelInformation>::Arguments CommunicationArgs;
      AggregatesMapList aggregatesMaps_;
      RedistributeInfoList redistributes_;
      ParallelMatrixHierarchy matrices_;
      ParallelInformationHierarchy parallelInformation_;
 
      bool built_;
 
      int maxlevels_;
 
      typename MatrixOperator::field_type prolongDamp_;
 
      template<class Matrix, bool print>
      struct MatrixStats
      {
 
        static void stats(const Matrix& matrix)
        {
          DUNE_UNUSED_PARAMETER(matrix);
        }
      };
 
      template<class Matrix>
      struct MatrixStats<Matrix,true>
      {
        struct calc
        {
          typedef typename Matrix::size_type size_type;
          typedef typename Matrix::row_type matrix_row;
 
          calc()
          {
            min=std::numeric_limits<size_type>::max();
            max=0;
            sum=0;
          }
 
          void operator()(const matrix_row& row)
          {
            min=std::min(min, row.size());
            max=std::max(max, row.size());
            sum += row.size();
          }
 
          size_type min;
          size_type max;
          size_type sum;
        };
        static void stats(const Matrix& matrix)
        {
          calc c= for_each(matrix.begin(), matrix.end(), calc());
          dinfo<<"Matrix row: min="<<c.min<<" max="<<c.max
               <<" average="<<static_cast<double>(c.sum)/matrix.N()
               <<std::endl;
        }
      };
    };
 
    template<class T>
    class CoarsenCriterion : public T
    {
    public:
      typedef T AggregationCriterion;
 
      CoarsenCriterion(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
                       double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
        : AggregationCriterion(Dune::Amg::Parameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate))
      {}
 
      CoarsenCriterion(const Dune::Amg::Parameters& parms)
        : AggregationCriterion(parms)
      {}
 
    };
 
    template<typename M, typename C1>
    bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix,
                                        SequentialInformation& origSequentialInformationomm,
                                        SequentialInformation*& newComm,
                                        RedistributeInformation<SequentialInformation>& ri,
                                        int nparts, C1& criterion)
    {
      DUNE_UNUSED_PARAMETER(origMatrix);
      DUNE_UNUSED_PARAMETER(newMatrix);
      DUNE_UNUSED_PARAMETER(origSequentialInformationomm);
      DUNE_UNUSED_PARAMETER(newComm);
      DUNE_UNUSED_PARAMETER(ri);
      DUNE_UNUSED_PARAMETER(nparts);
      DUNE_UNUSED_PARAMETER(criterion);
      DUNE_THROW(NotImplemented, "Redistribution does not make sense in sequential code!");
    }
 
 
    template<typename M, typename C, typename C1>
    bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix, C& origComm, C*& newComm,
                                        RedistributeInformation<C>& ri,
                                        int nparts, C1& criterion)
    {
      Timer time;
#ifdef AMG_REPART_ON_COMM_GRAPH
      // Done not repartition the matrix graph, but a graph of the communication scheme.
      bool existentOnRedist=Dune::commGraphRepartition(origMatrix, origComm, nparts, newComm,
                                                       ri.getInterface(),
                                                       criterion.debugLevel()>1);
 
#else
      typedef Dune::Amg::MatrixGraph<const M> MatrixGraph;
      typedef Dune::Amg::PropertiesGraph<MatrixGraph,
          VertexProperties,
          EdgeProperties,
          IdentityMap,
          IdentityMap> PropertiesGraph;
      MatrixGraph graph(origMatrix);
      PropertiesGraph pgraph(graph);
      buildDependency(pgraph, origMatrix, criterion, false);
 
#ifdef DEBUG_REPART
      if(origComm.communicator().rank()==0)
        std::cout<<"Original matrix"<<std::endl;
      origComm.communicator().barrier();
      printGlobalSparseMatrix(origMatrix, origComm, std::cout);
#endif
      bool existentOnRedist=Dune::graphRepartition(pgraph, origComm, nparts,
                                                   newComm, ri.getInterface(),
                                                   criterion.debugLevel()>1);
#endif // if else AMG_REPART
 
      if(origComm.communicator().rank()==0  && criterion.debugLevel()>1)
        std::cout<<"Repartitioning took "<<time.elapsed()<<" seconds."<<std::endl;
 
      ri.setSetup();
 
#ifdef DEBUG_REPART
      ri.checkInterface(origComm.indexSet(), newComm->indexSet(), origComm.communicator());
#endif
 
      redistributeMatrix(const_cast<M&>(origMatrix), newMatrix, origComm, *newComm, ri);
 
#ifdef DEBUG_REPART
      if(origComm.communicator().rank()==0)
        std::cout<<"Original matrix"<<std::endl;
      origComm.communicator().barrier();
      if(newComm->communicator().size()>0)
        printGlobalSparseMatrix(newMatrix, *newComm, std::cout);
      origComm.communicator().barrier();
#endif
 
      if(origComm.communicator().rank()==0  && criterion.debugLevel()>1)
        std::cout<<"Redistributing matrix took "<<time.elapsed()<<" seconds."<<std::endl;
      return existentOnRedist;
 
    }
 
    template<typename M>
    bool repartitionAndDistributeMatrix(M& origMatrix, M& newMatrix,
                                        SequentialInformation& origComm,
                                        SequentialInformation& newComm,
                                        RedistributeInformation<SequentialInformation>& ri)
    {
      return true;
    }
 
    template<class M, class IS, class A>
    MatrixHierarchy<M,IS,A>::MatrixHierarchy(const MatrixOperator& fineOperator,
                                             const ParallelInformation& pinfo)
      : matrices_(const_cast<MatrixOperator&>(fineOperator)),
        parallelInformation_(const_cast<ParallelInformation&>(pinfo))
    {
      static_assert((static_cast<int>(MatrixOperator::category) ==
                       static_cast<int>(SolverCategory::sequential)
                     || static_cast<int>(MatrixOperator::category) ==
                       static_cast<int>(SolverCategory::overlapping)
                     || static_cast<int>(MatrixOperator::category) ==
                       static_cast<int>(SolverCategory::nonoverlapping)),
                    "MatrixOperator must be of category sequential or overlapping or nonoverlapping");
      if (static_cast<int>(MatrixOperator::category) != static_cast<int>(pinfo.getSolverCategory()))
        DUNE_THROW(ISTLError, "MatrixOperator and ParallelInformation must belong to the same category!");
 
    }
 
    template<class M, class IS, class A>
    template<typename O, typename T>
    void MatrixHierarchy<M,IS,A>::build(const T& criterion)
    {
      prolongDamp_ = criterion.getProlongationDampingFactor();
      typedef O OverlapFlags;
      typedef typename ParallelMatrixHierarchy::Iterator MatIterator;
      typedef typename ParallelInformationHierarchy::Iterator PInfoIterator;
 
      static const int noints=(Dune::Amg::MAX_PROCESSES/4096>0) ? (Dune::Amg::MAX_PROCESSES/4096) : 1;
 
      typedef bigunsignedint<sizeof(int)*8*noints> BIGINT;
      GalerkinProduct<ParallelInformation> productBuilder;
      MatIterator mlevel = matrices_.finest();
      MatrixStats<typename M::matrix_type,MINIMAL_DEBUG_LEVEL<=INFO_DEBUG_LEVEL>::stats(mlevel->getmat());
 
      PInfoIterator infoLevel = parallelInformation_.finest();
      BIGINT finenonzeros=countNonZeros(mlevel->getmat());
      finenonzeros = infoLevel->communicator().sum(finenonzeros);
      BIGINT allnonzeros = finenonzeros;
 
 
      int level = 0;
      int rank = 0;
 
      BIGINT unknowns = mlevel->getmat().N();
 
      unknowns = infoLevel->communicator().sum(unknowns);
      double dunknowns=unknowns.todouble();
      infoLevel->buildGlobalLookup(mlevel->getmat().N());
      redistributes_.push_back(RedistributeInfoType());
 
      for(; level < criterion.maxLevel(); ++level, ++mlevel) {
        assert(matrices_.levels()==redistributes_.size());
        rank = infoLevel->communicator().rank();
        if(rank==0 && criterion.debugLevel()>1)
          std::cout<<"Level "<<level<<" has "<<dunknowns<<" unknowns, "<<dunknowns/infoLevel->communicator().size()
                   <<" unknowns per proc (procs="<<infoLevel->communicator().size()<<")"<<std::endl;
 
        MatrixOperator* matrix=&(*mlevel);
        ParallelInformation* info =&(*infoLevel);
 
        if((
#if HAVE_PARMETIS
             criterion.accumulate()==successiveAccu
#else
             false
#endif
             || (criterion.accumulate()==atOnceAccu
                 && dunknowns < 30*infoLevel->communicator().size()))
           && infoLevel->communicator().size()>1 &&
           dunknowns/infoLevel->communicator().size() <= criterion.coarsenTarget())
        {
          // accumulate to fewer processors
          Matrix* redistMat= new Matrix();
          ParallelInformation* redistComm=0;
          std::size_t nodomains = (std::size_t)std::ceil(dunknowns/(criterion.minAggregateSize()
                                                                    *criterion.coarsenTarget()));
          if( nodomains<=criterion.minAggregateSize()/2 ||
              dunknowns <= criterion.coarsenTarget() )
            nodomains=1;
 
          bool existentOnNextLevel =
            repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
                                           redistComm, redistributes_.back(), nodomains,
                                           criterion);
          BIGINT unknownsRedist = redistMat->N();
          unknownsRedist = infoLevel->communicator().sum(unknownsRedist);
          dunknowns= unknownsRedist.todouble();
          if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1)
            std::cout<<"Level "<<level<<" (redistributed) has "<<dunknowns<<" unknowns, "<<dunknowns/redistComm->communicator().size()
                     <<" unknowns per proc (procs="<<redistComm->communicator().size()<<")"<<std::endl;
          MatrixArgs args(*redistMat, *redistComm);
          mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
          assert(mlevel.isRedistributed());
          infoLevel.addRedistributed(redistComm);
          infoLevel->freeGlobalLookup();
 
          if(!existentOnNextLevel)
            // We do not hold any data on the redistributed partitioning
            break;
 
          // Work on the redistributed Matrix from now on
          matrix = &(mlevel.getRedistributed());
          info = &(infoLevel.getRedistributed());
          info->buildGlobalLookup(matrix->getmat().N());
        }
 
        rank = info->communicator().rank();
        if(dunknowns <= criterion.coarsenTarget())
          // No further coarsening needed
          break;
 
        typedef PropertiesGraphCreator<MatrixOperator> GraphCreator;
        typedef typename GraphCreator::PropertiesGraph PropertiesGraph;
        typedef typename GraphCreator::GraphTuple GraphTuple;
 
        typedef typename PropertiesGraph::VertexDescriptor Vertex;
 
        std::vector<bool> excluded(matrix->getmat().N(), false);
 
        GraphTuple graphs = GraphCreator::create(*matrix, excluded, *info, OverlapFlags());
 
        AggregatesMap* aggregatesMap=new AggregatesMap(std::get<1>(graphs)->maxVertex()+1);
 
        aggregatesMaps_.push_back(aggregatesMap);
 
        Timer watch;
        watch.reset();
        int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
 
        std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
          aggregatesMap->buildAggregates(matrix->getmat(), *(std::get<1>(graphs)), criterion, level==0);
 
        if(rank==0 && criterion.debugLevel()>2)
          std::cout<<" Have built "<<noAggregates<<" aggregates totally ("<<isoAggregates<<" isolated aggregates, "<<
          oneAggregates<<" aggregates of one vertex,  and skipped "<<
          skippedAggregates<<" aggregates)."<<std::endl;
#ifdef TEST_AGGLO
        {
          // calculate size of local matrix in the distributed direction
          int start, end, overlapStart, overlapEnd;
          int procs=info->communicator().rank();
          int n = UNKNOWNS/procs; // number of unknowns per process
          int bigger = UNKNOWNS%procs; // number of process with n+1 unknows
 
          // Compute owner region
          if(rank<bigger) {
            start = rank*(n+1);
            end   = (rank+1)*(n+1);
          }else{
            start = bigger + rank * n;
            end   = bigger + (rank + 1) * n;
          }
 
          // Compute overlap region
          if(start>0)
            overlapStart = start - 1;
          else
            overlapStart = start;
 
          if(end<UNKNOWNS)
            overlapEnd = end + 1;
          else
            overlapEnd = end;
 
          assert((UNKNOWNS)*(overlapEnd-overlapStart)==aggregatesMap->noVertices());
          for(int j=0; j< UNKNOWNS; ++j)
            for(int i=0; i < UNKNOWNS; ++i)
            {
              if(i>=overlapStart && i<overlapEnd)
              {
                int no = (j/2)*((UNKNOWNS)/2)+i/2;
                (*aggregatesMap)[j*(overlapEnd-overlapStart)+i-overlapStart]=no;
              }
            }
        }
#endif
        if(criterion.debugLevel()>1 && info->communicator().rank()==0)
          std::cout<<"aggregating finished."<<std::endl;
 
        BIGINT gnoAggregates=noAggregates;
        gnoAggregates = info->communicator().sum(gnoAggregates);
        double dgnoAggregates = gnoAggregates.todouble();
#ifdef TEST_AGGLO
        BIGINT gnoAggregates=((UNKNOWNS)/2)*((UNKNOWNS)/2);
#endif
 
        if(criterion.debugLevel()>2 && rank==0)
          std::cout << "Building "<<dgnoAggregates<<" aggregates took "<<watch.elapsed()<<" seconds."<<std::endl;
 
        if(dgnoAggregates==0 || dunknowns/dgnoAggregates<criterion.minCoarsenRate())
        {
          if(rank==0)
          {
            if(dgnoAggregates>0)
              std::cerr << "Stopped coarsening because of rate breakdown "<<dunknowns<<"/"<<dgnoAggregates
                        <<"="<<dunknowns/dgnoAggregates<<"<"
                        <<criterion.minCoarsenRate()<<std::endl;
            else
              std::cerr<< "Could not build any aggregates. Probably no connected nodes."<<std::endl;
          }
          aggregatesMap->free();
          delete aggregatesMap;
          aggregatesMaps_.pop_back();
 
          if(criterion.accumulate() && mlevel.isRedistributed() && info->communicator().size()>1) {
            // coarse level matrix was already redistributed, but to more than 1 process
            // Therefore need to delete the redistribution. Further down it will
            // then be redistributed to 1 process
            delete &(mlevel.getRedistributed().getmat());
            mlevel.deleteRedistributed();
            delete &(infoLevel.getRedistributed());
            infoLevel.deleteRedistributed();
            redistributes_.back().resetSetup();
          }
 
          break;
        }
        unknowns =  noAggregates;
        dunknowns = dgnoAggregates;
 
        CommunicationArgs commargs(info->communicator(),info->getSolverCategory());
        parallelInformation_.addCoarser(commargs);
 
        ++infoLevel; // parallel information on coarse level
 
        typename PropertyMapTypeSelector<VertexVisitedTag,PropertiesGraph>::Type visitedMap =
          get(VertexVisitedTag(), *(std::get<1>(graphs)));
 
        watch.reset();
        int aggregates = IndicesCoarsener<ParallelInformation,OverlapFlags>
                         ::coarsen(*info,
                                   *(std::get<1>(graphs)),
                                   visitedMap,
                                   *aggregatesMap,
                                   *infoLevel,
                                   noAggregates);
        GraphCreator::free(graphs);
 
        if(criterion.debugLevel()>2) {
          if(rank==0)
            std::cout<<"Coarsening of index sets took "<<watch.elapsed()<<" seconds."<<std::endl;
        }
 
        watch.reset();
 
        infoLevel->buildGlobalLookup(aggregates);
        AggregatesPublisher<Vertex,OverlapFlags,ParallelInformation>::publish(*aggregatesMap,
                                                                              *info,
                                                                              infoLevel->globalLookup());
 
 
        if(criterion.debugLevel()>2) {
          if(rank==0)
            std::cout<<"Communicating global aggregate numbers took "<<watch.elapsed()<<" seconds."<<std::endl;
        }
 
        watch.reset();
        std::vector<bool>& visited=excluded;
 
        typedef std::vector<bool>::iterator Iterator;
        typedef IteratorPropertyMap<Iterator, IdentityMap> VisitedMap2;
        Iterator end = visited.end();
        for(Iterator iter= visited.begin(); iter != end; ++iter)
          *iter=false;
 
        VisitedMap2 visitedMap2(visited.begin(), Dune::IdentityMap());
 
        typename MatrixOperator::matrix_type* coarseMatrix;
 
        coarseMatrix = productBuilder.build(*(std::get<0>(graphs)), visitedMap2,
                                            *info,
                                            *aggregatesMap,
                                            aggregates,
                                            OverlapFlags());
        dverb<<"Building of sparsity pattern took "<<watch.elapsed()<<std::endl;
        watch.reset();
        info->freeGlobalLookup();
 
        delete std::get<0>(graphs);
        productBuilder.calculate(matrix->getmat(), *aggregatesMap, *coarseMatrix, *infoLevel, OverlapFlags());
 
        if(criterion.debugLevel()>2) {
          if(rank==0)
            std::cout<<"Calculation entries of Galerkin product took "<<watch.elapsed()<<" seconds."<<std::endl;
        }
 
        BIGINT nonzeros = countNonZeros(*coarseMatrix);
        allnonzeros = allnonzeros + infoLevel->communicator().sum(nonzeros);
        MatrixArgs args(*coarseMatrix, *infoLevel);
 
        matrices_.addCoarser(args);
        redistributes_.push_back(RedistributeInfoType());
      } // end level loop
 
 
      infoLevel->freeGlobalLookup();
 
      built_=true;
      AggregatesMap* aggregatesMap=new AggregatesMap(0);
      aggregatesMaps_.push_back(aggregatesMap);
 
      if(criterion.debugLevel()>0) {
        if(level==criterion.maxLevel()) {
          BIGINT unknownsLevel = mlevel->getmat().N();
          unknownsLevel = infoLevel->communicator().sum(unknownsLevel);
          double dunknownsLevel = unknownsLevel.todouble();
          if(rank==0 && criterion.debugLevel()>1) {
            std::cout<<"Level "<<level<<" has "<<dunknownsLevel<<" unknowns, "<<dunknownsLevel/infoLevel->communicator().size()
                     <<" unknowns per proc (procs="<<infoLevel->communicator().size()<<")"<<std::endl;
          }
        }
      }
 
      if(criterion.accumulate() && !redistributes_.back().isSetup() &&
         infoLevel->communicator().size()>1) {
#if HAVE_MPI && !HAVE_PARMETIS
        if(criterion.accumulate()==successiveAccu &&
           infoLevel->communicator().rank()==0)
          std::cerr<<"Successive accumulation of data on coarse levels only works with ParMETIS installed."
                   <<"  Fell back to accumulation to one domain on coarsest level"<<std::endl;
#endif
 
        // accumulate to fewer processors
        Matrix* redistMat= new Matrix();
        ParallelInformation* redistComm=0;
        int nodomains = 1;
 
        repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
                                       redistComm, redistributes_.back(), nodomains,criterion);
        MatrixArgs args(*redistMat, *redistComm);
        BIGINT unknownsRedist = redistMat->N();
        unknownsRedist = infoLevel->communicator().sum(unknownsRedist);
 
        if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1) {
          double dunknownsRedist = unknownsRedist.todouble();
          std::cout<<"Level "<<level<<" redistributed has "<<dunknownsRedist<<" unknowns, "<<dunknownsRedist/redistComm->communicator().size()
                   <<" unknowns per proc (procs="<<redistComm->communicator().size()<<")"<<std::endl;
        }
        mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
        infoLevel.addRedistributed(redistComm);
        infoLevel->freeGlobalLookup();
      }
 
      int levels = matrices_.levels();
      maxlevels_ = parallelInformation_.finest()->communicator().max(levels);
      assert(matrices_.levels()==redistributes_.size());
      if(hasCoarsest() && rank==0 && criterion.debugLevel()>1)
        std::cout<<"operator complexity: "<<allnonzeros.todouble()/finenonzeros.todouble()<<std::endl;
 
    }
 
    template<class M, class IS, class A>
    const typename MatrixHierarchy<M,IS,A>::ParallelMatrixHierarchy&
    MatrixHierarchy<M,IS,A>::matrices() const
    {
      return matrices_;
    }
 
    template<class M, class IS, class A>
    const typename MatrixHierarchy<M,IS,A>::ParallelInformationHierarchy&
    MatrixHierarchy<M,IS,A>::parallelInformation() const
    {
      return parallelInformation_;
    }
 
    template<class M, class IS, class A>
    void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const
    {
      int levels=aggregatesMaps().size();
      int maxlevels=parallelInformation_.finest()->communicator().max(levels);
      std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
      // We need an auxiliary vector for the consecutive prolongation.
      std::vector<std::size_t> tmp;
      std::vector<std::size_t> *coarse, *fine;
 
      // make sure the allocated space suffices.
      tmp.reserve(size);
      data.reserve(size);
 
      // Correctly assign coarse and fine for the first prolongation such that
      // we end up in data in the end.
      if(levels%2==0) {
        coarse=&tmp;
        fine=&data;
      }else{
        coarse=&data;
        fine=&tmp;
      }
 
      // Number the unknowns on the coarsest level consecutively for each process.
      if(levels==maxlevels) {
        const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
        std::size_t m=0;
 
        for(typename AggregatesMap::const_iterator iter = map.begin(); iter != map.end(); ++iter)
          if(*iter< AggregatesMap::ISOLATED)
            m=std::max(*iter,m);
 
        coarse->resize(m+1);
        std::size_t i=0;
        srand((unsigned)std::clock());
        std::set<size_t> used;
        for(typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
            ++iter, ++i)
        {
          std::pair<std::set<std::size_t>::iterator,bool> ibpair
            = used.insert(static_cast<std::size_t>((((double)rand())/(RAND_MAX+1.0)))*coarse->size());
 
          while(!ibpair.second)
            ibpair = used.insert(static_cast<std::size_t>((((double)rand())/(RAND_MAX+1.0))*coarse->size()));
          *iter=*(ibpair.first);
        }
      }
 
      typename ParallelInformationHierarchy::Iterator pinfo = parallelInformation().coarsest();
      --pinfo;
 
      // Now consecutively project the numbers to the finest level.
      for(typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
          aggregates != aggregatesMaps().rend(); ++aggregates,--levels) {
 
        fine->resize((*aggregates)->noVertices());
        fine->assign(fine->size(), 0);
        Transfer<typename AggregatesMap::AggregateDescriptor, std::vector<std::size_t>, ParallelInformation>
        ::prolongateVector(*(*aggregates), *coarse, *fine, static_cast<std::size_t>(1), *pinfo);
        --pinfo;
        std::swap(coarse, fine);
      }
 
      // Assertion to check that we really projected to data on the last step.
      assert(coarse==&data);
    }
 
    template<class M, class IS, class A>
    const typename MatrixHierarchy<M,IS,A>::AggregatesMapList&
    MatrixHierarchy<M,IS,A>::aggregatesMaps() const
    {
      return aggregatesMaps_;
    }
    template<class M, class IS, class A>
    const typename MatrixHierarchy<M,IS,A>::RedistributeInfoList&
    MatrixHierarchy<M,IS,A>::redistributeInformation() const
    {
      return redistributes_;
    }
 
    template<class M, class IS, class A>
    MatrixHierarchy<M,IS,A>::~MatrixHierarchy()
    {
      typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
      typedef typename ParallelMatrixHierarchy::Iterator Iterator;
      typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
 
      AggregatesMapIterator amap = aggregatesMaps_.rbegin();
      InfoIterator info = parallelInformation_.coarsest();
      for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest;  --level, --info, ++amap) {
        (*amap)->free();
        delete *amap;
        delete &level->getmat();
        if(level.isRedistributed())
          delete &(level.getRedistributed().getmat());
      }
      delete *amap;
    }
 
    template<class M, class IS, class A>
    template<class V, class TA>
    void MatrixHierarchy<M,IS,A>::coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const
    {
      assert(hierarchy.levels()==1);
      typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
      typedef typename RedistributeInfoList::const_iterator RIter;
      RIter redist = redistributes_.begin();
 
      Iterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
      int level=0;
      if(redist->isSetup())
        hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
      Dune::dvverb<<"Level "<<level<<" has "<<matrices_.finest()->getmat().N()<<" unknowns!"<<std::endl;
 
      while(matrix != coarsest) {
        ++matrix; ++level; ++redist;
        Dune::dvverb<<"Level "<<level<<" has "<<matrix->getmat().N()<<" unknowns!"<<std::endl;
 
        hierarchy.addCoarser(matrix->getmat().N());
        if(redist->isSetup())
          hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
 
      }
 
    }
 
    template<class M, class IS, class A>
    template<class S, class TA>
    void MatrixHierarchy<M,IS,A>::coarsenSmoother(Hierarchy<S,TA>& smoothers,
                                                  const typename SmootherTraits<S>::Arguments& sargs) const
    {
      assert(smoothers.levels()==0);
      typedef typename ParallelMatrixHierarchy::ConstIterator MatrixIterator;
      typedef typename ParallelInformationHierarchy::ConstIterator PinfoIterator;
      typedef typename AggregatesMapList::const_iterator AggregatesIterator;
 
      typename ConstructionTraits<S>::Arguments cargs;
      cargs.setArgs(sargs);
      PinfoIterator pinfo = parallelInformation_.finest();
      AggregatesIterator aggregates = aggregatesMaps_.begin();
      int level=0;
      for(MatrixIterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
          matrix != coarsest; ++matrix, ++pinfo, ++aggregates, ++level) {
        cargs.setMatrix(matrix->getmat(), **aggregates);
        cargs.setComm(*pinfo);
        smoothers.addCoarser(cargs);
      }
      if(maxlevels()>levels()) {
        // This is not the globally coarsest level and therefore smoothing is needed
        cargs.setMatrix(matrices_.coarsest()->getmat(), **aggregates);
        cargs.setComm(*pinfo);
        smoothers.addCoarser(cargs);
        ++level;
      }
    }
 
    template<class M, class IS, class A>
    template<class F>
    void MatrixHierarchy<M,IS,A>::recalculateGalerkin(const F& copyFlags)
    {
      typedef typename AggregatesMapList::iterator AggregatesMapIterator;
      typedef typename ParallelMatrixHierarchy::Iterator Iterator;
      typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
 
      AggregatesMapIterator amap = aggregatesMaps_.begin();
      BaseGalerkinProduct productBuilder;
      InfoIterator info = parallelInformation_.finest();
      typename RedistributeInfoList::iterator riIter = redistributes_.begin();
      Iterator level = matrices_.finest(), coarsest=matrices_.coarsest();
      if(level.isRedistributed()) {
        info->buildGlobalLookup(level->getmat().N());
        redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()),
                                  const_cast<Matrix&>(level.getRedistributed().getmat()),
                                  *info,info.getRedistributed(), *riIter);
        info->freeGlobalLookup();
      }
 
      for(; level!=coarsest; ++amap) {
        const Matrix& fine = (level.isRedistributed() ? level.getRedistributed() : *level).getmat();
        ++level;
        ++info;
        ++riIter;
        productBuilder.calculate(fine, *(*amap), const_cast<Matrix&>(level->getmat()), *info, copyFlags);
        if(level.isRedistributed()) {
          info->buildGlobalLookup(level->getmat().N());
          redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()),
                                    const_cast<Matrix&>(level.getRedistributed().getmat()), *info,
                                    info.getRedistributed(), *riIter);
          info->freeGlobalLookup();
        }
      }
    }
 
    template<class M, class IS, class A>
    std::size_t MatrixHierarchy<M,IS,A>::levels() const
    {
      return matrices_.levels();
    }
 
    template<class M, class IS, class A>
    std::size_t MatrixHierarchy<M,IS,A>::maxlevels() const
    {
      return maxlevels_;
    }
 
    template<class M, class IS, class A>
    bool MatrixHierarchy<M,IS,A>::hasCoarsest() const
    {
      return levels()==maxlevels() &&
             (!matrices_.coarsest().isRedistributed() ||matrices_.coarsest()->getmat().N()>0);
    }
 
    template<class M, class IS, class A>
    bool MatrixHierarchy<M,IS,A>::isBuilt() const
    {
      return built_;
    }
 
    template<class T, class A>
    Hierarchy<T,A>::Hierarchy()
      : finest_(0), coarsest_(0), nonAllocated_(0), allocator_(), levels_(0)
    {}
 
    template<class T, class A>
    Hierarchy<T,A>::Hierarchy(MemberType& first)
      : allocator_()
    {
      finest_ = allocator_.allocate(1,0);
      finest_->element_ = &first;
      finest_->redistributed_ = nullptr;
      nonAllocated_ = finest_;
      coarsest_ = finest_;
      coarsest_->coarser_ = coarsest_->finer_ = nullptr;
      levels_ = 1;
    }
 
    template<class T, class A>
    Hierarchy<T,A>::Hierarchy(MemberType* first)
      : allocator_()
    {
      finest_ = allocator_.allocate(1,0);
      finest_->element_ = first;
      finest_->redistributed_ = nullptr;
      nonAllocated_ = nullptr;
      coarsest_ = finest_;
      coarsest_->coarser_ = coarsest_->finer_ = nullptr;
      levels_ = 1;
    }
    template<class T, class A>
    Hierarchy<T,A>::~Hierarchy()
    {
      while(coarsest_) {
        Element* current = coarsest_;
        coarsest_ = coarsest_->finer_;
        if(current != nonAllocated_) {
          if(current->redistributed_)
            ConstructionTraits<T>::deconstruct(current->redistributed_);
          ConstructionTraits<T>::deconstruct(current->element_);
        }
        allocator_.deallocate(current, 1);
        current=nullptr;
        //coarsest_->coarser_ = nullptr;
      }
    }
 
    template<class T, class A>
    Hierarchy<T,A>::Hierarchy(const Hierarchy& other)
    : nonAllocated_(), allocator_(other.allocator_),
      levels_(other.levels_)
    {
      if(!other.finest_)
      {
        finest_=coarsest_=nonAllocated_=nullptr;
        return;
      }
      finest_=allocator_.allocate(1,0);
      Element* finer_         = nullptr;
      Element* current_      = finest_;
      Element* otherCurrent_ = other.finest_;
 
      while(otherCurrent_)
      {
        T* t=new T(*(otherCurrent_->element_));
        current_->element_=t;
        current_->finer_=finer_;
        if(otherCurrent_->redistributed_)
          current_->redistributed_ = new T(*otherCurrent_->redistributed_);
        else
          current_->redistributed_= nullptr;
        finer_=current_;
        if(otherCurrent_->coarser_)
        {
          current_->coarser_=allocator_.allocate(1,0);
          current_=current_->coarser_;
        }else
          current_->coarser_=nullptr;
        otherCurrent_=otherCurrent_->coarser_;
      }
      coarsest_=current_;
    }
 
    template<class T, class A>
    std::size_t Hierarchy<T,A>::levels() const
    {
      return levels_;
    }
 
    template<class T, class A>
    void Hierarchy<T,A>::addRedistributedOnCoarsest(Arguments& args)
    {
      coarsest_->redistributed_ = ConstructionTraits<MemberType>::construct(args);
    }
 
    template<class T, class A>
    void Hierarchy<T,A>::addCoarser(Arguments& args)
    {
      if(!coarsest_) {
        assert(!finest_);
        coarsest_ = allocator_.allocate(1,0);
        coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
        finest_ = coarsest_;
        coarsest_->finer_ = nullptr;
      }else{
        coarsest_->coarser_ = allocator_.allocate(1,0);
        coarsest_->coarser_->finer_ = coarsest_;
        coarsest_ = coarsest_->coarser_;
        coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
      }
      coarsest_->redistributed_ = nullptr;
      coarsest_->coarser_=nullptr;
      ++levels_;
    }
 
 
    template<class T, class A>
    void Hierarchy<T,A>::addFiner(Arguments& args)
    {
      if(!finest_) {
        assert(!coarsest_);
        finest_ = allocator_.allocate(1,0);
        finest_->element = ConstructionTraits<T>::construct(args);
        coarsest_ = finest_;
        coarsest_->coarser_ = coarsest_->finer_ = nullptr;
      }else{
        finest_->finer_ = allocator_.allocate(1,0);
        finest_->finer_->coarser_ = finest_;
        finest_ = finest_->finer_;
        finest_->finer = nullptr;
        finest_->element = ConstructionTraits<T>::construct(args);
      }
      ++levels_;
    }
 
    template<class T, class A>
    typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::finest()
    {
      return Iterator(finest_);
    }
 
    template<class T, class A>
    typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::coarsest()
    {
      return Iterator(coarsest_);
    }
 
    template<class T, class A>
    typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::finest() const
    {
      return ConstIterator(finest_);
    }
 
    template<class T, class A>
    typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::coarsest() const
    {
      return ConstIterator(coarsest_);
    }
  } // namespace Amg
} // namespace Dune
 
#endif