mpi.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2012 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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/init_finalize.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/mpi_tags.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/utilities.h>
 
#include <complex>
#include <limits>
#include <map>
#include <numeric>
#include <vector>
 
 
 
#ifdef DEAL_II_WITH_MPI
#  define DEAL_II_MPI_CONST_CAST(expr) (expr)
#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 types::global_dof_index total_size);
 
  namespace MPI
  {
    template <typename T>
    constexpr bool is_mpi_type = is_same_as_any_of<T,
                                                   char,
                                                   signed short,
                                                   signed int,
                                                   signed long,
                                                   signed long long,
                                                   signed char,
                                                   unsigned char,
                                                   unsigned short,
                                                   unsigned int,
                                                   unsigned long int,
                                                   unsigned long long,
                                                   float,
                                                   double,
                                                   long double,
                                                   bool,
                                                   std::complex<float>,
                                                   std::complex<double>,
                                                   std::complex<long double>,
                                                   wchar_t>::value;
 
    unsigned int
    n_mpi_processes(const MPI_Comm mpi_communicator);
 
    unsigned int
    this_mpi_process(const MPI_Comm mpi_communicator);
 
    std::vector<unsigned int>
    mpi_processes_within_communicator(const MPI_Comm comm_large,
                                      const MPI_Comm comm_small);
 
    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);
      }
 
      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;
      };
 
      CollectiveMutex();
 
      ~CollectiveMutex();
 
      void
      lock(const MPI_Comm comm);
 
      void
      unlock(const MPI_Comm comm);
 
    private:
      bool locked;
 
      MPI_Request request;
    };
 
 
 
    template <typename T>
    class Future
    {
    public:
      template <typename W, typename G>
      Future(W &&wait_operation, G &&get_and_cleanup_operation);
 
      Future(const Future &) = delete;
 
      Future(Future &&) noexcept = default;
 
      ~Future();
 
      Future &
      operator=(const Future &) = delete;
 
      Future &
      operator=(Future &&) noexcept = default;
 
      void
      wait();
 
      T
      get();
 
    private:
      std::function<void()> wait_function;
      std::function<T()>    get_and_cleanup_function;
 
      bool is_done;
 
      bool get_was_called;
    };
 
    std::vector<IndexSet>
    create_ascending_partitioning(
      const MPI_Comm                comm,
      const types::global_dof_index locally_owned_size);
 
    IndexSet
    create_evenly_distributed_partitioning(
      const MPI_Comm                comm,
      const types::global_dof_index 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
 
 
    std::unique_ptr<MPI_Datatype, void (*)(MPI_Datatype *)>
    create_mpi_data_type_n_bytes(const std::size_t n_bytes);
 
    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);
 
    template <typename T>
    T
    logical_or(const T &t, const MPI_Comm mpi_communicator);
 
    template <typename T, typename U>
    void
    logical_or(const T &values, const MPI_Comm mpi_communicator, U &results);
 
    template <typename T>
    void
    logical_or(const ArrayView<const T> &values,
               const MPI_Comm            mpi_communicator,
               const ArrayView<T>       &results);
 
    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 InitFinalize
    {
    public:
      MPI_InitFinalize(
        int               &argc,
        char            **&argv,
        const unsigned int max_num_threads = numbers::invalid_unsigned_int);
 
      ~MPI_InitFinalize() = default;
    };
 
    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);
 
    template <typename T>
    T
    scatter(const MPI_Comm        comm,
            const std::vector<T> &objects_to_send,
            const unsigned int    root_process = 0);
 
    template <typename T>
    T
    broadcast(const MPI_Comm     comm,
              const T           &object_to_send,
              const unsigned int root_process = 0);
 
    template <typename T>
    void
    broadcast(T                 *buffer,
              const size_t       count,
              const unsigned int root,
              const MPI_Comm     comm);
 
    template <typename T>
    T
    reduce(const T                                      &local_value,
           const MPI_Comm                                comm,
           const std::function<T(const T &, const T &)> &combiner,
           const unsigned int                            root_process = 0);
 
 
    template <typename T, typename = std::enable_if_t<is_mpi_type<T> == true>>
    std::pair<T, T>
    partial_and_total_sum(const T &value, const MPI_Comm comm);
 
 
    template <typename T>
    T
    all_reduce(const T                                      &local_value,
               const MPI_Comm                                comm,
               const std::function<T(const T &, const T &)> &combiner);
 
 
    template <typename T>
    Future<void>
    isend(const T           &object,
          MPI_Comm           communicator,
          const unsigned int target_rank,
          const unsigned int mpi_tag = 0);
 
 
    template <typename T>
    Future<T>
    irecv(MPI_Comm           communicator,
          const unsigned int source_rank,
          const unsigned int mpi_tag = 0);
 
 
    std::vector<unsigned int>
    compute_index_owner(const IndexSet &owned_indices,
                        const IndexSet &indices_to_look_up,
                        const MPI_Comm  comm);
 
    std::pair<std::vector<unsigned int>, std::map<unsigned int, IndexSet>>
    compute_index_owner_and_requesters(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);
 
 
 
    /* --------------------------- inline functions ------------------------- */
 
    namespace internal
    {
      namespace MPIDataTypes
      {
#ifdef DEAL_II_WITH_MPI
        inline MPI_Datatype
        mpi_type_id(const bool *)
        {
          return MPI_CXX_BOOL;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const char *)
        {
          return MPI_CHAR;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const signed char *)
        {
          return MPI_SIGNED_CHAR;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const wchar_t *)
        {
          return MPI_WCHAR;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const short *)
        {
          return MPI_SHORT;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const int *)
        {
          return MPI_INT;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const long int *)
        {
          return MPI_LONG;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const long long int *)
        {
          return MPI_LONG_LONG;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const unsigned char *)
        {
          return MPI_UNSIGNED_CHAR;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const unsigned short *)
        {
          return MPI_UNSIGNED_SHORT;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const unsigned int *)
        {
          return MPI_UNSIGNED;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const unsigned long int *)
        {
          return MPI_UNSIGNED_LONG;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const unsigned long long int *)
        {
          return MPI_UNSIGNED_LONG_LONG;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const float *)
        {
          return MPI_FLOAT;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const double *)
        {
          return MPI_DOUBLE;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const long double *)
        {
          return MPI_LONG_DOUBLE;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const std::complex<float> *)
        {
          return MPI_COMPLEX;
        }
 
 
 
        inline MPI_Datatype
        mpi_type_id(const std::complex<double> *)
        {
          return MPI_DOUBLE_COMPLEX;
        }
#endif
      } // namespace MPIDataTypes
    }   // namespace internal
 
 
 
