gate.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/base_header.hpp>
#include <kernel/util/dist.hpp>
#include <kernel/util/exception.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/global/synch_vec.hpp>
#include <kernel/global/synch_scal.hpp>
 
#include <vector>
 
namespace FEAT
{
  namespace Global
  {
    template<typename LocalVector_, typename Mirror_>
    class Gate
    {
    public:
      typedef typename LocalVector_::DataType DataType;
      typedef typename LocalVector_::IndexType IndexType;
      typedef LocalVector_ LocalVectorType;
      typedef Mirror_ MirrorType;
      typedef LAFEM::DenseVector<DataType, IndexType> BufferVectorType;
 
      typedef SynchScalarTicket<DataType> ScalarTicketType;
      typedef SynchVectorTicket<LocalVector_, Mirror_> VectorTicketType;
 
    public:
      const Dist::Comm* _comm;
      std::vector<int> _ranks;
      std::vector<Mirror_> _mirrors;
      LocalVector_ _freqs;
 
      template <typename LocalVector2_, typename Mirror2_>
      using GateType = Gate<LocalVector2_, Mirror2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using GateTypeByDI = Gate<typename LocalVector_::template ContainerType<DataType2_, IndexType2_>, typename Mirror_::template MirrorType<DataType2_, IndexType2_> >;
 
    public:
      explicit Gate() :
        _comm(nullptr)
      {
      }
 
      explicit Gate(const Dist::Comm& comm) :
        _comm(&comm)
      {
      }
 
      Gate(Gate&& other) :
        _comm(other._comm),
        _ranks(std::forward<std::vector<int>>(other._ranks)),
        _mirrors(std::forward<std::vector<Mirror_>>(other._mirrors)),
        _freqs(std::forward<LocalVector_>(other._freqs))
      {
      }
 
      Gate& operator=(Gate&& other)
      {
        if(this == &other)
        {
          return *this;
        }
 
        _comm = other._comm;
        _ranks = std::forward<std::vector<int>>(other._ranks);
        _mirrors = std::forward<std::vector<Mirror_>>(other._mirrors);
        _freqs = std::forward<LocalVector_>(other._freqs);
 
        return *this;
      }
 
      virtual ~Gate()
      {
      }
 
      const Dist::Comm* get_comm() const
      {
        return this->_comm;
      }
 
      void set_comm(const Dist::Comm* comm_)
      {
        this->_comm = comm_;
      }
 
      const std::vector<int> get_ranks() const
      {
        return this->_ranks;
      }
 
      const std::vector<Mirror_>& get_mirrors() const
      {
        return this->_mirrors;
      }
 
      const LocalVector_& get_freqs() const
      {
        return this->_freqs;
      }
 
      template<typename LVT2_, typename MT2_>
      void convert(const Gate<LVT2_, MT2_>& other)
      {
        if((void*)this == (void*)&other)
          return;
 
        this->_comm = other._comm;
        this->_ranks = other._ranks;
 
        this->_mirrors.resize(other._mirrors.size());
        for(std::size_t i(0); i < other._mirrors.size(); ++i)
        {
          this->_mirrors.at(i).convert(other._mirrors.at(i));
        }
 
        this->_freqs.convert(other._freqs);
      }
 
      template<typename LVT2_>
      void convert(const Gate<LVT2_, Mirror_>& other, LocalVector_&& vector, LAFEM::CloneMode mode = LAFEM::CloneMode::Shallow)
      {
        if((void*)this == (void*)&other)
          return;
 
        this->_ranks.clear();
        this->_mirrors.clear();
        this->_mirrors.resize(other._mirrors.size());
        this->_freqs.clear(); // will be rebuild by compile function
 
        this->_comm = other._comm;
        this->_ranks = other._ranks;
 
        // shallow-clone mirrors
        for(std::size_t i(0); i < _mirrors.size(); ++i)
        {
          _mirrors.at(i).clone(other._mirrors.at(i), mode);
        }
 
        // compile this gate
        compile(std::forward<LocalVector_>(vector));
      }
 
      std::size_t bytes() const
      {
        size_t temp(0);
        for (auto& i : _mirrors)
        {
          temp += i.bytes();
        }
        temp += _freqs.bytes();
        temp += _ranks.size() * sizeof(int);
 
        return temp;
      }
 
      void push(int rank, Mirror_&& mirror)
      {
        XASSERT(this->_comm != nullptr);
        XASSERT(rank < this->_comm->size());
 
        // push rank and tags
        _ranks.push_back(rank);
 
        // push mirror
        _mirrors.push_back(std::move(mirror));
      }
 
