alg_dof_parti.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/adjacency/dynamic_graph.hpp>
#include <kernel/global/gate.hpp>
#include <kernel/global/vector.hpp>
#include <kernel/global/filter.hpp>
#include <kernel/global/matrix.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/unit_filter.hpp>
#include <kernel/lafem/vector_mirror.hpp>
#include <kernel/lafem/matrix_mirror.hpp>
 
#include <algorithm>
#include <vector>
 
namespace FEAT
{
  namespace Global
  {
    template<typename LocalVector_, typename Mirror_>
    class AlgDofParti;
 
    template<typename DT_, typename IT_>
    class AlgDofParti<LAFEM::DenseVector<DT_, IT_>, LAFEM::VectorMirror<DT_, IT_>>
    {
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
 
      typedef LAFEM::DenseVector<DataType, IndexType> LocalVectorType;
      typedef LAFEM::VectorMirror<DataType, IndexType> MirrorType;
      typedef Global::Vector<LocalVectorType, MirrorType> GlobalVectorType;
      typedef Global::Gate<LocalVectorType, MirrorType> GateType;
 
    protected:
      const Dist::Comm* _comm;
      IndexType _glob_dof_offset;
      IndexType _glob_dof_count;
      std::vector<IndexType> _glob_dof_idx;
      MirrorType _owned_mirror;
      std::vector<std::pair<int, MirrorType>> _neighbors_owner;
      std::vector<std::pair<int, MirrorType>> _neighbors_donee;
 
      std::vector<int> _all_glob_dof_offset;
      std::vector<int> _all_glob_dof_counts;
 
    public:
      AlgDofParti() :
        _comm(nullptr),
        _glob_dof_offset(0),
        _glob_dof_count(0)
      {
      }
 
      void clear()
      {
        _comm = nullptr;
        _glob_dof_offset = IndexType(0);
        _glob_dof_count = IndexType(0);
        _glob_dof_idx.clear();
        _owned_mirror.clear();
        _neighbors_owner.clear();
        _neighbors_donee.clear();
        _all_glob_dof_offset.clear();
        _all_glob_dof_counts.clear();
      }
 
      const Dist::Comm* get_comm() const
      {
        return _comm;
      }
 
      IndexType get_num_owned_dofs() const
      {
        return IndexType(this->_owned_mirror.num_indices());
      }
 
      IndexType get_num_local_dofs() const
      {
        return IndexType(this->_owned_mirror.size());
      }
 
      IndexType get_num_global_dofs() const
      {
        return this->_glob_dof_count;
      }
 
      IndexType get_global_dof_offset() const
      {
        return this->_glob_dof_offset;
      }
 
      const std::vector<IndexType> get_global_dof_indices() const
      {
        return this->_glob_dof_idx;
      }
 
      IndexType map_local_to_global_index(const IndexType local_dof_idx) const
      {
        ASSERT(local_dof_idx < this->_owned_mirror.size());
        return this->_glob_dof_idx.at(local_dof_idx);
      }
 
      IndexType map_owned_to_local_index(const IndexType owned_dof_idx) const
      {
        ASSERT(owned_dof_idx < this->_owned_mirror.num_indices());
        return this->_owned_mirror.indices()[owned_dof_idx];
      }
 
      const MirrorType& get_owned_mirror() const
      {
        return this->_owned_mirror;
      }
 
      IndexType get_num_owner_neighbors() const
      {
        return IndexType(this->_neighbors_owner.size());
      }
 
      IndexType get_num_donee_neighbors() const
      {
        return IndexType(this->_neighbors_donee.size());
      }
 
      int get_owner_rank(IndexType i) const
      {
        return this->_neighbors_owner.at(i).first;
      }
 
      int get_donee_rank(IndexType i) const
      {
        return this->_neighbors_donee.at(i).first;
      }
 
      const MirrorType& get_owner_mirror(IndexType i) const
      {
        return this->_neighbors_owner.at(i).second;
      }
 
      const MirrorType& get_donee_mirror(IndexType i) const
      {
        return this->_neighbors_donee.at(i).second;
      }
 
      const std::vector<int>& get_all_global_dof_offsets() const
      {
        XASSERTM(!this->_all_glob_dof_offset.empty(), "You did not ask to assemble the global dof offsets");
        return this->_all_glob_dof_offset;
      }
 
      const std::vector<int>& get_all_global_dof_counts() const
      {
        XASSERTM(!this->_all_glob_dof_counts.empty(), "You did not ask to assemble the global dof counts");
        return this->_all_glob_dof_counts;
      }
 
      void assemble_allgather(bool yes_i_really_want_this = false)
      {
        XASSERTM(yes_i_really_want_this, "You probably don't want to do this!");
 
        std::size_t num_ranks = std::size_t(this->_comm->size());
 
        int my_nod = int(this->get_num_owned_dofs());
        int my_off = int(this->get_global_dof_offset());
 
        // resize vectors
        this->_all_glob_dof_offset.resize(num_ranks, 0);
        this->_all_glob_dof_counts.resize(num_ranks, 0);
 
        // allgather offsets + counts
        this->_comm->allgather(&my_off, 1, this->_all_glob_dof_offset.data(), 1);
        this->_comm->allgather(&my_nod, 1, this->_all_glob_dof_counts.data(), 1);
      }
 
