aggregates.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_AGGREGATES_HH
#define DUNE_AMG_AGGREGATES_HH


#include "parameters.hh"
#include "graph.hh"
#include "properties.hh"
#include "combinedfunctor.hh"

#include <dune/common/timer.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/poolallocator.hh>
#include <dune/common/sllist.hh>
#include <dune/common/unused.hh>
#include <dune/common/ftraits.hh>

#include <utility>
#include <set>
#include <algorithm>
#include <complex>
#include <limits>
#include <ostream>
#include <tuple>

namespace Dune
{
  namespace Amg
  {

    template<class T>
    class AggregationCriterion : public T
    {

    public:
      typedef T DependencyPolicy;

      AggregationCriterion()
        : T()
      {}

      AggregationCriterion(const Parameters& parms)
        : T(parms)
      {}
      void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
      {
        this->setMaxDistance(diameter-1);
        std::size_t csize=1;

        for(; dim>0; dim--) {
          csize*=diameter;
          this->setMaxDistance(this->maxDistance()+diameter-1);
        }
        this->setMinAggregateSize(csize);
        this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5));
      }

      void setDefaultValuesAnisotropic(std::size_t dim,std::size_t diameter=2)
      {
        setDefaultValuesIsotropic(dim, diameter);
        this->setMaxDistance(this->maxDistance()+dim-1);
      }
    };

    template<class T>
    std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion)
    {
      os<<"{ maxdistance="<<criterion.maxDistance()<<" minAggregateSize="
      <<criterion.minAggregateSize()<< " maxAggregateSize="<<criterion.maxAggregateSize()
      <<" connectivity="<<criterion.maxConnectivity()<<" debugLevel="<<criterion.debugLevel()<<"}";
      return os;
    }

    template<class M, class N>
    class SymmetricMatrixDependency : public Dune::Amg::Parameters
    {
    public:
      typedef M Matrix;

      typedef N Norm;

      typedef typename Matrix::row_type Row;

      typedef typename Matrix::ConstColIterator ColIter;

      void init(const Matrix* matrix);

      void initRow(const Row& row, int index);

      void examine(const ColIter& col);

      template<class G>
      void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);

      bool isIsolated();


      SymmetricMatrixDependency(const Parameters& parms)
        : Parameters(parms)
      {}
      SymmetricMatrixDependency()
        : Parameters()
      {}

    protected:
      const Matrix* matrix_;
      typedef typename Matrix::field_type field_type;
      typedef typename FieldTraits<field_type>::real_type real_type;
      real_type maxValue_;
      Norm norm_;
      int row_;
      real_type diagonal_;
      std::vector<real_type> vals_;
      typename std::vector<real_type>::iterator valIter_;

    };


    template<class M, class N>
    inline void SymmetricMatrixDependency<M,N>::init(const Matrix* matrix)
    {
      matrix_ = matrix;
    }

    template<class M, class N>
    inline void SymmetricMatrixDependency<M,N>::initRow(const Row& row, int index)
    {
      vals_.assign(row.size(), 0.0);
      assert(vals_.size()==row.size());
      valIter_=vals_.begin();

      maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
      diagonal_=norm_(row[index]);
      row_ = index;
    }

    template<class M, class N>
    inline void SymmetricMatrixDependency<M,N>::examine(const ColIter& col)
    {
      // skip positive offdiagonals if norm preserves sign of them.
      real_type eij = norm_(*col);
      if(!N::is_sign_preserving || eij<0)  // || eji<0)
      {
        *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
        maxValue_ = std::max(maxValue_, *valIter_);
      }else
        *valIter_ =0;
      ++valIter_;
    }

    template<class M, class N>
    template<class G>
    inline void SymmetricMatrixDependency<M,N>::examine(G&, const typename G::EdgeIterator& edge, const ColIter&)
    {
      if(*valIter_ > alpha() * maxValue_) {
        edge.properties().setDepends();
        edge.properties().setInfluences();
      }
      ++valIter_;
    }

    template<class M, class N>
    inline bool SymmetricMatrixDependency<M,N>::isIsolated()
    {
      if(diagonal_==0)
        DUNE_THROW(Dune::ISTLError, "No diagonal entry for row "<<row_<<".");
      valIter_=vals_.begin();
      return maxValue_  < beta();
    }

    template<class M, class N>
    class Dependency : public Parameters
    {
    public:
      typedef M Matrix;

      typedef N Norm;

      typedef typename Matrix::row_type Row;

      typedef typename Matrix::ConstColIterator ColIter;

      void init(const Matrix* matrix);

      void initRow(const Row& row, int index);

      void examine(const ColIter& col);

      template<class G>
      void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);

      bool isIsolated();

      Dependency(const Parameters& parms)
        : Parameters(parms)
      {}

      Dependency()
        : Parameters()
      {}

    protected:
      const Matrix* matrix_;
      typedef typename Matrix::field_type field_type;
      typedef typename FieldTraits<field_type>::real_type real_type;
      real_type maxValue_;
      Norm norm_;
      int row_;
      real_type diagonal_;
    };

    template<class M, class N>
    class SymmetricDependency : public Parameters
    {
    public:
      typedef M Matrix;

      typedef N Norm;

      typedef typename Matrix::row_type Row;

      typedef typename Matrix::ConstColIterator ColIter;

      void init(const Matrix* matrix);

      void initRow(const Row& row, int index);

      void examine(const ColIter& col);

      template<class G>
      void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);

      bool isIsolated();


      SymmetricDependency(const Parameters& parms)
        : Parameters(parms)
      {}
      SymmetricDependency()
        : Parameters()
      {}

    protected:
      const Matrix* matrix_;
      typedef typename Matrix::field_type field_type;
      typedef typename FieldTraits<field_type>::real_type real_type;
      real_type maxValue_;
      Norm norm_;
      int row_;
      real_type diagonal_;
    };

    template<int N>
    class Diagonal
    {
    public:
      enum { /* @brief We preserve the sign.*/
        is_sign_preserving = true
      };

      template<class M>
      typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
      {
        typedef typename M::field_type field_type;
        typedef typename FieldTraits<field_type>::real_type real_type;
        static_assert( std::is_convertible<field_type, real_type >::value,
                  "use of diagonal norm in AMG not implemented for complex field_type");
        return m[N][N];
        // possible implementation for complex types: return signed_abs(m[N][N]);
      }

    private:

      template<typename T>
      static T signed_abs(const T & v)
      {
        return v;
      }

      template<typename T>
      static T signed_abs(const std::complex<T> & v)
      {
        // return sign * abs_value
        // in case of complex numbers this extends to using the csgn function to determine the sign
        return csgn(v) * std::abs(v);
      }

      template<typename T>
      static T csgn(const T & v)
      {
        return (T(0) < v) - (v < T(0));
      }

      template<typename T>
      static T csgn(std::complex<T> a)
      {
        return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
      }

    };

    class FirstDiagonal : public Diagonal<0>
    {};

    struct RowSum
    {

      enum { /* @brief We preserve the sign.*/
        is_sign_preserving = false
      };
      template<class M>
      typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
      {
        return m.infinity_norm();
      }
    };

    struct FrobeniusNorm
    {

      enum { /* @brief We preserve the sign.*/
        is_sign_preserving = false
      };
      template<class M>
      typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
      {
        return m.frobenius_norm();
      }
    };
    struct AlwaysOneNorm
    {

      enum { /* @brief We preserve the sign.*/
        is_sign_preserving = false
      };
      template<class M>
      typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
      {
        return 1;
      }
    };
    template<class M, class Norm>
    class SymmetricCriterion : public AggregationCriterion<SymmetricDependency<M,Norm> >
    {
    public:
      SymmetricCriterion(const Parameters& parms)
        : AggregationCriterion<SymmetricDependency<M,Norm> >(parms)
      {}
      SymmetricCriterion()
      {}
    };


    template<class M, class Norm>
    class UnSymmetricCriterion : public AggregationCriterion<Dependency<M,Norm> >
    {
    public:
      UnSymmetricCriterion(const Parameters& parms)
        : AggregationCriterion<Dependency<M,Norm> >(parms)
      {}
      UnSymmetricCriterion()
      {}
    };
    // forward declaration
    template<class G> class Aggregator;


    template<class V>
    class AggregatesMap
    {
    public:

      static const V UNAGGREGATED;

      static const V ISOLATED;
      typedef V VertexDescriptor;

      typedef V AggregateDescriptor;

      typedef PoolAllocator<VertexDescriptor,100> Allocator;

      typedef SLList<VertexDescriptor,Allocator> VertexList;

      class DummyEdgeVisitor
      {
      public:
        template<class EdgeIterator>
        void operator()(const EdgeIterator& edge) const
        {
          DUNE_UNUSED_PARAMETER(edge);
        }
      };


      AggregatesMap();

      AggregatesMap(std::size_t noVertices);

      ~AggregatesMap();

      template<class M, class G, class C>
      std::tuple<int,int,int,int> buildAggregates(const M& matrix, G& graph, const C& criterion,
                                                  bool finestLevel);

      template<bool reset, class G, class F, class VM>
      std::size_t breadthFirstSearch(const VertexDescriptor& start,
                                     const AggregateDescriptor& aggregate,
                                     const G& graph,
                                     F& aggregateVisitor,
                                     VM& visitedMap) const;

      template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
      std::size_t breadthFirstSearch(const VertexDescriptor& start,
                                     const AggregateDescriptor& aggregate,
                                     const G& graph, L& visited, F1& aggregateVisitor,
                                     F2& nonAggregateVisitor,
                                     VM& visitedMap) const;

      void allocate(std::size_t noVertices);

      std::size_t noVertices() const;

      void free();

      AggregateDescriptor& operator[](const VertexDescriptor& v);

      const AggregateDescriptor& operator[](const VertexDescriptor& v) const;

      typedef const AggregateDescriptor* const_iterator;

      const_iterator begin() const
      {
        return aggregates_;
      }

      const_iterator end() const
      {
        return aggregates_+noVertices();
      }

      typedef AggregateDescriptor* iterator;

      iterator begin()
      {
        return aggregates_;
      }

      iterator end()
      {
        return aggregates_+noVertices();
      }
    private:
      AggregatesMap(const AggregatesMap<V>&) = delete;
      AggregatesMap<V>& operator=(const AggregatesMap<V>&) = delete;

      AggregateDescriptor* aggregates_;

      std::size_t noVertices_;
    };

    template<class G, class C>
    void buildDependency(G& graph,
                         const typename C::Matrix& matrix,
                         C criterion,
                         bool finestLevel);

    template<class G, class S>
    class Aggregate
    {

    public:

      /***
       * @brief The type of the matrix graph we work with.
       */
      typedef G MatrixGraph;
      typedef typename MatrixGraph::VertexDescriptor Vertex;

      typedef PoolAllocator<Vertex,100> Allocator;

      typedef S VertexSet;

      typedef typename VertexSet::const_iterator const_iterator;

      typedef std::size_t* SphereMap;

      Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
                VertexSet& connectivity, std::vector<Vertex>& front_);

      void invalidate()
      {
        --id_;
      }

      void reconstruct(const Vertex& vertex);

      void seed(const Vertex& vertex);

      void add(const Vertex& vertex);

      void add(std::vector<Vertex>& vertex);
      void clear();

      typename VertexSet::size_type size();
      typename VertexSet::size_type connectSize();

      int id();

      const_iterator begin() const;

      const_iterator end() const;

    private:
      VertexSet vertices_;

      int id_;

      MatrixGraph& graph_;

      AggregatesMap<Vertex>& aggregates_;

      VertexSet& connected_;

      std::vector<Vertex>& front_;
    };

    template<class G>
    class Aggregator
    {
    public:

      typedef G MatrixGraph;

      typedef typename MatrixGraph::VertexDescriptor Vertex;

      typedef typename MatrixGraph::VertexDescriptor AggregateDescriptor;

      Aggregator();

      ~Aggregator();

      template<class M, class C>
      std::tuple<int,int,int,int> build(const M& m, G& graph,
                                        AggregatesMap<Vertex>& aggregates, const C& c,
                                        bool finestLevel);
    private:
      typedef PoolAllocator<Vertex,100> Allocator;

      typedef SLList<Vertex,Allocator> VertexList;

      typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;

      typedef std::size_t* SphereMap;

      MatrixGraph* graph_;

      Aggregate<MatrixGraph,VertexSet>* aggregate_;

      std::vector<Vertex> front_;

      VertexSet connected_;

      int size_;

      class Stack
      {
      public:
        static const Vertex NullEntry;

        Stack(const MatrixGraph& graph,
              const Aggregator<G>& aggregatesBuilder,
              const AggregatesMap<Vertex>& aggregates);
        ~Stack();
        Vertex pop();
      private:
        enum { N = 1300000 };

        const MatrixGraph& graph_;
        const Aggregator<G>& aggregatesBuilder_;
        const AggregatesMap<Vertex>& aggregates_;
        int size_;
        Vertex maxSize_;
        typename MatrixGraph::ConstVertexIterator begin_;
        typename MatrixGraph::ConstVertexIterator end_;

        Vertex* vals_;

      };

      friend class Stack;

      template<class V>
      void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate,
                                    const AggregatesMap<Vertex>& aggregates,
                                    V& visitor) const;

      template<class V>
      class AggregateVisitor
      {
      public:
        typedef V Visitor;
        AggregateVisitor(const AggregatesMap<Vertex>& aggregates, const AggregateDescriptor& aggregate,
                         Visitor& visitor);

        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);

      private:
        const AggregatesMap<Vertex>& aggregates_;
        AggregateDescriptor aggregate_;
        Visitor* visitor_;
      };

      class Counter
      {
      public:
        Counter();
        int value();

      protected:
        void increment();
        void decrement();

      private:
        int count_;
      };


      class FrontNeighbourCounter : public Counter
      {
      public:
        FrontNeighbourCounter(const MatrixGraph& front);

        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);

      private:
        const MatrixGraph& graph_;
      };

      int noFrontNeighbours(const Vertex& vertex) const;

      class TwoWayCounter : public Counter
      {
      public:
        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
      };

      int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
                            const AggregatesMap<Vertex>& aggregates) const;

      class OneWayCounter : public Counter
      {
      public:
        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
      };

      int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
                            const AggregatesMap<Vertex>& aggregates) const;

      class ConnectivityCounter : public Counter
      {
      public:
        ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates);

        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);

      private:
        const VertexSet& connected_;
        const AggregatesMap<Vertex>& aggregates_;

      };

      double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
      bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate,
                     const AggregatesMap<Vertex>& aggregates) const;

      bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList,
                     const AggregatesMap<Vertex>& aggregates) const;

      class DependencyCounter : public Counter
      {
      public:
        DependencyCounter();

        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
      };

      class FrontMarker
      {
      public:
        FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph);

        void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);

      private:
        std::vector<Vertex>& front_;
        MatrixGraph& graph_;
      };

      void unmarkFront();

      int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;

      std::pair<int,int> neighbours(const Vertex& vertex,
                                    const AggregateDescriptor& aggregate,
                                    const AggregatesMap<Vertex>& aggregates) const;
      int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;

      bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;

      std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates);

      Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;

      void nonisoNeighbourAggregate(const Vertex& vertex,
                                    const AggregatesMap<Vertex>& aggregates,
                                    SLList<Vertex>& neighbours) const;

      template<class C>
      void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
      template<class C>
      void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
    };