#ifdef DEAL_II_WITH_MPI
    template <typename T>
    inline const MPI_Datatype mpi_type_id_for_type =
      internal::MPIDataTypes::mpi_type_id(
        static_cast<std::remove_cv_t<std::remove_reference_t<T>> *>(nullptr));
#endif
 
#ifndef DOXYGEN
    namespace internal
    {
      // declaration for an internal function that lives in mpi.templates.h
      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);
    } // namespace internal
 
 
    template <typename T>
    template <typename W, typename G>
    Future<T>::Future(W &&wait_operation, G &&get_and_cleanup_operation)
      : wait_function(wait_operation)
      , get_and_cleanup_function(get_and_cleanup_operation)
      , is_done(false)
      , get_was_called(false)
    {}
 
 
 
    template <typename T>
    Future<T>::~Future()
    {
      // If there is a clean-up function, and if it has not been
      // called yet, then do so. Note that we may not have a
      // clean-up function (not even an empty one) if the current
      // object has been moved from, into another object, and as
      // a consequence the std::function objects are now empty
      // even though they were initialized in the constructor.
      // (A std::function object whose object is a an empty lambda
      // function, [](){}, is not an empty std::function object.)
      if ((get_was_called == false) && get_and_cleanup_function)
        get();
    }
 
 
 
    template <typename T>
    void
    Future<T>::wait()
    {
      if (is_done == false)
        {
          wait_function();
 
          is_done = true;
        }
    }
 
 
    template <typename T>
    T
    Future<T>::get()
    {
      Assert(get_was_called == false,
             ExcMessage(
               "You can't call get() more than once on a Future object."));
      get_was_called = true;
 
      wait();
      return get_and_cleanup_function();
    }
 
 
 
    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, unsigned int N>
    void
    logical_or(const T (&values)[N],
               const MPI_Comm mpi_communicator,
               T (&results)[N])
    {
      static_assert(std::is_integral_v<T>,
                    "The MPI_LOR operation only allows integral data types.");
 
      internal::all_reduce(MPI_LOR,
                           ArrayView<const T>(values, N),
                           mpi_communicator,
                           ArrayView<T>(results, 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.empty(),
             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_expected_incoming_messages =
        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.empty() && n_expected_incoming_messages == 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,
                                                   /*allow_compression=*/false);
              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_expected_incoming_messages; ++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,
                                   /*allow_compression=*/false);
          }
      }
 
      // 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, typename>
    std::pair<T, T>
    partial_and_total_sum(const T &value, const MPI_Comm comm)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      return {0, value};