      void assemble_by_gate(const GateType& gate)
      {
        // set our communicator
        this->_comm = gate.get_comm();
        XASSERTM(this->_comm != nullptr, "need a communicator here!");
 
        // get the communicator
        const Dist::Comm& comm = *this->_comm;
 
        // get my rank and neighbor ranks
        const int my_rank = comm.rank();
        const std::vector<int>& neighbor_ranks = gate._ranks;
 
        // get number of local DOFs on our patch
        const IndexType num_local_dofs = IndexType(gate._freqs.size());
 
        // As a very first step, we have to decide which process will become the
        // owner of each DOF. In this implementation, each DOF is owned by by the
        // process with the lowest rank, which simplifies a lot of things.
 
        // So, first create a vector, which stores the rank of the owner process
        // of each of our local DOFs and format the vector to our rank, i.e. at
        // startup assume that we own all local DOFs.
        std::vector<int> dof_owners(num_local_dofs, my_rank);
 
        // Now loop over all our neighbor processes
        for(std::size_t i(0); i < neighbor_ranks.size(); ++i)
        {
          // Check whether the neighbors rank is less than our rank, as otherwise
          // that particular neighbor can not own any of our local DOFs.
          const int neighbor_rank = neighbor_ranks.at(i);
          if(neighbor_rank < my_rank)
          {
            // loop over all DOFs in the mirror of that neighbor
            const MirrorType& mirror = gate._mirrors.at(i);
            const Index num_indices = mirror.num_indices();
            const IndexType* idx = mirror.indices();
            for(Index j(0); j < num_indices; ++j)
            {
              // get a reference to the current owner rank of that DOF and
              // update the DOF owner if that neighbor has a lower rank
              // than the previous "owner candidate"
              int& dof_owner = dof_owners.at(idx[j]);
              if(neighbor_rank < dof_owner)
                dof_owner = neighbor_rank;
            }
          }
        }
 
        // Okay, at this point we (and all other processes) know the owners of
        // each of our local DOFs. Now, we have to generate a vector of all
        // local DOF indices of the DOFs that we own:
        std::vector<IndexType> owned_dofs;
        owned_dofs.reserve(num_local_dofs);
 
        // Furthermore, we need the 'inverse' information, i.e. for each of our
        // local DOFs, we have to store the "owned DOF index" if we own this
        // particular DOF. We initialize the vector with ~0, which stands for
        // "not my DOF":
        std::vector<IndexType> own_dof_idx(num_local_dofs, ~IndexType(0));
 
        // Loop over all local DOFs
        for(IndexType i(0); i < num_local_dofs; ++i)
        {
          // Am I the owner of this DOF?
          if(dof_owners[i] == my_rank)
          {
            // Yes, that's my DOF, so let's give it a new 'my owned DOF' index:
            own_dof_idx.at(i) = IndexType(owned_dofs.size());
            // And remember that I own this DOF:
            owned_dofs.push_back(i);
          }
        }
 
        XASSERTM(!owned_dofs.empty(), "this process has no DOFs!");
 
        // get number of my owned DOFS
        const IndexType num_owned_dofs = IndexType(owned_dofs.size());
 
        // Next, we have to determine the global offset of our first owned DOF, which
        // is easily obtained by an exclusive-scan over each processes owned DOF count:
        comm.exscan(&num_owned_dofs, &this->_glob_dof_offset, 1, Dist::op_sum);
 
        // Furthermore, we also perform an allreduce to obtain the total number of
        // global DOFs, which is generally helpful and can be used for sanity checks.
        comm.allreduce(&num_owned_dofs, &this->_glob_dof_count, 1, Dist::op_sum);
 
        // Next, we have to set up the mirror for all of our owned DOFs, which we can
        // easily generate from our "owned_dofs" vector:
        this->_owned_mirror = _vidx2mirror(num_local_dofs, owned_dofs);
 
        // Now we also have to create the mirrors for each of our neighbor processes:
        for(std::size_t i(0); i < neighbor_ranks.size(); ++i)
        {
          // get the neighbor rank and its mirror
          const int neighbor_rank = neighbor_ranks.at(i);
          const MirrorType& mirror = gate._mirrors.at(i);
 
          // There are two cases now:
          // 1) The neighbor process has a lower rank, so it is a (potential) owner
          //    of at least one of our local DOFs.
          // 2) The neighbor process has a greater rank, so we are an owner of at least
          //    one of the neighbor process's local DOFs.
          if(neighbor_rank < my_rank)
          {
            // This is a potential "owner-neighbor", which may own one or several of
            // our local DOFs. We call a helper function, which will create a vector
            // of all *local* DOF indices, which are owned by that neighbor process.
            // Note that it is perfectly legit if the vector is empty, as this simply
            // means that all our local DOFs, which are shared with that particular
            // neighbor, are owned by some *other* neighbor process(es) and not by
            // the particular neighbor that we are currently considering.
            std::vector<IndexType> v = _owner_vidx(mirror, dof_owners, neighbor_rank);
            if(!v.empty())
            {
              // That neighbor owns at least one of our local DOFs, so create a mirror
              // from the index-vector and push it into our list of "owner-neighbors".
              // Note that the indices in the mirror are *local* DOF indices.
              _neighbors_owner.emplace_back(std::make_pair(neighbor_rank, _vidx2mirror(num_local_dofs, v)));
            }
          }
          else
          {
            // This is a potential "donee-neighbor", which shares one or several of
            // our local DOFs, which we might own. We call a helper function, which
            // will create a vector of all our *owned* DOF indices, which are shared
            // with that particular neighbor process.
            // Note that it is perfectly legit if the vector is empty, as this simply
            // means that all our local DOFs, which are shared with that particular
            // neighbor, are owned by some *other* neighbor process(es) and not by
            // this process (i.e. us).
            std::vector<IndexType> v = _donee_vidx(mirror, own_dof_idx);
            if(!v.empty())
            {
              // We own at least of one of the neighbor's local DOFs, so create a mirror
              // from the index-vector and push it into our list of "donee-neighbors".
              // Note that the indices in the mirror are *owned* DOF indices.
              _neighbors_donee.emplace_back(std::make_pair(neighbor_rank, _vidx2mirror(num_owned_dofs, v)));
            }
          }
        }
 
        // Finally, we have to determine the global DOF index for each of our local DOFs,
        // so that we can perform a "local-to-global" DOF index lookup. This information
        // is required for the assembly of matrices later on, so we have to store this in
        // a member-variable vector.
 
        // Let's create the vector of the required length, initialize all its indices to
        // ~0, which will can be used for a sanity check later on to ensure that we know
        // the global indices of all our local DOFs:
        this->_glob_dof_idx.resize(num_local_dofs, ~IndexType(0));
 
        // Next, let's loop over all DOFs that we own and assign their corresponding
        // global DOF indices, starting with our global DOF offset, which we have already
        // determined before:
        for(std::size_t i(0); i < owned_dofs.size(); ++i)
          this->_glob_dof_idx[owned_dofs[i]] = this->_glob_dof_offset + IndexType(i);
 
        // Finally, we also need to query the global DOF indices of all our local DOFs,
        // which are owned by some neighbor process (i.e. we are the donee for these DOFs).
        // Of course, we also have to send the global DOF indices of all DOFs that we own
        // to all of our "donee-neighbors".
 
        // Allocate index buffers and post receives for all of our owner-neighbors:
        const IndexType num_neigh_owner = IndexType(_neighbors_owner.size());
        Dist::RequestVector recv_reqs(num_neigh_owner);
        std::vector<std::vector<IndexType>> recv_bufs(num_neigh_owner);
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          const std::size_t num_idx = _neighbors_owner.at(i).second.num_indices();
          recv_bufs.at(i).resize(num_idx);
          recv_reqs[i] = comm.irecv(recv_bufs.at(i).data(), num_idx, _neighbors_owner.at(i).first);
        }
 
