mpi.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_mpi_h
#define dealii_mpi_h

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

#include <vector>
#include <map>

#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
typedef int MPI_Comm;
typedef int MPI_Datatype;
typedef int MPI_Op;
#  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

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;

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

    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

    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[0]), 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
    }

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


DEAL_II_NAMESPACE_CLOSE

#endif