mpi.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_mpi_h
#define dealii_mpi_h

#include <deal.II/base/config.h>

#include <deal.II/base/array_view.h>
#include <deal.II/base/mpi_tags.h>
#include <deal.II/base/numbers.h>

#include <map>
#include <numeric>
#include <set>
#include <vector>

#if !defined(DEAL_II_WITH_MPI) && !defined(DEAL_II_WITH_PETSC)
// without MPI, we would still like to use
// some constructs with MPI data
// types. Therefore, create some dummies
using MPI_Comm     = int;
using MPI_Request  = int;
using MPI_Datatype = int;
using MPI_Op       = int;
#  ifndef MPI_COMM_WORLD
#    define MPI_COMM_WORLD 0
#  endif
#  ifndef MPI_COMM_SELF
#    define MPI_COMM_SELF 0
#  endif
#  ifndef MPI_REQUEST_NULL
#    define MPI_REQUEST_NULL 0
#  endif
#  ifndef MPI_MIN
#    define MPI_MIN 0
#  endif
#  ifndef MPI_MAX
#    define MPI_MAX 0
#  endif
#  ifndef MPI_SUM
#    define MPI_SUM 0
#  endif
#endif



#ifdef DEAL_II_WITH_MPI
#  if DEAL_II_MPI_VERSION_GTE(3, 0)

#    define DEAL_II_MPI_CONST_CAST(expr) (expr)

#  else

#    include <type_traits>

#    define DEAL_II_MPI_CONST_CAST(expr)     \
      const_cast<typename std::remove_const< \
        typename std::remove_pointer<decltype(expr)>::type>::type *>(expr)

#  endif
#endif



DEAL_II_NAMESPACE_OPEN


// Forward type declarations to allow MPI sums over tensorial types
#ifndef DOXYGEN
template <int rank, int dim, typename Number>
class Tensor;
template <int rank, int dim, typename Number>
class SymmetricTensor;
template <typename Number>
class SparseMatrix;
class IndexSet;
#endif

namespace Utilities
{
  IndexSet
  create_evenly_distributed_partitioning(const unsigned int my_partition_id,
                                         const unsigned int n_partitions,
                                         const IndexSet::size_type total_size);

  namespace MPI
  {
    unsigned int
    n_mpi_processes(const MPI_Comm &mpi_communicator);

    unsigned int
    this_mpi_process(const MPI_Comm &mpi_communicator);

    std::vector<unsigned int>
    compute_point_to_point_communication_pattern(
      const MPI_Comm &                 mpi_comm,
      const std::vector<unsigned int> &destinations);

    unsigned int
    compute_n_point_to_point_communications(
      const MPI_Comm &                 mpi_comm,
      const std::vector<unsigned int> &destinations);

    MPI_Comm
    duplicate_communicator(const MPI_Comm &mpi_communicator);

    void
    free_communicator(MPI_Comm &mpi_communicator);

    class DuplicatedCommunicator
    {
    public:
      explicit DuplicatedCommunicator(const MPI_Comm &communicator)
        : comm(duplicate_communicator(communicator))
      {}

      DuplicatedCommunicator(const DuplicatedCommunicator &) = delete;

      ~DuplicatedCommunicator()
      {
        free_communicator(comm);
      }

      const MPI_Comm &operator*() const
      {
        return comm;
      }


      DuplicatedCommunicator &
      operator=(const DuplicatedCommunicator &) = delete;

    private:
      MPI_Comm comm;
    };

    class CollectiveMutex
    {
    public:
      class ScopedLock
      {
      public:
        explicit ScopedLock(CollectiveMutex &mutex, const MPI_Comm &comm)
          : mutex(mutex)
          , comm(comm)
        {
          mutex.lock(comm);
        }

        ~ScopedLock()
        {
          mutex.unlock(comm);
        }

      private:
        CollectiveMutex &mutex;
        const MPI_Comm comm;
      };

      explicit CollectiveMutex();

      ~CollectiveMutex();

      void
      lock(MPI_Comm comm);

      void
      unlock(MPI_Comm comm);

    private:
      bool locked;

      MPI_Request request;
    };



#ifdef DEAL_II_WITH_MPI
    int
    create_group(const MPI_Comm & comm,
                 const MPI_Group &group,
                 const int        tag,
                 MPI_Comm *       new_comm);
#endif

    std::vector<IndexSet>
    create_ascending_partitioning(const MPI_Comm &          comm,
                                  const IndexSet::size_type local_size);

    IndexSet
    create_evenly_distributed_partitioning(
      const MPI_Comm &          comm,
      const IndexSet::size_type total_size);

#ifdef DEAL_II_WITH_MPI