#  else
      if (Utilities::MPI::n_mpi_processes(comm) == 1)
        return {0, value};
      else
        {
          T prefix = {};
 
          // First obtain every process's prefix sum:
          int ierr =
            MPI_Exscan(&value,
                       &prefix,
                       1,
                       Utilities::MPI::mpi_type_id_for_type<decltype(value)>,
                       MPI_SUM,
                       comm);
          AssertThrowMPI(ierr);
 
          // Then we also need the total sum. We could obtain it by
          // calling Utilities::MPI::sum(), but it is cheaper if we
          // broadcast it from the last process, which can compute it
          // from its own prefix sum plus its own value.
          const T sum = Utilities::MPI::broadcast(
            comm, prefix + value, Utilities::MPI::n_mpi_processes(comm) - 1);
 
          return {prefix, sum};
        }
#  endif
    }
 
 
 
    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
      int ierr = MPI_Allgather(
        &n_local_data, 1, MPI_INT, size_all_data.data(), 1, MPI_INT, comm);
      AssertThrowMPI(ierr);
 
      // 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());
 
      ierr = MPI_Allgatherv(buffer.data(),
                            n_local_data,
                            MPI_CHAR,
                            received_unrolled_buffer.data(),
                            size_all_data.data(),
                            rdispls.data(),
                            MPI_CHAR,
                            comm);
      AssertThrowMPI(ierr);
 
      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
    }
 
 
 
    template <typename T>
    T
    scatter(const MPI_Comm        comm,
            const std::vector<T> &objects_to_send,
            const unsigned int    root_process)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      (void)root_process;
 
      AssertDimension(objects_to_send.size(), 1);
 
      return objects_to_send[0];
#  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);
      AssertThrow(
        (my_rank != root_process && objects_to_send.empty()) ||
          objects_to_send.size() == n_procs,
        ExcMessage(
          "The number of objects to be scattered must correspond to the number processes."));
 
      std::vector<char> send_buffer;
      std::vector<int>  send_counts;
      std::vector<int>  send_displacements;
 
      if (my_rank == root_process)
        {
          send_counts.resize(n_procs, 0);
          send_displacements.resize(n_procs + 1, 0);
 
          for (unsigned int i = 0; i < n_procs; ++i)
            {
              const auto packed_data = Utilities::pack(objects_to_send[i]);
              send_buffer.insert(send_buffer.end(),
                                 packed_data.begin(),
                                 packed_data.end());
              send_counts[i] = packed_data.size();
            }
 
          for (unsigned int i = 0; i < n_procs; ++i)
            send_displacements[i + 1] = send_displacements[i] + send_counts[i];
        }
 
      int n_local_data;
      int ierr = MPI_Scatter(send_counts.data(),
                             1,
                             MPI_INT,
                             &n_local_data,
                             1,
                             MPI_INT,
                             root_process,
                             comm);
      AssertThrowMPI(ierr);
 
      std::vector<char> recv_buffer(n_local_data);
 
      ierr = MPI_Scatterv(send_buffer.data(),
                          send_counts.data(),
                          send_displacements.data(),
                          MPI_CHAR,
                          recv_buffer.data(),
                          n_local_data,
                          MPI_CHAR,
                          root_process,
                          comm);
      AssertThrowMPI(ierr);
 
      return Utilities::unpack<T>(recv_buffer);
