mpi.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2016 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 <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
typedef int MPI_Comm;
const int MPI_COMM_SELF = 0;
typedef int MPI_Datatype;
typedef int MPI_Op;
namespace MPI
{
  static const unsigned int UNSIGNED = 0;
  static const unsigned int LONG_DOUBLE = 0;
  static const unsigned int LONG_DOUBLE_COMPLEX = 0;
  static const unsigned int MAX = 0;
  static const unsigned int MIN = 0;
  static const unsigned int SUM = 0;
}
static const int MPI_MIN = 0;
static const int MPI_MAX = 0;
static const int MPI_SUM = 0;
#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;

//Forward type declaration to allow MPI sums over Vector<number> type
template <typename Number> class Vector;


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

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

    template <typename T, unsigned int N>
    inline
    void sum (const T (&values)[N],
              const MPI_Comm &mpi_communicator,
              T (&sums)[N]);

    template <typename T>
    inline
    void sum (const std::vector<T> &values,
              const MPI_Comm &mpi_communicator,
              std::vector<T> &sums);

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


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

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

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

    template <typename T, unsigned int N>
    inline
    void max (const T (&values)[N],
              const MPI_Comm &mpi_communicator,
              T (&maxima)[N]);

    template <typename T>
    inline
    void max (const std::vector<T> &values,
              const MPI_Comm &mpi_communicator,
              std::vector<T> &maxima);

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

    template <typename T, unsigned int N>
    inline
    void min (const T (&values)[N],
              const MPI_Comm &mpi_communicator,
              T (&minima)[N]);

    template <typename T>
    inline
    void min (const std::vector<T> &values,
              const MPI_Comm &mpi_communicator,
              std::vector<T> &minima);


    struct MinMaxAvg
    {
      // Note: We assume a POD property of this struct in the MPI calls in
      // min_max_avg
      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 ();

    namespace internal
    {
#ifdef DEAL_II_WITH_MPI

      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 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;
      }
#endif

      template <typename T>
      inline
      T all_reduce (const MPI_Op &mpi_op,
                    const T &t,
                    const MPI_Comm &mpi_communicator)
      {
#ifdef DEAL_II_WITH_MPI
        if (job_supports_mpi())
          {
            T output;
            MPI_Allreduce (const_cast<void *>(static_cast<const void *>(&t)),
                           &output, 1, internal::mpi_type_id(&t), mpi_op,
                           mpi_communicator);
            return output;
          }
        else
#endif
          {
            (void)mpi_op;
            (void)mpi_communicator;
            return t;
          }
      }

      template <typename T, unsigned int N>
      inline
      void all_reduce (const MPI_Op &mpi_op,
                       const T (&values)[N],
                       const MPI_Comm &mpi_communicator,
                       T (&output)[N])
      {
#ifdef DEAL_II_WITH_MPI
        if (job_supports_mpi())
          {
            MPI_Allreduce ((&values[0] != &output[0]
                            ?
                            const_cast<void *>(static_cast<const void *>(&values[0]))
                            :
                            MPI_IN_PLACE),
                           &output[0], N, internal::mpi_type_id(values), mpi_op,
                           mpi_communicator);
          }
        else
#endif
          {
            (void)mpi_op;
            (void)mpi_communicator;
            for (unsigned int i=0; i<N; ++i)
              output[i] = values[i];
          }
      }

      template <typename T>
      inline
      void all_reduce (const MPI_Op &mpi_op,
                       const std::vector<T> &values,
                       const MPI_Comm       &mpi_communicator,
                       std::vector<T>       &output)
      {
#ifdef DEAL_II_WITH_MPI
        if (job_supports_mpi())
          {
            output.resize (values.size());
            MPI_Allreduce ((&values[0] != &output[0]
                            ?
                            const_cast<void *>(static_cast<const void *>(&values[0]))
                            :
                            MPI_IN_PLACE),
                           &output[0], values.size(), internal::mpi_type_id((T *)0), mpi_op,
                           mpi_communicator);
          }
        else
#endif
          {
            (void)mpi_op;
            (void)mpi_communicator;
            output = values;
          }
      }

      template <typename T>
      inline
      void all_reduce (const MPI_Op    &mpi_op,
                       const Vector<T> &values,
                       const MPI_Comm  &mpi_communicator,
                       Vector<T>  &output)
      {
#ifdef DEAL_II_WITH_MPI
        if (job_supports_mpi())
          {
            if (values.begin() != output.begin())
              output.reinit (values.size());

            MPI_Allreduce ((values.begin() != output.begin()
                            ?
                            const_cast<void *>(static_cast<const void *>(values.begin()))
                            :
                            MPI_IN_PLACE),
                           output.begin(), values.size(), internal::mpi_type_id((T *)0), mpi_op,
                           mpi_communicator);
          }
        else
#endif
          {
            (void)mpi_op;
            (void)mpi_communicator;
            output = values;
          }
      }



    }


    template <typename T>
    inline
    T sum (const T &t,
           const MPI_Comm &mpi_communicator)
    {
      return internal::all_reduce(MPI_SUM, t, mpi_communicator);
    }


    template <typename T, unsigned int N>
    inline
    void sum (const T (&values)[N],
              const MPI_Comm &mpi_communicator,
              T (&sums)[N])
    {
      internal::all_reduce(MPI_SUM, values, mpi_communicator, sums);
    }


    template <typename T>
    inline
    void sum (const std::vector<T> &values,
              const MPI_Comm       &mpi_communicator,
              std::vector<T>       &sums)
    {
      internal::all_reduce(MPI_SUM, values, mpi_communicator, sums);
    }

    template <typename T>
    inline
    void sum (const Vector<T> &values,
              const MPI_Comm &mpi_communicator,
              Vector<T> &sums)
    {
      internal::all_reduce(MPI_SUM, values, mpi_communicator, sums);
    }


    template <int rank, int dim, typename Number>
    inline
    Tensor<rank,dim,Number>
    sum (const Tensor<rank,dim,Number> &local,
         const MPI_Comm &mpi_communicator)
    {
      const unsigned int n_entries = Tensor<rank,dim,Number>::n_independent_components;
      Number entries[ Tensor<rank,dim,Number>::n_independent_components ];

      for (unsigned int i=0; i< n_entries; ++i)
        entries[i] = local[ local.unrolled_to_component_indices(i) ];

      Number global_entries[ Tensor<rank,dim,Number>::n_independent_components ];
      ::Utilities::MPI::sum( entries, mpi_communicator, global_entries );

      Tensor<rank,dim,Number> global;
      for (unsigned int i=0; i< n_entries; ++i)
        global[ global.unrolled_to_component_indices(i) ] = global_entries[i];

      return global;
    }

    template <int rank, int dim, typename Number>
    inline
    SymmetricTensor<rank,dim,Number>
    sum (const SymmetricTensor<rank,dim,Number> &local,
         const MPI_Comm &mpi_communicator)
    {
      const unsigned int n_entries = SymmetricTensor<rank,dim,Number>::n_independent_components;
      Number entries[ SymmetricTensor<rank,dim,Number>::n_independent_components ];

      for (unsigned int i=0; i< n_entries; ++i)
        entries[i] = local[ local.unrolled_to_component_indices(i) ];

      Number global_entries[ SymmetricTensor<rank,dim,Number>::n_independent_components ];
      ::Utilities::MPI::sum( entries, mpi_communicator, global_entries );

      SymmetricTensor<rank,dim,Number> global;
      for (unsigned int i=0; i< n_entries; ++i)
        global[ global.unrolled_to_component_indices(i) ] = global_entries[i];

      return global;
    }

    template <typename T>
    inline
    T max (const T &t,
           const MPI_Comm &mpi_communicator)
    {
      return internal::all_reduce(MPI_MAX, t, mpi_communicator);
    }


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


    template <typename T>
    inline
    void max (const std::vector<T> &values,
              const MPI_Comm       &mpi_communicator,
              std::vector<T>       &maxima)
    {
      internal::all_reduce(MPI_MAX, values, mpi_communicator, maxima);
    }


    template <typename T>
    inline
    T min (const T &t,
           const MPI_Comm &mpi_communicator)
    {
      return internal::all_reduce(MPI_MIN, t, mpi_communicator);
    }


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


    template <typename T>
    inline
    void min (const std::vector<T> &values,
              const MPI_Comm       &mpi_communicator,
              std::vector<T>       &minima)
    {
      internal::all_reduce(MPI_MIN, values, mpi_communicator, minima);
    }


    inline
    bool job_supports_mpi ()
    {
#ifdef DEAL_II_WITH_MPI
      int MPI_has_been_started = 0;
      MPI_Initialized(&MPI_has_been_started);

      return (MPI_has_been_started > 0);
#else
      return false;
#endif
    }
  } // end of namespace MPI
} // end of namespace Utilities


DEAL_II_NAMESPACE_CLOSE

#endif