indicescoarsener.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_AMG_INDICESCOARSENER_HH
#define DUNE_AMG_INDICESCOARSENER_HH

#include <dune/common/parallel/indicessyncer.hh>
#include <dune/common/unused.hh>
#include <vector>
#include "renumberer.hh"

#if HAVE_MPI
#include <dune/istl/owneroverlapcopy.hh>
#endif

#include "pinfo.hh"

namespace Dune
{
  namespace Amg
  {

    template<typename T, typename E>
    class IndicesCoarsener
    {};


#if HAVE_MPI

    template<typename T, typename E>
    class ParallelIndicesCoarsener
    {
    public:
      typedef E ExcludedAttributes;

      typedef T ParallelInformation;

      typedef typename ParallelInformation::ParallelIndexSet ParallelIndexSet;

      typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;

      typedef typename ParallelIndexSet::LocalIndex LocalIndex;

      typedef typename LocalIndex::Attribute Attribute;

      typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;

      template<typename Graph, typename VM>
      static typename Graph::VertexDescriptor
      coarsen(ParallelInformation& fineInfo,
              Graph& fineGraph,
              VM& visitedMap,
              AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
              ParallelInformation& coarseInfo,
              typename Graph::VertexDescriptor noAggregates);

    private:
      template<typename G, typename I>
      class ParallelAggregateRenumberer : public AggregateRenumberer<G>
      {
        typedef typename G::VertexDescriptor Vertex;

        typedef I GlobalLookupIndexSet;

        typedef typename GlobalLookupIndexSet::IndexPair IndexPair;

        typedef typename IndexPair::GlobalIndex GlobalIndex;

      public:
        ParallelAggregateRenumberer(AggregatesMap<Vertex>& aggregates, const I& lookup)
          :  AggregateRenumberer<G>(aggregates),  isPublic_(false), lookup_(lookup),
            globalIndex_(std::numeric_limits<GlobalIndex>::max())
        {}


        void operator()(const typename G::ConstEdgeIterator& edge)
        {
          AggregateRenumberer<G>::operator()(edge);
          const IndexPair* pair= lookup_.pair(edge.target());
          if(pair!=0) {
            globalIndex(pair->global());
            attribute(pair->local().attribute());
            isPublic(pair->local().isPublic());
          }
        }

        Vertex operator()(const GlobalIndex& global)
        {
          DUNE_UNUSED_PARAMETER(global);
          Vertex current = this->number_;
          this->operator++();
          return current;
        }

        bool isPublic()
        {
          return isPublic_;
        }

        void isPublic(bool b)
        {
          isPublic_ = isPublic_ || b;
        }

        void reset()
        {
          globalIndex_ = std::numeric_limits<GlobalIndex>::max();
          isPublic_=false;
        }

        void attribute(const Attribute& attribute)
        {
          attribute_=attribute;
        }

        Attribute attribute()
        {
          return attribute_;
        }

        const GlobalIndex& globalIndex() const
        {
          return globalIndex_;
        }

        void globalIndex(const GlobalIndex& global)
        {
          globalIndex_ = global;
        }

      private:
        bool isPublic_;
        Attribute attribute_;
        const GlobalLookupIndexSet& lookup_;
        GlobalIndex globalIndex_;
      };

      template<typename Graph, typename VM, typename I>
      static void buildCoarseIndexSet(const ParallelInformation& pinfo,
                                      Graph& fineGraph,
                                      VM& visitedMap,
                                      AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                      ParallelIndexSet& coarseIndices,
                                      ParallelAggregateRenumberer<Graph,I>& renumberer);

      template<typename Graph,typename I>
      static void buildCoarseRemoteIndices(const RemoteIndices& fineRemote,
                                           const AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                           ParallelIndexSet& coarseIndices,
                                           RemoteIndices& coarseRemote,
                                           ParallelAggregateRenumberer<Graph,I>& renumberer);

    };

    template<typename G, typename L, typename E>
    class IndicesCoarsener<OwnerOverlapCopyCommunication<G,L>,E>
      : public ParallelIndicesCoarsener<OwnerOverlapCopyCommunication<G,L>,E>
    {};


#endif