#  endif
    }
 
 
    template <typename T>
    void
    broadcast(T                 *buffer,
              const size_t       count,
              const unsigned int root,
              const MPI_Comm     comm)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)buffer;
      (void)count;
      (void)root;
      (void)comm;
#  else
      Assert(root < n_mpi_processes(comm),
             ExcMessage("Invalid root rank specified."));
 
      // MPI_Bcast's count is a signed int, so send at most 2^31 in each
      // iteration:
      const size_t max_send_count = std::numeric_limits<signed int>::max();
 
      size_t total_sent_count = 0;
      while (total_sent_count < count)
        {
          const size_t current_count =
            std::min(count - total_sent_count, max_send_count);
 
          const int ierr = MPI_Bcast(buffer + total_sent_count,
                                     current_count,
                                     mpi_type_id_for_type<decltype(*buffer)>,
                                     root,
                                     comm);
          AssertThrowMPI(ierr);
          total_sent_count += current_count;
        }
#  endif
    }
 
 
 
    template <typename T>
    T
    broadcast(const MPI_Comm     comm,
              const T           &object_to_send,
              const unsigned int root_process)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)comm;
      (void)root_process;
      return object_to_send;
#  else
      const auto n_procs = ::Utilities::MPI::n_mpi_processes(comm);
      AssertIndexRange(root_process, n_procs);
      (void)n_procs;
 
      if constexpr (is_mpi_type<T>)
        {
          T   object = object_to_send;
          int ierr =
            MPI_Bcast(&object, 1, mpi_type_id_for_type<T>, root_process, comm);
          AssertThrowMPI(ierr);
 
          return object;
        }
      else
        {
          std::vector<char> buffer;
          std::size_t       buffer_size = numbers::invalid_size_type;
 
          // On the root process, pack the data and determine what the
          // buffer size needs to be.
          if (this_mpi_process(comm) == root_process)
            {
              buffer      = Utilities::pack(object_to_send, false);
              buffer_size = buffer.size();
            }
 
          // Exchange the size of buffer
          int ierr = MPI_Bcast(&buffer_size,
                               1,
                               mpi_type_id_for_type<decltype(buffer_size)>,
                               root_process,
                               comm);
          AssertThrowMPI(ierr);
 
          // If not on the root process, correctly size the buffer to
          // receive the data, then do exactly that.
          if (this_mpi_process(comm) != root_process)
            buffer.resize(buffer_size);
 
          broadcast(buffer.data(), buffer_size, root_process, comm);
 
          if (Utilities::MPI::this_mpi_process(comm) == root_process)
            return object_to_send;
          else
            return Utilities::unpack<T>(buffer, false);
        }