        // Allocate index buffers, fill them with the global DOF indices of our owned
        // DOFs and send them to the donee-neighbors:
        const IndexType num_neigh_donee = IndexType(_neighbors_donee.size());
        Dist::RequestVector send_reqs(num_neigh_donee);
        std::vector<std::vector<IndexType>> send_bufs(num_neigh_donee);
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          // Get the mirror for this donee-neighbor:
          const MirrorType& mirror = _neighbors_donee.at(i).second;
          const Index num_idx = mirror.num_indices();
          const IndexType* mir_idx = mirror.indices();
          // Allocate buffer of required size and translate the indices of our mirror,
          // which are "owned DOF indices", to global DOF indices, where:
          //   global_dof_index := global_dof_offset + owned_dof_index
          send_bufs.at(i).resize(num_idx);
          IndexType* jdx = send_bufs.at(i).data();
          for(IndexType j(0); j < num_idx; ++j)
            jdx[j] = this->_glob_dof_offset + mir_idx[j];
          send_reqs[i] = comm.isend(send_bufs.at(i).data(), num_idx, _neighbors_donee.at(i).first);
        }
 
        // Now process all receive requests from our owner-neighbors:
        for(std::size_t i(0u); recv_reqs.wait_any(i); )
        {
          const MirrorType& mirror = _neighbors_owner.at(i).second;
          const Index num_idx = mirror.num_indices();
          const IndexType* mir_idx = mirror.indices();
          const IndexType* jdx = recv_bufs.at(i).data();
          // Scatter the global DOF indices, which our friendly owner-neighbor has
          // provided us with, into our global DOF index vector:
          for(IndexType j(0); j < num_idx; ++j)
            this->_glob_dof_idx[mir_idx[j]] = jdx[j];
        }
 
        // wait for all sends to finish
        send_reqs.wait_all();
 
#ifdef DEBUG
        // Last but not least: debug-mode sanity check
        // We now should know the global DOF indices for all of our local DOFs:
        for(IndexType gdi : this->_glob_dof_idx)
        {
          XASSERTM(gdi != ~IndexType(0), "invalid global DOF index");
        }
#endif // DEBUG
      }
 
    private:
      // auxiliary function: create a mirror from a vector of indices
      static MirrorType _vidx2mirror(const IndexType size, const std::vector<IndexType>& vidx)
      {
        const IndexType nidx = IndexType(vidx.size());
        MirrorType mirror(size, nidx);
        IndexType* idx = mirror.indices();
        for(IndexType i(0); i < nidx; ++i)
          idx[i] = vidx[i];
        return mirror;
      }
 
      // auxiliary function: create an index vector of all *local DOFs*
      // which are owned by the process with the given neighbor rank:
      static std::vector<IndexType> _owner_vidx(const MirrorType& old_mir,
        const std::vector<int>& dof_owners, const int neighbor_rank)
      {
        std::vector<IndexType> v;
        v.reserve(old_mir.num_indices());
 
        const IndexType* old_idx = old_mir.indices();
        for(IndexType i(0); i < old_mir.num_indices(); ++i)
        {
          // is this neighbor the owner?
          if(neighbor_rank == dof_owners[old_idx[i]])
            v.push_back(old_idx[i]);
        }
        return v;
      }
 