    template <class Iterator, typename Number = long double>
    std::pair<Number, typename numbers::NumberTraits<Number>::real_type>
    mean_and_standard_deviation(const Iterator  begin,
                                const Iterator  end,
                                const MPI_Comm &comm);
#endif

    template <typename T>
    T
    sum(const T &t, const MPI_Comm &mpi_communicator);

    template <typename T, typename U>
    void
    sum(const T &values, const MPI_Comm &mpi_communicator, U &sums);

    template <typename T>
    void
    sum(const ArrayView<const T> &values,
        const MPI_Comm &          mpi_communicator,
        const ArrayView<T> &      sums);

    template <int rank, int dim, typename Number>
    SymmetricTensor<rank, dim, Number>
    sum(const SymmetricTensor<rank, dim, Number> &local,
        const MPI_Comm &                          mpi_communicator);

    template <int rank, int dim, typename Number>
    Tensor<rank, dim, Number>
    sum(const Tensor<rank, dim, Number> &local,
        const MPI_Comm &                 mpi_communicator);

    template <typename Number>
    void
    sum(const SparseMatrix<Number> &local,
        const MPI_Comm &            mpi_communicator,
        SparseMatrix<Number> &      global);

    template <typename T>
    T
    max(const T &t, const MPI_Comm &mpi_communicator);

    template <typename T, typename U>
    void
    max(const T &values, const MPI_Comm &mpi_communicator, U &maxima);

    template <typename T>
    void
    max(const ArrayView<const T> &values,
        const MPI_Comm &          mpi_communicator,
        const ArrayView<T> &      maxima);

    template <typename T>
    T
    min(const T &t, const MPI_Comm &mpi_communicator);

    template <typename T, typename U>
    void
    min(const T &values, const MPI_Comm &mpi_communicator, U &minima);

    template <typename T>
    void
    min(const ArrayView<const T> &values,
        const MPI_Comm &          mpi_communicator,
        const ArrayView<T> &      minima);

    struct MinMaxAvg
    {
      double sum;

      double min;

      double max;

      unsigned int min_index;

      unsigned int max_index;

      double avg;
    };

    MinMaxAvg
    min_max_avg(const double my_value, const MPI_Comm &mpi_communicator);

    std::vector<MinMaxAvg>
    min_max_avg(const std::vector<double> &my_value,
                const MPI_Comm &           mpi_communicator);


    void
    min_max_avg(const ArrayView<const double> &my_values,
                const ArrayView<MinMaxAvg> &   result,
                const MPI_Comm &               mpi_communicator);


    class MPI_InitFinalize
    {
    public:
      MPI_InitFinalize(
        int &              argc,
        char **&           argv,
        const unsigned int max_num_threads = numbers::invalid_unsigned_int);

      ~MPI_InitFinalize();

      static void
      register_request(MPI_Request &request);

      static void
      unregister_request(MPI_Request &request);

    private:
      static std::set<MPI_Request *> requests;
    };

    bool
    job_supports_mpi();

    template <typename T>
    std::map<unsigned int, T>
    some_to_some(const MPI_Comm &                 comm,
                 const std::map<unsigned int, T> &objects_to_send);

    template <typename T>
    std::vector<T>
    all_gather(const MPI_Comm &comm, const T &object_to_send);

    template <typename T>
    std::vector<T>
    gather(const MPI_Comm &   comm,
           const T &          object_to_send,
           const unsigned int root_process = 0);

    std::vector<unsigned int>
    compute_index_owner(const IndexSet &owned_indices,
                        const IndexSet &indices_to_look_up,
                        const MPI_Comm &comm);

    template <typename T>
    std::vector<T>
    compute_set_union(const std::vector<T> &vec, const MPI_Comm &comm);

    template <typename T>
    std::set<T>
    compute_set_union(const std::set<T> &set, const MPI_Comm &comm);

#ifndef DOXYGEN
    // declaration for an internal function that lives in mpi.templates.h
    namespace internal
    {
      template <typename T>
      void
      all_reduce(const MPI_Op &            mpi_op,
                 const ArrayView<const T> &values,
                 const MPI_Comm &          mpi_communicator,
                 const ArrayView<T> &      output);
    }

    // Since these depend on N they must live in the header file
    template <typename T, unsigned int N>
    void
    sum(const T (&values)[N], const MPI_Comm &mpi_communicator, T (&sums)[N])
    {
      internal::all_reduce(MPI_SUM,
                           ArrayView<const T>(values, N),
                           mpi_communicator,
                           ArrayView<T>(sums, N));
    }

    template <typename T, unsigned int N>
    void
    max(const T (&values)[N], const MPI_Comm &mpi_communicator, T (&maxima)[N])
    {
      internal::all_reduce(MPI_MAX,
                           ArrayView<const T>(values, N),
                           mpi_communicator,
                           ArrayView<T>(maxima, N));
    }

