mpi.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2019 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/numbers.h>

#include <map>
#include <numeric>
#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_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_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
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;

namespace Utilities
{
  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);

#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);

#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);

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

      ~MPI_InitFinalize();
    };

    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);

#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() == 0,
             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

      std::vector<unsigned int> send_to(objects_to_send.size());
      {
        unsigned int i = 0;
        for (const auto &m : objects_to_send)
          send_to[i++] = m.first;
      }
      AssertDimension(send_to.size(), objects_to_send.size());

      const auto receive_from =
        Utilities::MPI::compute_point_to_point_communication_pattern(comm,
                                                                     send_to);

      // 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)
          {
            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,
                                       21,
                                       comm,
                                       &buffer_send_requests[i]);
            AssertThrowMPI(ierr);
            ++i;
          }
      }

      // Receiving buffers
      std::map<unsigned int, T> received_objects;
      {
        std::vector<char> buffer;
        // We do this on a first come/first served basis
        for (unsigned int i = 0; i < receive_from.size(); ++i)
          {
            // Probe what's going on. Take data from the first available sender
            MPI_Status status;
            int        ierr = MPI_Probe(MPI_ANY_SOURCE, 21, 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, rank, 21, 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.
      MPI_Waitall(send_to.size(),
                  buffer_send_requests.data(),
                  MPI_STATUSES_IGNORE);

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

    template <typename T>
    std::vector<T>
    all_gather(const MPI_Comm &comm, const T &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);

      Assert(root_process < n_procs, ExcIndexRange(root_process, 0, 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