#ifndef DOXYGEN

    template<class M, class N>
    inline void SymmetricDependency<M,N>::init(const Matrix* matrix)
    {
      matrix_ = matrix;
    }

    template<class M, class N>
    inline void SymmetricDependency<M,N>::initRow(const Row& row, int index)
    {
      DUNE_UNUSED_PARAMETER(row);
      maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
      row_ = index;
      diagonal_ = norm_(matrix_->operator[](row_)[row_]);
    }

    template<class M, class N>
    inline void SymmetricDependency<M,N>::examine(const ColIter& col)
    {
      real_type eij = norm_(*col);
      typename Matrix::ConstColIterator opposite_entry =
        matrix_->operator[](col.index()).find(row_);
      if ( opposite_entry == matrix_->operator[](col.index()).end() )
      {
        // Consider this a weak connection we disregard.
        return;
      }
      real_type eji = norm_(*opposite_entry);

      // skip positive offdiagonals if norm preserves sign of them.
      if(!N::is_sign_preserving || eij<0 || eji<0)
        maxValue_ = std::max(maxValue_,
                             eij /diagonal_ * eji/
                             norm_(matrix_->operator[](col.index())[col.index()]));
    }

    template<class M, class N>
    template<class G>
    inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
    {
      real_type eij = norm_(*col);
      typename Matrix::ConstColIterator opposite_entry =
        matrix_->operator[](col.index()).find(row_);

      if ( opposite_entry == matrix_->operator[](col.index()).end() )
      {
        // Consider this as a weak connection we disregard.
        return;
      }
      real_type eji = norm_(*opposite_entry);
      // skip positve offdiagonals if norm preserves sign of them.
      if(!N::is_sign_preserving || (eij<0 || eji<0))
        if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
           eij/ diagonal_ > alpha() * maxValue_) {
          edge.properties().setDepends();
          edge.properties().setInfluences();
          typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
          other.setInfluences();
          other.setDepends();
        }
    }

    template<class M, class N>
    inline bool SymmetricDependency<M,N>::isIsolated()
    {
      return maxValue_  < beta();
    }


    template<class M, class N>
    inline void Dependency<M,N>::init(const Matrix* matrix)
    {
      matrix_ = matrix;
    }

    template<class M, class N>
    inline void Dependency<M,N>::initRow(const Row& row, int index)
    {
      DUNE_UNUSED_PARAMETER(row);
      maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
      row_ = index;
      diagonal_ = norm_(matrix_->operator[](row_)[row_]);
    }

    template<class M, class N>
    inline void Dependency<M,N>::examine(const ColIter& col)
    {
      maxValue_ = std::max(maxValue_,
                           -norm_(*col));
    }

    template<class M, class N>
    template<class G>
    inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
    {
      if(-norm_(*col) >= maxValue_ * alpha()) {
        edge.properties().setDepends();
        typedef typename G::EdgeDescriptor ED;
        ED e= graph.findEdge(edge.target(), edge.source());
        if(e!=std::numeric_limits<ED>::max())
        {
          typename G::EdgeProperties& other = graph.getEdgeProperties(e);
          other.setInfluences();
        }
      }
    }

    template<class M, class N>
    inline bool Dependency<M,N>::isIsolated()
    {
      return maxValue_  < beta() * diagonal_;
    }

    template<class G,class S>
    Aggregate<G,S>::Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
                              VertexSet& connected, std::vector<Vertex>& front)
      : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
        connected_(connected), front_(front)
    {}

    template<class G,class S>
    void Aggregate<G,S>::reconstruct(const Vertex& vertex)
    {
      /*
         vertices_.push_back(vertex);
         typedef typename VertexList::const_iterator iterator;
         iterator begin = vertices_.begin();
         iterator end   = vertices_.end();*/
      throw "Not yet implemented";

      //      while(begin!=end){
      //for();
      //      }

    }

    template<class G,class S>
    inline void Aggregate<G,S>::seed(const Vertex& vertex)
    {
      dvverb<<"Connected cleared"<<std::endl;
      connected_.clear();
      vertices_.clear();
      connected_.insert(vertex);
      dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
      ++id_ ;
      add(vertex);
    }


    template<class G,class S>
    inline void Aggregate<G,S>::add(const Vertex& vertex)
    {
      vertices_.insert(vertex);
      aggregates_[vertex]=id_;
      if(front_.size())
        front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));


      typedef typename MatrixGraph::ConstEdgeIterator iterator;
      const iterator end = graph_.endEdges(vertex);
      for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
        dvverb << " Inserting "<<aggregates_[edge.target()];
        connected_.insert(aggregates_[edge.target()]);
        dvverb <<" size="<<connected_.size();
        if