      // auxiliary function: create an index vector all *owned DOFs*
      // which are local DOFs of the process whose mirror is given:
      static std::vector<IndexType> _donee_vidx(const MirrorType& old_mir,
        const std::vector<IndexType>& own_dof_idx)
      {
        std::vector<IndexType> v;
        v.reserve(old_mir.num_indices());
 
        const IndexType* old_idx = old_mir.indices();
        for(IndexType i(0); i < old_mir.num_indices(); ++i)
        {
          // do we own the DOF in this mirror?
          if(own_dof_idx[old_idx[i]] != ~IndexType(0))
            v.push_back(own_dof_idx[old_idx[i]]);
        }
        return v;
      }
    }; // class AlgDofParti
 
    template<typename LocalVector_, typename Mirror_>
    class AlgDofPartiVector
    {
    public:
      typedef typename LocalVector_::DataType DataType;
      typedef typename LocalVector_::IndexType IndexType;
 
      typedef LocalVector_ LocalVectorType;
      typedef Mirror_ MirrorType;
 
      typedef Global::Vector<LocalVector_, Mirror_> GlobalVectorType;
 
      typedef Global::AlgDofParti<LocalVector_, Mirror_> AlgDofPartiType;
 
      typedef LAFEM::DenseVector<DataType, IndexType> BufferVectorType;
      typedef LAFEM::DenseVector<DataType, IndexType> OwnedVectorType;
 
    protected:
      const AlgDofPartiType* _alg_dof_parti;
      OwnedVectorType _vector;
 
    public:
      AlgDofPartiVector() :
        _alg_dof_parti(nullptr),
        _vector()
      {
      }
 
      explicit AlgDofPartiVector(const AlgDofPartiType* adp_in) :
        _alg_dof_parti(adp_in),
        _vector(adp_in->get_num_owned_dofs())
      {
      }
 
      void clear()
      {
        this->_alg_dof_parti = nullptr;
        this->_vector.clear();
      }
 
      OwnedVectorType& owned()
      {
        return this->_vector;
      }
 
      const OwnedVectorType& owned() const
      {
        return this->_vector;
      }
 
      const AlgDofPartiType* get_alg_dof_parti() const
      {
        return this->_alg_dof_parti;
      }
 
      const Dist::Comm* get_comm() const
      {
        XASSERT(this->_alg_dof_parti != nullptr);
        return this->_alg_dof_parti->get_comm();
      }
 
      void allgather(BufferVectorType& vec_full) const
      {
        XASSERTM(!this->_alg_dof_parti->get_all_global_dof_counts().empty(),
          "You did not ask to assemble the required allgather data");
 
        XASSERT(vec_full.size() == this->_alg_dof_parti->get_num_global_dofs());
        this->_alg_dof_parti->get_comm()->allgatherv(this->_vector.elements(),
          this->_alg_dof_parti->get_num_owned_dofs(),
          vec_full.elements(),
          this->_alg_dof_parti->get_all_global_dof_counts().data(),
          this->_alg_dof_parti->get_all_global_dof_offsets().data());
      }
 
      void upload(const GlobalVectorType& global_vector)
      {
        this->upload(global_vector.local());
      }
 
      void upload(const LocalVectorType& local_vector)
      {
        XASSERT(this->_alg_dof_parti != nullptr);
        XASSERT(this->_vector.size() == this->_alg_dof_parti->get_num_owned_dofs());
        XASSERT(local_vector.size() == this->_alg_dof_parti->get_num_local_dofs());
 
        // Uploading is simple: we have a mirror that will extract all our owned DOFs
        // from our local vector, so we just have to call this mirror's gather function:
        this->_alg_dof_parti->get_owned_mirror().gather(this->_vector, local_vector);
      }
 
      void download(GlobalVectorType& global_vector) const
      {
        this->download(global_vector.local());
      }
 
      void download(LocalVectorType& local_vector) const
      {
        XASSERT(this->_alg_dof_parti != nullptr);
        XASSERT(this->_vector.size() == this->_alg_dof_parti->get_num_owned_dofs());
        XASSERT(local_vector.size() == this->_alg_dof_parti->get_num_local_dofs());
 
        // get our algebraic dof partitioning object
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // get the communicator
        const Dist::Comm& comm = *this->get_comm();
 
        // get the number of owner- and donee-neighbors
        const IndexType num_neigh_owner = adp.get_num_owner_neighbors();
        const IndexType num_neigh_donee = adp.get_num_donee_neighbors();
 
        // The basic idea is simple:
        // 1) Send the shared DOFs to all of our donee-neighbors
        // 2) Receive the shared DOFs from all of our owner-neighbors
        // 3) Scatter our own owned DOFs into our local vector
 
        // Create receive buffers and post receives for all owner-neighbors
        Dist::RequestVector recv_reqs(num_neigh_owner);
        std::vector<BufferVectorType> recv_bufs(num_neigh_owner);
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          // create a vector buffer
          // note: owner mirrors relate to local DOF indices
          recv_bufs.at(i) = adp.get_owner_mirror(i).create_buffer(local_vector);
          // post receive from owner neighbor
          recv_reqs[i] = comm.irecv(recv_bufs.at(i).elements(), recv_bufs.at(i).size(), adp.get_owner_rank(i));
        }
 
        // Create send buffers, fill them with our owned DOFs and send this
        // to our donee-neighbors
        Dist::RequestVector send_reqs(num_neigh_donee);
        std::vector<BufferVectorType> send_bufs(num_neigh_donee);
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          // get mirror, create buffer and gather the shared DOFs
          // note: donee mirrors relate to owned DOF indices
          const MirrorType& mir = adp.get_donee_mirror(i);
          send_bufs.at(i) = mir.create_buffer(this->_vector);
          mir.gather(send_bufs.at(i), this->_vector);
          // post send to donee neighbor
          send_reqs[i] = comm.isend(send_bufs.at(i).elements(), send_bufs.at(i).size(), adp.get_donee_rank(i));
        }
 
        // format the output vector
        local_vector.format();
 
        // scatter our own owned DOFs into our output vector
        adp.get_owned_mirror().scatter_axpy(local_vector, this->_vector);
 
        // process receives from our owner-neighbors
        for(std::size_t idx(0u); recv_reqs.wait_any(idx); )
        {
          // scatter received DOFs into our output vector
          adp.get_owner_mirror(IndexType(idx)).scatter_axpy(local_vector, recv_bufs.at(idx));
        }
 
        // at this point, all receives should have finished
        XASSERT(recv_reqs.is_null());
 
        // wait for all sends to finish
        send_reqs.wait_all();
      }
    }; // class AlgDofPartiVector
 
    template<typename LocalMatrix_, typename Mirror_>
    class AlgDofPartiMatrix
    {
    public:
      typedef typename LocalMatrix_::DataType DataType;
      typedef typename LocalMatrix_::IndexType IndexType;
 
      typedef LocalMatrix_ LocalMatrixType;
      typedef Mirror_ MirrorType;
      typedef Global::Matrix<LocalMatrixType, MirrorType, MirrorType> GlobalMatrixType;
 
      typedef typename LocalMatrixType::VectorTypeR LocalVectorType;
 
      typedef LAFEM::DenseVector<DataType, IndexType> BufferVectorType;
      typedef LAFEM::MatrixMirrorBuffer<DataType, IndexType> BufferMatrixType;
 
      typedef LAFEM::SparseMatrixCSR<DataType, IndexType> OwnedMatrixType;
 
      typedef Global::AlgDofParti<LocalVectorType, MirrorType> AlgDofPartiType;
      typedef Global::AlgDofPartiVector<LocalVectorType, MirrorType> AlgDofPartiVectorType;
 
    protected:
      const AlgDofPartiType* _alg_dof_parti;
      OwnedMatrixType _matrix;
 
      std::vector<BufferMatrixType> _donee_bufs, _owner_bufs;
      std::vector<MirrorType> _data_donee_mirs, _data_owner_mirs;
      std::vector<std::pair<IndexType,IndexType>> _my_data_mir;
 
    public:
      AlgDofPartiMatrix() :
        _alg_dof_parti(nullptr),
        _matrix()
      {
      }
 
      explicit AlgDofPartiMatrix(const AlgDofPartiType* adp_in) :
        _alg_dof_parti(adp_in),
        _matrix()
      {
      }
 
      void clear()
      {
        _alg_dof_parti = nullptr;
        _matrix.clear();
        _donee_bufs.clear();
        _owner_bufs.clear();
        _data_donee_mirs.clear();
        _data_owner_mirs.clear();
        _my_data_mir.clear();
      }
 
      OwnedMatrixType& owned()
      {
        return this->_matrix;
      }
 
      const OwnedMatrixType& owned() const
      {
        return this->_matrix;
      }
 
      const AlgDofPartiType* get_alg_dof_parti() const
      {
        return this->_alg_dof_parti;
      }
 
      const Dist::Comm* get_comm() const
      {
        XASSERT(this->_alg_dof_parti != nullptr);
        return this->_alg_dof_parti->get_comm();
      }
 
      void apply(AlgDofPartiVectorType& vec_r, const AlgDofPartiVectorType& vec_x) const
      {
        XASSERTM(!this->_alg_dof_parti->get_all_global_dof_counts().empty(),
          "You did not ask to assemble the required allgather data");
 
        // create full vector
        LocalVectorType vec_full(this->_alg_dof_parti->get_num_global_dofs());
 
        // gather full vector
        vec_x.allgather(vec_full);
 
        // apply our matrix
        this->_matrix.apply(vec_r.owned(), vec_full);
      }
 
      void upload(const GlobalMatrixType& mat_a)
      {
        upload(mat_a.local());
      }
 
      void upload(const LocalMatrixType& mat_a)
      {
        upload_symbolic(mat_a);
        upload_numeric(mat_a);
      }
 
      void upload_symbolic(const GlobalMatrixType& mat_a)
      {
        upload_symbolic(mat_a.local());
      }
 
      void upload_symbolic(const LocalMatrixType& mat_a)
      {
        XASSERT(this->_alg_dof_parti != nullptr);
 
        this->_assemble_buffers(mat_a);
        this->_assemble_structure(mat_a);
        this->_assemble_data_mirrors(mat_a);
      }
 
      void upload_numeric(const GlobalMatrixType& mat_a)
      {
        upload_numeric(mat_a.local());
      }
 
      void upload_numeric(const LocalMatrixType& mat_a)
      {
        XASSERT(this->_alg_dof_parti != nullptr);
 
        // get our partitioning
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // get our communicator
        const Dist::Comm& comm = *this->get_comm();
 
        // format our internal matrix
        this->_matrix.format();
 
        const IndexType num_neigh_owner = adp.get_num_owner_neighbors();
        const IndexType num_neigh_donee = adp.get_num_donee_neighbors();
 
        Dist::RequestVector recv_reqs(num_neigh_donee);
        Dist::RequestVector send_reqs(num_neigh_owner);
 
        std::vector<BufferVectorType> donee_vbufs(num_neigh_donee);
        std::vector<BufferVectorType> owner_vbufs(num_neigh_owner);
 
        // create receive buffers and post receives
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          // create buffer
          BufferVectorType& buf = donee_vbufs.at(i);
          buf = BufferVectorType(this->_data_donee_mirs.at(i).num_indices());
 
          // post receive
          recv_reqs[i] = comm.irecv(buf.elements(), buf.size(), adp.get_donee_rank(i));
        }
 
        // post sends
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          // create buffer
          BufferVectorType& buf = owner_vbufs.at(i);
          buf = BufferVectorType(this->_data_owner_mirs.at(i).num_indices());
 
          // gather from mirror
          this->_gather_data(buf, mat_a, this->_data_owner_mirs.at(i));
 
          // post send
          send_reqs[i] = comm.isend(buf.elements(), buf.size(), adp.get_owner_rank(i));
        }
 
        // upload our own data
        this->_upload_own(mat_a);
 
        // process all pending receives
        for(std::size_t idx(0u); recv_reqs.wait_any(idx); )
        {
          // scatter from buffer
          this->_scatter_data(this->_matrix, donee_vbufs.at(idx), this->_data_donee_mirs.at(idx));
        }
 
        // wait for all sends to finish
        send_reqs.wait_all();
      }
 
      template<typename LocalFilter_>
      void filter_matrix(const Global::Filter<LocalFilter_, MirrorType>& filter)
      {
        filter_matrix(filter.local());
      }
 
      void filter_matrix(const LAFEM::UnitFilter<DataType, IndexType>& filter)
      {
        // empty filter?
        if(filter.used_elements() <= IndexType(0))
          return;
 
        // get our partitioning
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // get global owned DOF offset and count
        const IndexType off_glob = adp.get_global_dof_offset();
        const IndexType num_glob = adp.get_num_global_dofs();
 
        // get local filter indices
        const Index num_idx = filter.used_elements();
        const IndexType* fil_idx = filter.get_indices();
 
        // get owned matrix arrays
        const IndexType* row_ptr = this->_matrix.row_ptr();
        const IndexType* col_idx = this->_matrix.col_ind();
        DataType* val = this->_matrix.val();
 
        // loop over all filter indices
        for(IndexType i(0); i < num_idx; ++i)
        {
          // translate local to global filter DOF index
          const IndexType gidx = adp.map_local_to_global_index(fil_idx[i]);
 
          // determine whether we own this DOF
          if((gidx < off_glob) || (off_glob + num_glob <= gidx))
            continue; // not my DOF
 
          // okay, that's my DOF, so filter the row
          const IndexType row = gidx - off_glob;
          for(IndexType j(row_ptr[row]); j < row_ptr[row+1]; ++j)
          {
            val[j] = DataType(col_idx[j] == gidx ? 1 : 0);
          }
        }
      }
 
    protected:
      static BufferMatrixType _asm_mat_buf(
        const LocalMatrixType& local_matrix, const MirrorType& mirror,
        const std::vector<IndexType>& gdi, const IndexType num_glob_dofs)
      {
        const Index num_idx = mirror.num_indices();
        const IndexType* mir_idx = mirror.indices();
 
        const IndexType* row_ptr = local_matrix.row_ptr();
        const IndexType* col_idx = local_matrix.col_ind();
 
        // count number of non-zeros in matrix buffer
        IndexType nnze = IndexType(0);
        for(IndexType i(0); i < num_idx; ++i)
        {
          nnze += (row_ptr[mir_idx[i]+1] - row_ptr[mir_idx[i]]);
        }
 
        // allocate matrix buffer of appropriate dimensions
        BufferMatrixType buf(num_idx, num_glob_dofs, nnze, 1);
        IndexType* buf_ptr = buf.row_ptr();
        IndexType* buf_idx = buf.col_ind();
 
        // set up buffer structure
        buf_ptr[0] = IndexType(0);
        for(IndexType i(0); i < num_idx; ++i)
        {
          const IndexType row = mir_idx[i];
          IndexType k = buf_ptr[i];
          // map local DOFs to global DOFs
          for(IndexType j(row_ptr[row]); j < row_ptr[row+1]; ++j, ++k)
            buf_idx[k] = gdi[col_idx[j]];
          buf_ptr[i+1] = k;
        }
 
        // Note: we leave the buffer indices unsorted
 
        // okay
        return buf;
      }
 
      void _assemble_buffers(const LocalMatrixType& local_matrix)
      {
        // get our partitioning
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // get out communicator
        const Dist::Comm& comm = *this->get_comm();
 
        // get number of neighbors and allocate buffer vectors
        const IndexType num_neigh_owner = adp.get_num_owner_neighbors();
        const IndexType num_neigh_donee = adp.get_num_donee_neighbors();
        _owner_bufs.resize(num_neigh_owner);
        _donee_bufs.resize(num_neigh_donee);
 
        // create the owner buffers by using the auxiliary helper function.
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          this->_owner_bufs.at(i) = _asm_mat_buf(local_matrix, adp.get_owner_mirror(i),
            adp.get_global_dof_indices(), adp.get_num_global_dofs());
        }
 
        // We have assembled our owner-neighbor matrix buffers, however, we cannot
        // assemble the donee-neighbor matrix buffers on our own. Instead, each owner
        // neighbor has to receive its donee-buffers from its donee-neighbors, because
        // only the donee knows the matrix layout of those buffers.
        // Yes, this is brain-twisting, but believe me, it cannot be realized otherwise.
 
        // Note:
        // The following code is effectively just copy-&-paste from the SynchMatrix::init()
        // function, as the same buffer problem also arises when converting global type-0
        // matrices to local type-1 matrices.
 
        Dist::RequestVector send_reqs(num_neigh_owner);
        Dist::RequestVector recv_reqs(num_neigh_donee);
 
        // receive buffer dimensions vector
        std::vector<std::array<IndexType,4>> send_dims(num_neigh_owner);
        std::vector<std::array<IndexType,4>> recv_dims(num_neigh_donee);
 
        // post send-buffer dimension receives
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          recv_reqs[i] = comm.irecv(recv_dims.at(i).data(), std::size_t(4), adp.get_donee_rank(i));
        }
 
        // send owner-buffer dimensions
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          const BufferMatrixType& sbuf = this->_owner_bufs.at(i);
          send_dims.at(i)[0] = IndexType(sbuf.rows());
          send_dims.at(i)[1] = IndexType(sbuf.columns());
          send_dims.at(i)[2] = IndexType(sbuf.entries_per_nonzero());
          send_dims.at(i)[3] = IndexType(sbuf.used_elements());
          send_reqs[i] = comm.isend(send_dims.at(i).data(), std::size_t(4), adp.get_owner_rank(i));
        }
 
        // wait for all receives to finish
        recv_reqs.wait_all();
 
        // allocate donee-buffers
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          // get the receive buffer dimensions
          IndexType nrows = recv_dims.at(i)[0];
          IndexType ncols = recv_dims.at(i)[1];
          IndexType nepnz = recv_dims.at(i)[2];
          IndexType nnze  = recv_dims.at(i)[3];
 
          // allocate receive buffer
          this->_donee_bufs.at(i) = BufferMatrixType(nrows, ncols, nnze, nepnz);
        }
 
        // post donee-buffer row-pointer array receives
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          recv_reqs[i] = comm.irecv(this->_donee_bufs.at(i).row_ptr(),
            this->_donee_bufs.at(i).rows()+std::size_t(1), adp.get_donee_rank(i));
        }
 
        // wait for all previous sends to finish
        send_reqs.wait_all();
 
        // post owner-buffer row-pointer array sends
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          send_reqs[i] = comm.isend(this->_owner_bufs.at(i).row_ptr(),
            this->_owner_bufs.at(i).rows()+std::size_t(1), adp.get_owner_rank(i));
        }
 
        // wait for all previous receives to finish
        recv_reqs.wait_all();
 
        // post donee-buffer column-index array receives
        for(IndexType i(0); i < num_neigh_donee; ++i)
        {
          recv_reqs[i] = comm.irecv(this->_donee_bufs.at(i).col_ind(),
            this->_donee_bufs.at(i).used_elements(), adp.get_donee_rank(i));
        }
 
        // wait for all previous sends to finish
        send_reqs.wait_all();
 
        // post owner-buffer column-index array sends
        for(IndexType i(0); i < num_neigh_owner; ++i)
        {
          send_reqs[i] = comm.isend(this->_owner_bufs.at(i).col_ind(),
            this->_owner_bufs.at(i).used_elements(), adp.get_owner_rank(i));
        }
 
        // wait for all receives and sends to finish
        recv_reqs.wait_all();
        send_reqs.wait_all();
      }
 
      void _assemble_structure(const LocalMatrixType& local_matrix)
      {
        // get our partitioning
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // get our matrix dimensions:
        // * each row of our matrix corresponds to one owned DOF
        // * each column corresponds to one global DOF
        const IndexType num_rows = adp.get_num_owned_dofs();
        const IndexType num_cols = adp.get_num_global_dofs();
 
        // Unfortunately, assembling the matrix structure is not that easy,
        // because it is a union of our owned matrix structure and all of
        // our donee matrix buffers. We don't have a chance to determine how
        // many duplicate (i,j) pairs we will encounter here, so we use the
        // DynamicGraph class here to keep things simple.
 
        // create a dynamic graph for our matrix
        Adjacency::DynamicGraph dynamic_graph(num_rows, num_cols);
 
        // First, let us add all local matrix rows which correspond to the
        // owned DOFs of this process, which are given by the "owned mirror"
        const IndexType* loc_row_ptr = local_matrix.row_ptr();
        const IndexType* loc_col_idx = local_matrix.col_ind();
        const IndexType* own_mir_idx = adp.get_owned_mirror().indices();
        for(IndexType i(0); i < num_rows; ++i)
        {
          // get the local DOF index of our i-th owned DOF
          const IndexType lrow = own_mir_idx[i];
          for(IndexType j(loc_row_ptr[lrow]); j < loc_row_ptr[lrow + 1]; ++j)
          {
            // translate the local column indices into global DOF indices
            dynamic_graph.insert(i, adp.map_local_to_global_index(loc_col_idx[j]));
          }
        }
 
        // process all donee-neighbor buffers
        const IndexType num_neigh_donee = adp.get_num_donee_neighbors();
        for(IndexType ineigh(0); ineigh < num_neigh_donee; ++ineigh)
        {
          // get the mirror and the matrix buffer of that donee neighbor
          // note: the donee mirror indices are owned DOF indices
          const MirrorType& mir = adp.get_donee_mirror(ineigh);
          const BufferMatrixType& buf = this->_donee_bufs.at(ineigh);
          XASSERT(mir.num_indices() == buf.rows());
          const Index num_idx = mir.num_indices();
          const IndexType* mir_idx = mir.indices();
          const IndexType* buf_ptr = buf.row_ptr();
          const IndexType* buf_idx = buf.col_ind();
 
          // loop over all matrix buffer rows
          for(IndexType i(0); i < num_idx; ++i)
          {
            // the owned DOF index of our the buffer matrix row
            const IndexType row = mir_idx[i];
 
            // loop over all buffer indices
            for(IndexType j(buf_ptr[i]); j < buf_ptr[i + 1]; ++j)
            {
              // insert entry into our local matrix
              dynamic_graph.insert(row, buf_idx[j]);
            }
          }
        }
 
        // render into a standard graph
        Adjacency::Graph graph(Adjacency::RenderType::as_is, dynamic_graph);
 
        // create our dof-partitoned matrix object
        this->_matrix = OwnedMatrixType(graph);
      }
 
      static MirrorType _asm_donee_data_mir(const LocalMatrixType& local_matrix,
        const MirrorType& mirror, const BufferMatrixType& buffer)
      {
        // allocate mirror
        const Index buf_rows = buffer.rows();
        const Index num_idx = buffer.used_elements();
        MirrorType dat_mir(local_matrix.used_elements(), num_idx);
 
        const IndexType* row_idx = mirror.indices();
        const IndexType* row_ptr = local_matrix.row_ptr();
        const IndexType* col_idx = local_matrix.col_ind();
        const IndexType* buf_ptr = buffer.row_ptr();
        const IndexType* buf_idx = buffer.col_ind();
        IndexType* dat_idx = dat_mir.indices();
 
        // loop over all buffer matrix rows
        for(IndexType i(0); i < buf_rows; ++i)
        {
          // get local matrix row index
          const IndexType row = row_idx[i];
          for(IndexType j(buf_ptr[i]); j < buf_ptr[i + 1]; ++j)
          {
            // get global column index
            const IndexType col = buf_idx[j];
 
            // loop over the corresponding row of our matrix
            for(IndexType k(row_ptr[row]); k < row_ptr[row + 1]; ++k)
            {
              // is this the column we are looking for?
              if(col_idx[k] == col)
              {
                // Yup, store this index in the mirror
                dat_idx[j] = k;
                break;
              }
            }
          }
        }
 
        return dat_mir;
      }
 
      static MirrorType _asm_owner_data_mir(const LocalMatrixType& local_matrix,
        const MirrorType& mirror, const BufferMatrixType& buffer,
        const std::vector<IndexType>& glob_dof_idx)
      {
        // allocate mirror
        const Index buf_rows = buffer.rows();
        const Index num_idx = buffer.used_elements();
        MirrorType dat_mir(local_matrix.used_elements(), num_idx);
 
        const IndexType* row_idx = mirror.indices();
        const IndexType* row_ptr = local_matrix.row_ptr();
        const IndexType* col_idx = local_matrix.col_ind();
        const IndexType* buf_ptr = buffer.row_ptr();
        const IndexType* buf_idx = buffer.col_ind();
        IndexType* dat_idx = dat_mir.indices();
 
        // loop over all buffer matrix rows
        for(IndexType i(0); i < buf_rows; ++i)
        {
          // get local matrix row index
          const IndexType row = row_idx[i];
          for(IndexType j(buf_ptr[i]); j < buf_ptr[i + 1]; ++j)
          {
            // get global column index
            const IndexType col = buf_idx[j];
 
            // loop over the corresponding row of our matrix
            for(IndexType k(row_ptr[row]); k < row_ptr[row + 1]; ++k)
            {
              // is this the column we are looking for?
              if(glob_dof_idx[col_idx[k]] == col)
              {
                // Yup, store this index in the mirror
                dat_idx[j] = k;
                break;
              }
            }
          }
        }
 
        return dat_mir;
      }
 
      void _assemble_data_mirrors(const LocalMatrixType& local_matrix)
      {
        // get our partitioning
        const AlgDofPartiType& adp = *this->_alg_dof_parti;
 
        // assemble donee data mirrors
        this->_data_donee_mirs.resize(this->_donee_bufs.size());
        for(IndexType i(0); i < this->_donee_bufs.size(); ++i)
          this->_data_donee_mirs.at(i) = _asm_donee_data_mir(this->_matrix,
            adp.get_donee_mirror(i), this->_donee_bufs.at(i));
 
        // assemble owner data mirrors
        this->_data_owner_mirs.resize(this->_owner_bufs.size());
        for(IndexType i(0); i < this->_owner_bufs.size(); ++i)
          this->_data_owner_mirs.at(i) = _asm_owner_data_mir(local_matrix,
            adp.get_owner_mirror(i), this->_owner_bufs.at(i),
            adp.get_global_dof_indices());
 
        // get matrix arrays
        const IndexType* row_ptr_x = this->_matrix.row_ptr();
        const IndexType* col_idx_x = this->_matrix.col_ind();
        const IndexType* row_ptr_a = local_matrix.row_ptr();
        const IndexType* col_idx_a = local_matrix.col_ind();
 
        // reset my data mirror
        _my_data_mir.clear();
        _my_data_mir.reserve(this->_matrix.used_elements());
 
        const IndexType num_owned_dofs = adp.get_num_owned_dofs();
        const IndexType* loc_dof_idx = adp.get_owned_mirror().indices();
        const std::vector<IndexType>& glob_dof_idx = adp.get_global_dof_indices();
 
        XASSERT(num_owned_dofs == this->_matrix.rows());
 
        // loop over all our owned DOFS
        for(IndexType own_dof(0); own_dof < num_owned_dofs; ++own_dof)
        {
          // get the local DOF index for this owned DOF
          const IndexType loc_dof = loc_dof_idx[own_dof];
 
          // loop over all columns of our owned DOF matrix
          for(IndexType j(row_ptr_x[own_dof]); j < row_ptr_x[own_dof + 1]; ++j)
          {
            // get the column index, which is a global DOF index
            const IndexType col_x = col_idx_x[j];
 
            // loop over the row of the local matrix
            for(IndexType k(row_ptr_a[loc_dof]); k < row_ptr_a[loc_dof + 1]; ++k)
            {
              // get the local column index and map it to a global one
              const IndexType col_a = glob_dof_idx[col_idx_a[k]];
 
              // match?
              if(col_a == col_x)
              {
                // Okay, add this matrix entry to our mirror
                _my_data_mir.push_back(std::make_pair(j, k));
                break;
              }
            }
          }
        }
      }
 
      void _upload_own(const LocalMatrixType& local_matrix)
      {
        // get our data arrays
        const DataType* val_a = local_matrix.val();
        DataType* val_x = this->_matrix.val();
 
        // loop over our own data mirror
        for(auto ij : this->_my_data_mir)
          val_x[ij.first] = val_a[ij.second];
      }
 
      static void _gather_data(BufferVectorType& buffer, const LocalMatrixType& local_matrix,
        const MirrorType& data_mirror)
      {
        XASSERT(buffer.size() == data_mirror.num_indices());
        XASSERT(data_mirror.size() == local_matrix.used_elements());
 
        DataType* buf_val = buffer.elements();
        const DataType* mat_val = local_matrix.val();
        const IndexType* mir_idx = data_mirror.indices();
        const Index n = buffer.size();
        for(IndexType i(0); i < n; ++i)
          buf_val[i] = mat_val[mir_idx[i]];
      }
 
      static void _scatter_data(LocalMatrixType& local_matrix, const BufferVectorType& buffer,
        const MirrorType& data_mirror)
      {
        XASSERT(buffer.size() == data_mirror.num_indices());
        XASSERT(data_mirror.size() == local_matrix.used_elements());
 
        DataType* mat_val = local_matrix.val();
        const DataType* buf_val = buffer.elements();
        const IndexType* mir_idx = data_mirror.indices();
        const Index n = buffer.size();
        for(IndexType i(0); i < n; ++i)
          mat_val[mir_idx[i]] += buf_val[i];
      }
    }; // class AlgDofPartiMatrixCSR
  } // namespace Global
} // namespace FEAT