#  endif
    }
 
 
 
    template <typename T>
    Future<void>
    isend(const T           &object,
          MPI_Comm           communicator,
          const unsigned int target_rank,
          const unsigned int mpi_tag)
    {
#  ifndef DEAL_II_WITH_MPI
      Assert(false, ExcNeedsMPI());
      (void)object;
      (void)communicator;
      (void)target_rank;
      (void)mpi_tag;
      return Future<void>([]() {}, []() {});
#  else
      // Create a pointer to a send buffer into which we pack the object
      // to be sent. The buffer will be released by the Future object once
      // the send has been verified to have succeeded.
      //
      // Conceptually, we would like this send buffer to be a
      // std::unique_ptr object whose ownership is later handed over
      // to the cleanup function. That has the disadvantage that the
      // cleanup object is a non-copyable lambda capture, leading to
      // awkward semantics. Instead, we use a std::shared_ptr; we move
      // this shared pointer into the cleanup function, which means
      // that there is exactly one shared pointer who owns the buffer
      // at any given time, though the latter is not an important
      // optimization.
      std::shared_ptr<std::vector<char>> send_buffer =
        std::make_unique<std::vector<char>>(Utilities::pack(object, false));
 
      // Now start the send, and store the result in a request object that
      // we can then wait for later:
      MPI_Request request;
      const int   ierr =
        MPI_Isend(send_buffer->data(),
                  send_buffer->size(),
                  mpi_type_id_for_type<decltype(*send_buffer->data())>,
                  target_rank,
                  mpi_tag,
                  communicator,
                  &request);
      AssertThrowMPI(ierr);
 
      // Then return a std::future-like object that has a wait()
      // function one can use to wait for the communication to finish,
      // and that has a cleanup function to be called at some point
      // after that makes sure the send buffer gets deallocated. This
      // cleanup function takes over ownership of the send buffer.
      //
      // Note that the body of the lambda function of the clean-up
      // function could be left empty. If that were so, once the
      // lambda function object goes out of scope, the 'send_buffer'
      // member of the closure object goes out of scope as well and so
      // the send_buffer is destroyed. But we may want to release the
      // buffer itself as early as possible, and so we clear the
      // buffer when the Future::get() function is called.
      auto wait = [request]() mutable {
        const int ierr = MPI_Wait(&request, MPI_STATUS_IGNORE);
        AssertThrowMPI(ierr);
      };
      auto cleanup = [send_buffer = std::move(send_buffer)]() {
        send_buffer->clear();
      };
      return Future<void>(wait, cleanup);
#  endif
    }
 
 
 
    template <typename T>
    Future<T>
    irecv(MPI_Comm           communicator,
          const unsigned int source_rank,
          const unsigned int mpi_tag)
    {
#  ifndef DEAL_II_WITH_MPI
      Assert(false, ExcNeedsMPI());
      (void)communicator;
      (void)source_rank;
      (void)mpi_tag;
      return Future<void>([]() {}, []() { return T{}; });
#  else
      // Use a 'probe' operation for the 'wait' operation of the
      // Future this function returns. It will trigger whenever we get
      // the incoming message. Later, once we have received the message, we
      // can query its size and allocate a receiver buffer.
      //
      // Since we may be waiting for multiple messages from the same
      // incoming process (with possibly the same tag -- we can't
      // know), we must make sure that the 'probe' operation we have
      // here (and which we use to determine the buffer size) matches
      // the 'recv' operation with which we actually get the data
      // later on. This is exactly what the 'MPI_Mprobe' function and
      // its 'I'mmediate variant is there for, coupled with the
      // 'MPI_Mrecv' call that would put into the clean-up function
      // below.
      std::shared_ptr<MPI_Message> message = std::make_shared<MPI_Message>();
      std::shared_ptr<MPI_Status>  status  = std::make_shared<MPI_Status>();
 
      auto wait = [source_rank, mpi_tag, communicator, message, status]() {
        const int ierr = MPI_Mprobe(
          source_rank, mpi_tag, communicator, message.get(), status.get());
        AssertThrowMPI(ierr);
      };
 
 
      // Now also define the function that actually gets the data:
      auto get = [status, message]() {
        int number_amount;
        int ierr;
        ierr = MPI_Get_count(status.get(), MPI_CHAR, &number_amount);
        AssertThrowMPI(ierr);
 
        std::vector<char> receive_buffer(number_amount);
 
        // Then actually get the data, using the matching MPI_Mrecv to the above
        // MPI_Mprobe:
        ierr = MPI_Mrecv(receive_buffer.data(),
                         number_amount,
                         mpi_type_id_for_type<decltype(*receive_buffer.data())>,
                         message.get(),
                         status.get());
        AssertThrowMPI(ierr);
 
        // Return the unpacked object:
        return Utilities::unpack<T>(receive_buffer, false);
      };
 
      return Future<T>(wait, get);
#  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