    template <typename T, unsigned int N>
    void
    min(const T (&values)[N], const MPI_Comm &mpi_communicator, T (&minima)[N])
    {
      internal::all_reduce(MPI_MIN,
                           ArrayView<const T>(values, N),
                           mpi_communicator,
                           ArrayView<T>(minima, N));
    }

    template <typename T>
    std::map<unsigned int, T>
    some_to_some(const MPI_Comm &                 comm,
                 const std::map<unsigned int, T> &objects_to_send)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      Assert(objects_to_send.size() < 2,
             ExcMessage("Cannot send to more than one processor."));
      Assert(objects_to_send.find(0) != objects_to_send.end() ||
               objects_to_send.size() == 0,
             ExcMessage("Can only send to myself or to nobody."));
      return objects_to_send;
#  else
      const auto my_proc = this_mpi_process(comm);

      std::map<unsigned int, T> received_objects;

      std::vector<unsigned int> send_to;
      send_to.reserve(objects_to_send.size());
      for (const auto &m : objects_to_send)
        if (m.first == my_proc)
          received_objects[my_proc] = m.second;
        else
          send_to.emplace_back(m.first);

      const unsigned int n_point_point_communications =
        Utilities::MPI::compute_n_point_to_point_communications(comm, send_to);

      // Protect the following communication:
      static CollectiveMutex      mutex;
      CollectiveMutex::ScopedLock lock(mutex, comm);

      // If we have something to send, or we expect something from other
      // processors, we need to visit one of the two scopes below. Otherwise,
      // no other action is required by this mpi process, and we can safely
      // return.
      if (send_to.size() == 0 && n_point_point_communications == 0)
        return received_objects;

      const int mpi_tag =
        internal::Tags::compute_point_to_point_communication_pattern;

      // Sending buffers
      std::vector<std::vector<char>> buffers_to_send(send_to.size());
      std::vector<MPI_Request>       buffer_send_requests(send_to.size());
      {
        unsigned int i = 0;
        for (const auto &rank_obj : objects_to_send)
          if (rank_obj.first != my_proc)
            {
              const auto &rank   = rank_obj.first;
              buffers_to_send[i] = Utilities::pack(rank_obj.second);
              const int ierr     = MPI_Isend(buffers_to_send[i].data(),
                                         buffers_to_send[i].size(),
                                         MPI_CHAR,
                                         rank,
                                         mpi_tag,
                                         comm,
                                         &buffer_send_requests[i]);
              AssertThrowMPI(ierr);
              ++i;
            }
      }

      // Fill the output map
      {
        std::vector<char> buffer;
        // We do this on a first come/first served basis
        for (unsigned int i = 0; i < n_point_point_communications; ++i)
          {
            // Probe what's going on. Take data from the first available sender
            MPI_Status status;
            int        ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, comm, &status);
            AssertThrowMPI(ierr);

            // Length of the message
            int len;
            ierr = MPI_Get_count(&status, MPI_CHAR, &len);
            AssertThrowMPI(ierr);
            buffer.resize(len);

            // Source rank
            const unsigned int rank = status.MPI_SOURCE;

            // Actually receive the message
            ierr = MPI_Recv(buffer.data(),
                            len,
                            MPI_CHAR,
                            status.MPI_SOURCE,
                            status.MPI_TAG,
                            comm,
                            MPI_STATUS_IGNORE);
            AssertThrowMPI(ierr);
            Assert(received_objects.find(rank) == received_objects.end(),
                   ExcInternalError(
                     "I should not receive again from this rank"));
            received_objects[rank] = Utilities::unpack<T>(buffer);
          }
      }

      // Wait to have sent all objects.
      const int ierr = MPI_Waitall(send_to.size(),
                                   buffer_send_requests.data(),
                                   MPI_STATUSES_IGNORE);
      AssertThrowMPI(ierr);