    template<typename E>
    class IndicesCoarsener<SequentialInformation,E>
    {
    public:
      template<typename Graph, typename VM>
      static typename Graph::VertexDescriptor
      coarsen(const SequentialInformation & fineInfo,
              Graph& fineGraph,
              VM& visitedMap,
              AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
              SequentialInformation& coarseInfo,
              typename Graph::VertexDescriptor noAggregates);
    };

#if HAVE_MPI
    template<typename T, typename E>
    template<typename Graph, typename VM>
    inline typename Graph::VertexDescriptor
    ParallelIndicesCoarsener<T,E>::coarsen(ParallelInformation& fineInfo,
                                           Graph& fineGraph,
                                           VM& visitedMap,
                                           AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                           ParallelInformation& coarseInfo,
                                           typename Graph::VertexDescriptor noAggregates)
    {
      DUNE_UNUSED_PARAMETER(noAggregates);
      ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup());
      buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates,
                          coarseInfo.indexSet(), renumberer);
      buildCoarseRemoteIndices(fineInfo.remoteIndices(), aggregates, coarseInfo.indexSet(),
                               coarseInfo.remoteIndices(), renumberer);

      return renumberer;
    }

    template<typename T, typename E>
    template<typename Graph, typename VM, typename I>
    void ParallelIndicesCoarsener<T,E>::buildCoarseIndexSet(const ParallelInformation& pinfo,
                                                            Graph& fineGraph,
                                                            VM& visitedMap,
                                                            AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                                            ParallelIndexSet& coarseIndices,
                                                            ParallelAggregateRenumberer<Graph,I>& renumberer)
    {
      // fineGraph is the local subgraph corresponding to the vertices the process owns.
      // i.e. no overlap/copy vertices can be visited traversing the graph
      typedef typename Graph::ConstVertexIterator Iterator;
      typedef typename ParallelInformation::GlobalLookupIndexSet GlobalLookupIndexSet;

      Iterator end = fineGraph.end();
      const GlobalLookupIndexSet& lookup = pinfo.globalLookup();

      coarseIndices.beginResize();

      // Setup the coarse index set and renumber the aggregate consecutively
      // ascending from zero according to the minimum global index belonging
      // to the aggregate
      for(Iterator index = fineGraph.begin(); index != end; ++index) {
        if(aggregates[*index]!=AggregatesMap<typename Graph::VertexDescriptor>::ISOLATED)
          // Isolated vertices will not be represented on the next level.
          // These should only be there if skipIsolated is activiated in
          // the coarsening criterion as otherwise they will be aggregated
          // and should have real aggregate number in the map right now.
          if(!get(visitedMap, *index)) {
            // This vertex was not visited by breadthFirstSearch yet.
            typedef typename GlobalLookupIndexSet::IndexPair IndexPair;
            const IndexPair* pair= lookup.pair(*index);

            renumberer.reset(); // reset attribute and global index.
            if(pair!=0) {
              // vertex is in the index set. Note that not all vertices have
              // to be in the index set, just the ones where communication
              // will happen.
              assert(!ExcludedAttributes::contains(pair->local().attribute()));
              renumberer.attribute(pair->local().attribute());
              renumberer.isPublic(pair->local().isPublic());
              renumberer.globalIndex(pair->global());
            }

            // Reconstruct aggregate and mark vertices as visited
            aggregates.template breadthFirstSearch<false>(*index, aggregates[*index],
                                                          fineGraph, renumberer, visitedMap);

            typedef typename GlobalLookupIndexSet::IndexPair::GlobalIndex GlobalIndex;

            if(renumberer.globalIndex()!=std::numeric_limits<GlobalIndex>::max()) {
              // vertex is in the index set.
              //std::cout <<" Adding global="<< renumberer.globalIndex()<<" local="<<static_cast<std::size_t>(renumberer)<<std::endl;
              coarseIndices.add(renumberer.globalIndex(),
                                LocalIndex(renumberer, renumberer.attribute(),
                                           renumberer.isPublic()));
            }

            aggregates[*index] = renumberer;
            ++renumberer;
          }
      }

      coarseIndices.endResize();

      assert(static_cast<std::size_t>(renumberer) >= coarseIndices.size());

      // Reset the visited flags
      for(Iterator vertex=fineGraph.begin(); vertex != end; ++vertex)
        put(visitedMap, *vertex, false);
    }

    template<typename T, typename E>
    template<typename Graph, typename I>
    void ParallelIndicesCoarsener<T,E>::buildCoarseRemoteIndices(