      void compile(LocalVector_&& vector)
      {
        // initialize frequency vector
        _freqs = std::move(vector);
        _freqs.format(DataType(1));
 
        // loop over all mirrors
        for(std::size_t i(0); i < _mirrors.size(); ++i)
        {
          // sum up number of ranks per frequency entry, listed in different mirrors
          auto temp = _mirrors.at(i).create_buffer(_freqs);
          temp.format(DataType(1));
 
          // gather-axpy into frequency vector
          _mirrors.at(i).scatter_axpy(_freqs, temp);
        }
 
        // invert frequencies
        _freqs.component_invert(_freqs);
      }
 
      template<LAFEM::Perspective perspective_ = LAFEM::Perspective::native>
      Index get_num_local_dofs() const
      {
        return _freqs.template size<perspective_>();
      }
 
      template<LAFEM::Perspective perspective_ = LAFEM::Perspective::native>
      Index get_num_global_dofs() const
      {
        XASSERT(this->_comm != nullptr);
        if(this->_comm->size() <= 1)
          return this->template get_num_local_dofs<perspective_>();
 
        // get my rank
        const int my_rank = this->_comm->rank();
 
        // get local number of DOFs
        const Index loc_dofs = this->template get_num_local_dofs<perspective_>();
 
        // create a local mask vector and format it to 1
        std::vector<int> mask(std::size_t(loc_dofs), 1);
 
        // set all DOFs, which are shared with a lower rank neighbor, to 0
        for(std::size_t i(0); i < _mirrors.size(); ++i)
        {
          if(this->_ranks.at(i) < my_rank)
            this->_mirrors.at(i).template mask_scatter<perspective_>(this->_freqs, mask, 0);
        }
 
        // count the number of DOFs that are still 1 and thus not owned by a lower rank neighbor
        Index owned_dofs(0u);
        for(const auto& k : mask)
          owned_dofs += Index(k);
 
        // now sum up the number of owned DOFs over all processes
        Index global_dofs(0u);
        this->_comm->allreduce(&owned_dofs, &global_dofs, std::size_t(1), Dist::op_sum);
 
        // done!
        return global_dofs;
      }
 
      void from_1_to_0(LocalVector_& vector) const
      {
        if(!_ranks.empty())
        {
          vector.component_product(vector, _freqs);
        }
      }
 
      void sync_0(LocalVector_& vector) const
      {
        if(!_ranks.empty())
        {
          SynchVectorTicket<LocalVector_, Mirror_> ticket(vector, *_comm, _ranks, _mirrors);
          ticket.wait();
        }
      }
 
      VectorTicketType sync_0_async(LocalVector_& vector) const
      {
        return SynchVectorTicket<LocalVector_, Mirror_>(vector, *_comm, _ranks, _mirrors);
      }
 
      void sync_1(LocalVector_& vector) const
      {
        if(!_ranks.empty())
        {
          from_1_to_0(vector);
          SynchVectorTicket<LocalVector_, Mirror_> ticket(vector, *_comm, _ranks, _mirrors);
          ticket.wait();
        }
      }
 
      VectorTicketType sync_1_async(LocalVector_& vector) const
      {
        from_1_to_0(vector);
        return SynchVectorTicket<LocalVector_, Mirror_>(vector, *_comm, _ranks, _mirrors);
      }
 
      DataType dot(const LocalVector_& x, const LocalVector_& y) const
      {
        // This is if there is only one process
        if(_comm == nullptr || _comm->size() == 1)
        {
          return x.dot(y);
        }
        // Even if there are no neighbors, we still need to sum up globally
        else if(_ranks.empty())
        {
          return sum(x.dot(y));
        }
        // If there are neighbors, we have to use the frequencies and sum up globally
        else
        {
          return sum(_freqs.triple_dot(x, y));
        }
      }
 
      ScalarTicketType dot_async(const LocalVector_& x, const LocalVector_& y, bool sqrt = false) const
      {
        return sum_async(_freqs.triple_dot(x, y), sqrt);
      }
 
      DataType sum(DataType x) const
      {
        return sum_async(x).wait();
      }
 
      ScalarTicketType sum_async(DataType x, bool sqrt = false) const
      {
        return SynchScalarTicket<DataType>(x, *_comm, Dist::op_sum, sqrt);
      }
 
      DataType min(DataType x) const
      {
        return min_async(x).wait();
      }
 
      ScalarTicketType min_async(DataType x) const
      {
        return SynchScalarTicket<DataType>(x, *_comm, Dist::op_min, false);
      }
 
      DataType max(DataType x) const
      {
        return max_async(x).wait();
      }
 
      ScalarTicketType max_async(DataType x) const
      {
        return SynchScalarTicket<DataType>(x, *_comm, Dist::op_max, false);
      }
 
      DataType norm2(DataType x) const
      {
        return norm2_async(x).wait();
      }
 
      ScalarTicketType norm2_async(DataType x) const
      {
        return SynchScalarTicket<DataType>(x*x, *_comm, Dist::op_sum, true);
      }
    }; // class Gate<...>
  } // namespace Global
} // namespace FEAT