      return received_objects;
#  endif // deal.II with MPI
    }

    template <typename T>
    std::vector<T>
    all_gather(const MPI_Comm &comm, const T &object)
    {
      if (job_supports_mpi() == false)
        return {object};

#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      std::vector<T> v(1, object);
      return v;
#  else
      const auto n_procs = ::Utilities::MPI::n_mpi_processes(comm);

      std::vector<char> buffer = Utilities::pack(object);

      int n_local_data = buffer.size();

      // Vector to store the size of loc_data_array for every process
      std::vector<int> size_all_data(n_procs, 0);

      // Exchanging the size of each buffer
      MPI_Allgather(
        &n_local_data, 1, MPI_INT, size_all_data.data(), 1, MPI_INT, comm);

      // Now computing the displacement, relative to recvbuf,
      // at which to store the incoming buffer
      std::vector<int> rdispls(n_procs);
      rdispls[0] = 0;
      for (unsigned int i = 1; i < n_procs; ++i)
        rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];

      // Step 3: exchange the buffer:
      std::vector<char> received_unrolled_buffer(rdispls.back() +
                                                 size_all_data.back());

      MPI_Allgatherv(buffer.data(),
                     n_local_data,
                     MPI_CHAR,
                     received_unrolled_buffer.data(),
                     size_all_data.data(),
                     rdispls.data(),
                     MPI_CHAR,
                     comm);

      std::vector<T> received_objects(n_procs);
      for (unsigned int i = 0; i < n_procs; ++i)
        {
          std::vector<char> local_buffer(received_unrolled_buffer.begin() +
                                           rdispls[i],
                                         received_unrolled_buffer.begin() +
                                           rdispls[i] + size_all_data[i]);
          received_objects[i] = Utilities::unpack<T>(local_buffer);
        }

      return received_objects;
#  endif
    }

    template <typename T>
    std::vector<T>
    gather(const MPI_Comm &   comm,
           const T &          object_to_send,
           const unsigned int root_process)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      (void)root_process;
      std::vector<T> v(1, object_to_send);
      return v;
#  else
      const auto n_procs = ::Utilities::MPI::n_mpi_processes(comm);
      const auto my_rank = ::Utilities::MPI::this_mpi_process(comm);

      AssertIndexRange(root_process, n_procs);

      std::vector<char> buffer       = Utilities::pack(object_to_send);
      int               n_local_data = buffer.size();

      // Vector to store the size of loc_data_array for every process
      // only the root process needs to allocate memory for that purpose
      std::vector<int> size_all_data;
      if (my_rank == root_process)
        size_all_data.resize(n_procs, 0);

      // Exchanging the size of each buffer
      int ierr = MPI_Gather(&n_local_data,
                            1,
                            MPI_INT,
                            size_all_data.data(),
                            1,
                            MPI_INT,
                            root_process,
                            comm);
      AssertThrowMPI(ierr);

      // Now computing the displacement, relative to recvbuf,
      // at which to store the incoming buffer; only for root
      std::vector<int> rdispls;
      if (my_rank == root_process)
        {
          rdispls.resize(n_procs, 0);
          for (unsigned int i = 1; i < n_procs; ++i)
            rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
        }
      // exchange the buffer:
      std::vector<char> received_unrolled_buffer;
      if (my_rank == root_process)
        received_unrolled_buffer.resize(rdispls.back() + size_all_data.back());

      ierr = MPI_Gatherv(buffer.data(),
                         n_local_data,
                         MPI_CHAR,
                         received_unrolled_buffer.data(),
                         size_all_data.data(),
                         rdispls.data(),
                         MPI_CHAR,
                         root_process,
                         comm);
      AssertThrowMPI(ierr);

      std::vector<T> received_objects;

      if (my_rank == root_process)
        {
          received_objects.resize(n_procs);

          for (unsigned int i = 0; i < n_procs; ++i)
            {
              const std::vector<char> local_buffer(
                received_unrolled_buffer.begin() + rdispls[i],
                received_unrolled_buffer.begin() + rdispls[i] +
                  size_all_data[i]);
              received_objects[i] = Utilities::unpack<T>(local_buffer);
            }
        }
      return received_objects;
#  endif
    }


#  ifdef DEAL_II_WITH_MPI
    template <class Iterator, typename Number>
    std::pair<Number, typename numbers::NumberTraits<Number>::real_type>
    mean_and_standard_deviation(const Iterator  begin,
                                const Iterator  end,
                                const MPI_Comm &comm)
    {
      // below we do simple and straight-forward implementation. More elaborate
      // options are:
      // http://dx.doi.org/10.1145/2807591.2807644 section 3.1.2
      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online
      using Std        = typename numbers::NumberTraits<Number>::real_type;
      const Number sum = std::accumulate(begin, end, Number(0.));

      const auto size = Utilities::MPI::sum(std::distance(begin, end), comm);
      Assert(size > 0, ExcDivideByZero());
      const Number mean =
        Utilities::MPI::sum(sum, comm) / static_cast<Std>(size);
      Std sq_sum = 0.;
      std::for_each(begin, end, [&mean, &sq_sum](const Number &v) {
        sq_sum += numbers::NumberTraits<Number>::abs_square(v - mean);
      });
      sq_sum = Utilities::MPI::sum(sq_sum, comm);
      return std::make_pair(mean,
                            std::sqrt(sq_sum / static_cast<Std>(size - 1)));
    }
#  endif

#endif
  } // end of namespace MPI
} // end of namespace Utilities


DEAL_II_NAMESPACE_CLOSE

#endif