MpiWrapper.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_COMMON_MPIWRAPPER_HPP_
#define GEOS_COMMON_MPIWRAPPER_HPP_
 
#include "common/DataTypes.hpp"
#include "common/Span.hpp"
#include "common/TypesHelpers.hpp"
 
#include <numeric>
 
#if defined(GEOS_USE_MPI)
  #include <mpi.h>
#define MPI_PARAM( x ) x
#else
#define MPI_PARAM( x )
typedef int MPI_Comm;
 
#define MPI_COMM_NULL      ((MPI_Comm)0x04000000)
#define MPI_COMM_WORLD     ((MPI_Comm)0x44000000)
#define MPI_COMM_SELF      ((MPI_Comm)0x40000000)
 
 
typedef int MPI_Datatype;
#define MPI_CHAR           ((MPI_Datatype)0x4c000101)
#define MPI_SIGNED_CHAR    ((MPI_Datatype)0x4c000118)
#define MPI_UNSIGNED_CHAR  ((MPI_Datatype)0x4c000102)
#define MPI_BYTE           ((MPI_Datatype)0x4c00010d)
#define MPI_WCHAR          ((MPI_Datatype)0x4c00040e)
#define MPI_SHORT          ((MPI_Datatype)0x4c000203)
#define MPI_UNSIGNED_SHORT ((MPI_Datatype)0x4c000204)
#define MPI_INT            ((MPI_Datatype)0x4c000405)
#define MPI_UNSIGNED       ((MPI_Datatype)0x4c000406)
#define MPI_LONG           ((MPI_Datatype)0x4c000807)
#define MPI_UNSIGNED_LONG  ((MPI_Datatype)0x4c000808)
#define MPI_FLOAT          ((MPI_Datatype)0x4c00040a)
#define MPI_DOUBLE         ((MPI_Datatype)0x4c00080b)
#define MPI_LONG_DOUBLE    ((MPI_Datatype)0x4c00100c)
#define MPI_LONG_LONG_INT  ((MPI_Datatype)0x4c000809)
#define MPI_UNSIGNED_LONG_LONG ((MPI_Datatype)0x4c000819)
#define MPI_LONG_LONG      MPI_LONG_LONG_INT
 
typedef int MPI_Op;
 
#define MPI_MAX     (MPI_Op)(0x58000001)
#define MPI_MIN     (MPI_Op)(0x58000002)
#define MPI_SUM     (MPI_Op)(0x58000003)
#define MPI_PROD    (MPI_Op)(0x58000004)
#define MPI_LAND    (MPI_Op)(0x58000005)
#define MPI_BAND    (MPI_Op)(0x58000006)
#define MPI_LOR     (MPI_Op)(0x58000007)
#define MPI_BOR     (MPI_Op)(0x58000008)
#define MPI_LXOR    (MPI_Op)(0x58000009)
#define MPI_BXOR    (MPI_Op)(0x5800000a)
#define MPI_MINLOC  (MPI_Op)(0x5800000b)
#define MPI_MAXLOC  (MPI_Op)(0x5800000c)
#define MPI_REPLACE (MPI_Op)(0x5800000d)
#define MPI_NO_OP   (MPI_Op)(0x5800000e)
 
#define MPI_SUCCESS          0        /* Successful return code */
#define MPI_UNDEFINED (-32766)
#define MPI_STATUS_IGNORE (MPI_Status *)1
#define MPI_STATUSES_IGNORE (MPI_Status *)1
#define MPI_REQUEST_NULL  ((MPI_Request)0x2c000000)
typedef int MPI_Request;
 
typedef int MPI_Info;
#define MPI_INFO_NULL (MPI_Info)(0x60000000)
 
struct MPI_Status
{
  int junk;
};
 
#endif
 
#if defined(NDEBUG)
#define MPI_CHECK_ERROR( error ) ((void) error)
#else
#define MPI_CHECK_ERROR( error ) GEOS_ERROR_IF_NE( error, MPI_SUCCESS );
#endif
 
 
namespace geos
{
 
#ifdef GEOS_USE_MPI
extern MPI_Comm MPI_COMM_GEOS;
#else
extern int MPI_COMM_GEOS;
#endif
 
struct MpiWrapper
{
public:
 
  enum class Reduction
  {
    Max,  
    Min,  
    Sum,  
    Prod, 
    LogicalAnd, 
    LogicalOr, 
  };
 
  enum class PairReduction
  {
    Max, 
    Min, 
  };
 
  template< typename FIRST, typename SECOND >
  struct PairType
  {
    FIRST first;
    SECOND second;
  };
 
  MpiWrapper() = delete;
 
 
  static void barrier( MPI_Comm const & MPI_PARAM( comm )=MPI_COMM_GEOS );
 
  static int cartCoords( MPI_Comm comm, int rank, int maxdims, int coords[] );
 
  static int cartCreate( MPI_Comm comm_old, int ndims, const int dims[], const int periods[],
                         int reorder, MPI_Comm * comm_cart );
 
  static int cartRank( MPI_Comm comm, const int coords[] );
 
  static void commFree( MPI_Comm & comm );
 
  static int commRank( MPI_Comm const & MPI_PARAM( comm )=MPI_COMM_GEOS );
 
  static int commSize( MPI_Comm const & MPI_PARAM( comm )=MPI_COMM_GEOS );
 
  static bool commCompare( MPI_Comm const & comm1, MPI_Comm const & comm2 );
 
  static bool initialized();
 
  static int init( int * argc, char * * * argv );
 
  static void finalize();
 
  static MPI_Comm commDup( MPI_Comm const comm );
 
  static MPI_Comm commSplit( MPI_Comm const comm, int color, int key );
 
  static int test( MPI_Request * request, int * flag, MPI_Status * status );
 
  static int testAny( int count, MPI_Request array_of_requests[], int * idx, int * flags, MPI_Status array_of_statuses[] );
 
  static int testSome( int count, MPI_Request array_of_requests[], int * outcount, int array_of_indices[], MPI_Status array_of_statuses[] );
 
  static int testAll( int count, MPI_Request array_of_requests[], int * flags, MPI_Status array_of_statuses[] );
 
  static int check( MPI_Request * request, int * flag, MPI_Status * status );
 
  static int checkAny( int count, MPI_Request array_of_requests[], int * idx, int * flag, MPI_Status array_of_statuses[] );
 
  static int checkAll( int count, MPI_Request array_of_requests[], int * flag, MPI_Status array_of_statuses[] );
 
  static int wait( MPI_Request * request, MPI_Status * status );
 
  static int waitAny( int count, MPI_Request array_of_requests[], int * indx, MPI_Status array_of_statuses[] );
 
  static int waitSome( int count, MPI_Request array_of_requests[], int * outcount, int array_of_indices[], MPI_Status array_of_statuses[] );
 
  static int waitAll( int count, MPI_Request array_of_requests[], MPI_Status array_of_statuses[] );
 
  static double wtime( void );
 
 
  static int activeWaitAny( const int count,
                            MPI_Request array_of_requests[],
                            MPI_Status array_of_statuses[],
                            std::function< MPI_Request ( int ) > func );
 
  static int activeWaitSome( const int count,
                             MPI_Request array_of_requests[],
                             MPI_Status array_of_statuses[],
                             std::function< MPI_Request ( int ) > func );
 
  static int activeWaitSomeCompletePhase( const int participants,
                                          stdVector< std::tuple< MPI_Request *, MPI_Status *, std::function< MPI_Request ( int ) > > > const & phases );
 
  static int activeWaitOrderedCompletePhase( const int participants,
                                             stdVector< std::tuple< MPI_Request *, MPI_Status *, std::function< MPI_Request ( int ) > > > const & phases );
 
#if !defined(GEOS_USE_MPI)
  static std::map< int, std::pair< int, void * > > & getTagToPointersMap()
  {
    static std::map< int, std::pair< int, void * > > tagToPointers;
    return tagToPointers;
  }
#endif
 
  static int nodeCommSize();
 
  template< typename T_SEND, typename T_RECV >
  static int allgather( T_SEND const * sendbuf,
                        int sendcount,
                        T_RECV * recvbuf,
                        int recvcount,
                        MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T_SEND, typename T_RECV >
  static int allgatherv( T_SEND const * sendbuf,
                         int sendcount,
                         T_RECV * recvbuf,
                         int * recvcounts,
                         int * displacements,
                         MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void allGather( T const myValue, array1d< T > & allValues, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static int allGather( arrayView1d< T const > const & sendbuf,
                        array1d< T > & recvbuf,
                        MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static int allGatherv( arrayView1d< T const > const & sendbuf,
                         array1d< T > & recvbuf,
                         MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static T allReduce( T const & value, Reduction const op, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void allReduce( Span< T const > src, Span< T > dst, Reduction const op, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
  static void allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, Reduction const op, MPI_Comm const comm = MPI_COMM_GEOS );
 
  template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
  static void allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, int const count, Reduction const op, MPI_Comm const comm );
 
  template< typename FIRST, typename SECOND, PairReduction OP >
  static PairType< FIRST, SECOND > allReduce( PairType< FIRST, SECOND > const & pair,
                                              MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename FIRST, typename SECOND, typename CONTAINER, PairReduction OP >
  static PairType< FIRST, SECOND > allReduce( CONTAINER const & pairs,
                                              MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static int reduce( T const * sendbuf, T * recvbuf, int count, MPI_Op op, int root, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static T reduce( T const & value, Reduction const op, int root, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void reduce( Span< T const > src, Span< T > dst, Reduction const op, int root, MPI_Comm comm = MPI_COMM_GEOS );
 
 
  template< typename T >
  static int scan( T const * sendbuf, T * recvbuf, int count, MPI_Op op, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static int exscan( T const * sendbuf, T * recvbuf, int count, MPI_Op op, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static int bcast( T * buffer, int count, int root, MPI_Comm comm = MPI_COMM_GEOS );
 
 
  template< typename T >
  static void broadcast( T & value, int srcRank = 0, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename TS, typename TR >
  static int gather( TS const * const sendbuf,
                     int sendcount,
                     TR * const recvbuf,
                     int recvcount,
                     int root,
                     MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T, typename DST_CONTAINER,
            typename = std::enable_if_t<
              std::is_trivially_copyable_v< T > &&
              std::is_same_v< decltype(std::declval< DST_CONTAINER >().size()), std::size_t > &&
              std::is_same_v< decltype(std::declval< DST_CONTAINER >().data()), T * > > >
  static int gather( T const & value,
                     DST_CONTAINER & destValuesBuffer,
                     int root,
                     MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename TS, typename TR >
  static int gatherv( TS const * const sendbuf,
                      int sendcount,
                      TR * const recvbuf,
                      const int * recvcounts,
                      const int * displs,
                      int root,
                      MPI_Comm comm = MPI_COMM_GEOS );
 
 
  template< typename TS, typename TR >
  static int scatter( TS const * const sendbuf,
                      int sendcount,
                      TR * const recvbuf,
                      int recvcount,
                      int root,
                      MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename TS, typename TR >
  static int scatterv( TS const * const sendbuf,
                       const int * sendcounts,
                       const int * displs,
                       TR * const recvbuf,
                       int recvcount,
                       int root,
                       MPI_Comm comm = MPI_COMM_GEOS );
 
 
 
  static MPI_Op getMpiOp( Reduction const op );
 
  template< typename T >
  static int recv( array1d< T > & buf,
                   int MPI_PARAM( source ),
                   int tag,
                   MPI_Comm MPI_PARAM( comm ),
                   MPI_Status * MPI_PARAM( request ) );
 
  template< typename T >
  static int iSend( arrayView1d< T > const & buf,
                    int MPI_PARAM( dest ),
                    int tag,
                    MPI_Comm MPI_PARAM( comm ),
                    MPI_Request * MPI_PARAM( request ) );
 
  template< typename T >
  static int iRecv( T * const buf,
                    int count,
                    int source,
                    int tag,
                    MPI_Comm comm,
                    MPI_Request * request );
 
  template< typename T >
  static int send( T const * const buf,
                   int count,
                   int dest,
                   int tag,
                   MPI_Comm comm );
 
  template< typename T >
  static int iSend( T const * const buf,
                    int count,
                    int dest,
                    int tag,
                    MPI_Comm comm,
                    MPI_Request * request );
 
  template< typename U, typename T >
  static U prefixSum( T const value, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static T sum( T const & value, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void sum( Span< T const > src, Span< T > dst, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static T min( T const & value, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void min( Span< T const > src, Span< T > dst, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename FIRST, typename SECOND >
  static PairType< FIRST, SECOND > min( PairType< FIRST, SECOND > const & pair, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename FIRST, typename SECOND, typename CONTAINER >
  static PairType< FIRST, SECOND > min( CONTAINER const & pairs, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static T max( T const & value, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename T >
  static void max( Span< T const > src, Span< T > dst, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename FIRST, typename SECOND >
  static PairType< FIRST, SECOND > max( PairType< FIRST, SECOND > const & pair, MPI_Comm comm = MPI_COMM_GEOS );
 
  template< typename FIRST, typename SECOND, typename CONTAINER >
  static PairType< FIRST, SECOND > max( CONTAINER const & pairs, MPI_Comm comm = MPI_COMM_GEOS );
 
private:
 
  template< typename T >
  static int allReduce( T const * sendbuf, T * recvbuf, int count, MPI_Op op, MPI_Comm comm = MPI_COMM_GEOS );
};
 
namespace internal
{
 
struct ManagedResources
{
  // The list of managed MPI_Op instances
  std::set< MPI_Op > m_mpiOps;
 
  // The list of managed MPI_Type instances
  std::set< MPI_Datatype > m_mpiTypes;
 
  void finalize();
};
 
ManagedResources & getManagedResources();
 
template< typename T, typename ENABLE = void >
struct MpiTypeImpl {};
 
#define ADD_MPI_TYPE_MAP( T, MPI_T ) \
  template<> struct MpiTypeImpl< T > { static MPI_Datatype get() { return MPI_T; } }
 
ADD_MPI_TYPE_MAP( float, MPI_FLOAT );
ADD_MPI_TYPE_MAP( double, MPI_DOUBLE );
 
ADD_MPI_TYPE_MAP( char, MPI_CHAR );
ADD_MPI_TYPE_MAP( signed char, MPI_SIGNED_CHAR );
ADD_MPI_TYPE_MAP( unsigned char, MPI_UNSIGNED_CHAR );
 
ADD_MPI_TYPE_MAP( int, MPI_INT );
ADD_MPI_TYPE_MAP( long int, MPI_LONG );
ADD_MPI_TYPE_MAP( long long int, MPI_LONG_LONG );
 
ADD_MPI_TYPE_MAP( unsigned int, MPI_UNSIGNED );
ADD_MPI_TYPE_MAP( unsigned long int, MPI_UNSIGNED_LONG );
ADD_MPI_TYPE_MAP( unsigned long long int, MPI_UNSIGNED_LONG_LONG );
 
#undef ADD_MPI_TYPE_MAP
 
template< typename T >
struct MpiTypeImpl< T, std::enable_if_t< std::is_enum< T >::value > >
{
  static MPI_Datatype get() { return MpiTypeImpl< std::underlying_type_t< T > >::get(); }
};
 
template<>
struct MpiTypeImpl< bool * >
{
  static MPI_Datatype get()
  {
    // Return the appropriate MPI_Datatype for bool*
    return MPI_BYTE;
  }
};
 
 
template< typename T >
MPI_Datatype getMpiType()
{
  return MpiTypeImpl< T >::get();
}
 
template< typename FIRST, typename SECOND >
MPI_Datatype getMpiPairType()
{
  static_assert( "no default implementation, please add a template specialization and add it in the \"testMpiWrapper\" unit test." );
  return {};
}
template<> MPI_Datatype getMpiPairType< int, int >();
template<> MPI_Datatype getMpiPairType< long int, int >();
template<> MPI_Datatype getMpiPairType< long int, long int >();
template<> MPI_Datatype getMpiPairType< long long int, long long int >();
template<> MPI_Datatype getMpiPairType< float, int >();
template<> MPI_Datatype getMpiPairType< double, int >();
template<> MPI_Datatype getMpiPairType< double, long int >();
template<> MPI_Datatype getMpiPairType< double, long long int >();
template<> MPI_Datatype getMpiPairType< double, double >();
 
// It is advised to always use this custom operator for pairs as MPI_MAXLOC is not a true lexicographical comparator.
template< typename FIRST, typename SECOND, MpiWrapper::PairReduction OP >
MPI_Op getMpiPairReductionOp()
{
  static auto const createOpHolder = [] () {
    using PAIR_T = MpiWrapper::PairType< FIRST, SECOND >;
 
    auto const customOpFunc =
      []( void * invec, void * inoutvec, int * len, MPI_Datatype * )
    {
      for( int i = 0; i < *len; ++i )
      {
        PAIR_T & in = static_cast< PAIR_T * >(invec)[i];
        PAIR_T & inout = static_cast< PAIR_T * >(inoutvec)[i];
        if constexpr ( OP == MpiWrapper::PairReduction::Min )
        {
          if( std::tie( in.first, in.second ) < std::tie( inout.first, inout.second ) )
            inout = in;
        }
        else
        {
          if( std::tie( in.first, in.second ) > std::tie( inout.first, inout.second ) )
            inout = in;
        }
      }
    };
 
    MPI_Op mpiOp;
    GEOS_ERROR_IF_NE( MPI_Op_create( customOpFunc, 1, &mpiOp ), MPI_SUCCESS );
    // Resource registered to be destroyed at MpiWrapper::finalize().
    internal::getManagedResources().m_mpiOps.emplace( mpiOp );
    return mpiOp;
  };
  // Static storage to ensure the MPI operation is created only once and reused for all calls to this function.
  static MPI_Op mpiOp{ createOpHolder() };
  return mpiOp;
}
 
}
 
inline MPI_Op MpiWrapper::getMpiOp( Reduction const op )
{
  switch( op )
  {
    case Reduction::Sum:
    {
      return MPI_SUM;
    }
    case Reduction::Min:
    {
      return MPI_MIN;
    }
    case Reduction::Max:
    {
      return MPI_MAX;
    }
    case Reduction::Prod:
    {
      return MPI_PROD;
    }
    case Reduction::LogicalAnd:
    {
      return MPI_LAND;
    }
    case Reduction::LogicalOr:
    {
      return MPI_LOR;
    }
    default:
      GEOS_ERROR( "Unsupported reduction operation" );
      return MPI_NO_OP;
  }
}
 
template< typename T_SEND, typename T_RECV >
int MpiWrapper::allgather( T_SEND const * const sendbuf,
                           int sendcount,
                           T_RECV * const recvbuf,
                           int recvcount,
                           MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Allgather( sendbuf, sendcount, internal::getMpiType< T_SEND >(),
                        recvbuf, recvcount, internal::getMpiType< T_RECV >(),
                        comm );
#else
  static_assert( std::is_same< T_SEND, T_RECV >::value,
                 "MpiWrapper::allgather() for serial run requires send and receive buffers are of the same type" );
  GEOS_ERROR_IF_NE_MSG( sendcount, recvcount, "sendcount is not equal to recvcount." );
  std::copy( sendbuf, sendbuf + sendcount, recvbuf )
  return 0;
#endif
}
 
template< typename T_SEND, typename T_RECV >
int MpiWrapper::allgatherv( T_SEND const * const sendbuf,
                            int sendcount,
                            T_RECV * const recvbuf,
                            int * recvcounts,
                            int * displacements,
                            MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Allgatherv( sendbuf, sendcount, internal::getMpiType< T_SEND >(),
                         recvbuf, recvcounts, displacements, internal::getMpiType< T_RECV >(),
                         comm );
#else
  static_assert( std::is_same< T_SEND, T_RECV >::value,
                 "MpiWrapper::allgatherv() for serial run requires send and receive buffers are of the same type" );
  GEOS_ERROR_IF_NE_MSG( sendcount, recvcount, "sendcount is not equal to recvcount." );
  std::copy( sendbuf, sendbuf + sendcount, recvbuf )
  return 0;
#endif
}
 
 
template< typename T >
void MpiWrapper::allGather( T const myValue, array1d< T > & allValues, MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  int const mpiSize = commSize( comm );
  allValues.resize( mpiSize );
 
  MPI_Datatype const MPI_TYPE = internal::getMpiType< T >();
 
  MPI_Allgather( &myValue, 1, MPI_TYPE, allValues.data(), 1, MPI_TYPE, comm );
 
#else
  allValues.resize( 1 );
  allValues[0] = myValue;
#endif
}
 
template< typename T >
int MpiWrapper::allGather( arrayView1d< T const > const & sendValues,
                           array1d< T > & allValues,
                           MPI_Comm MPI_PARAM( comm ) )
{
  int const sendSize = LvArray::integerConversion< int >( sendValues.size() );
#ifdef GEOS_USE_MPI
  int const mpiSize = commSize( comm );
  allValues.resize( mpiSize * sendSize );
  return MPI_Allgather( sendValues.data(),
                        sendSize,
                        internal::getMpiType< T >(),
                        allValues.data(),
                        sendSize,
                        internal::getMpiType< T >(),
                        comm );
 
#else
  allValues.resize( sendSize );
  for( localIndex a=0; a<sendSize; ++a )
  {
    allValues[a] = sendValues[a];
  }
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::allGatherv( arrayView1d< T const > const & sendValues,
                            array1d< T > & allValues,
                            MPI_Comm MPI_PARAM( comm ) )
{
  int const sendSize = LvArray::integerConversion< int >( sendValues.size() );
#ifdef GEOS_USE_MPI
  int const mpiSize = commSize( comm );
  array1d< int > counts;
  allGather( sendSize, counts, comm );
  array1d< int > displs( mpiSize + 1 );
  std::partial_sum( counts.begin(), counts.end(), displs.begin() + 1 );
  allValues.resize( displs.back() );
  return MPI_Allgatherv( sendValues.data(),
                         sendSize,
                         internal::getMpiType< T >(),
                         allValues.data(),
                         counts.data(),
                         displs.data(),
                         internal::getMpiType< T >(),
                         comm );
 
#else
  allValues.resize( sendSize );
  for( localIndex a=0; a<sendSize; ++a )
  {
    allValues[a] = sendValues[a];
  }
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::allReduce( T const * const sendbuf,
                           T * const recvbuf,
                           int const count,
                           MPI_Op const MPI_PARAM( op ),
                           MPI_Comm const MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  MPI_Datatype const mpiType = internal::getMpiType< T >();
  return MPI_Allreduce( sendbuf == recvbuf ? MPI_IN_PLACE : sendbuf, recvbuf, count, mpiType, op, comm );
#else
  if( sendbuf != recvbuf )
  {
    memcpy( recvbuf, sendbuf, count * sizeof( T ) );
  }
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::reduce( T const * const sendbuf,
                        T * const recvbuf,
                        int const count,
                        MPI_Op const MPI_PARAM( op ),
                        int root,
                        MPI_Comm const MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  MPI_Datatype const mpiType = internal::getMpiType< T >();
  return MPI_Reduce( sendbuf == recvbuf ? MPI_IN_PLACE : sendbuf, recvbuf, count, mpiType, op, root, comm );
#else
  if( sendbuf != recvbuf )
  {
    memcpy( recvbuf, sendbuf, count * sizeof( T ) );
  }
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::scan( T const * const sendbuf,
                      T * const recvbuf,
                      int count,
                      MPI_Op MPI_PARAM( op ),
                      MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Scan( sendbuf, recvbuf, count, internal::getMpiType< T >(), op, comm );
#else
  memcpy( recvbuf, sendbuf, count*sizeof(T) );
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::exscan( T const * const MPI_PARAM( sendbuf ),
                        T * const recvbuf,
                        int count,
                        MPI_Op MPI_PARAM( op ),
                        MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Exscan( sendbuf, recvbuf, count, internal::getMpiType< T >(), op, comm );
#else
  memset( recvbuf, 0, count*sizeof(T) );
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::bcast( T * const MPI_PARAM( buffer ),
                       int MPI_PARAM( count ),
                       int MPI_PARAM( root ),
                       MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Bcast( buffer, count, internal::getMpiType< T >(), root, comm );
#else
  return 0;
#endif
 
}
 
template< typename T >
void MpiWrapper::broadcast( T & MPI_PARAM( value ), int MPI_PARAM( srcRank ), MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  MPI_Bcast( &value, 1, internal::getMpiType< T >(), srcRank, comm );
#endif
}
 
template<>
inline
void MpiWrapper::broadcast< string >( string & MPI_PARAM( value ),
                                      int MPI_PARAM( srcRank ),
                                      MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  int size = LvArray::integerConversion< int >( value.size() );
  broadcast( size, srcRank, comm );
  value.resize( size );
  MPI_Bcast( const_cast< char * >( value.data() ), size, internal::getMpiType< char >(), srcRank, comm );
#endif
}
 
template< typename TS, typename TR >
int MpiWrapper::gather( TS const * const sendbuf,
                        int sendcount,
                        TR * const recvbuf,
                        int recvcount,
                        int MPI_PARAM( root ),
                        MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Gather( sendbuf, sendcount, internal::getMpiType< TS >(),
                     recvbuf, recvcount, internal::getMpiType< TR >(),
                     root, comm );
#else
  static_assert( std::is_same< TS, TR >::value,
                 "MpiWrapper::gather() for serial run requires send and receive buffers are of the same type" );
  std::size_t const sendBufferSize = sendcount * sizeof(TS);
  std::size_t const recvBufferSize = recvcount * sizeof(TR);
  GEOS_ERROR_IF_NE_MSG( sendBufferSize, recvBufferSize, "size of send buffer and receive buffer are not equal" );
  memcpy( recvbuf, sendbuf, sendBufferSize );
  return 0;
#endif
}
 
template< typename T, typename DST_CONTAINER, typename >
int MpiWrapper::gather( T const & value,
                        DST_CONTAINER & destValuesBuffer,
                        int root,
                        MPI_Comm MPI_PARAM( comm ) )
{
  if( commRank() == 0 )
    GEOS_ERROR_IF_LT_MSG( destValuesBuffer.size(), size_t( commSize() ),
                          "Receive buffer is not large enough to contain the values to receive." );
#ifdef GEOS_USE_MPI
  return MPI_Gather( &value, sizeof( T ), internal::getMpiType< uint8_t >(),
                     destValuesBuffer.data(), sizeof( T ), internal::getMpiType< uint8_t >(),
                     root, comm );
#else
  memcpy( destValuesBuffer.data(), &value, sendBufferSize );
  return 0;
#endif
}
 
template< typename TS, typename TR >
int MpiWrapper::gatherv( TS const * const sendbuf,
                         int sendcount,
                         TR * const recvbuf,
                         const int * recvcounts,
                         const int * MPI_PARAM( displs ),
                         int MPI_PARAM( root ),
                         MPI_Comm MPI_PARAM( comm ) )
{
#ifdef GEOS_USE_MPI
  return MPI_Gatherv( sendbuf, sendcount, internal::getMpiType< TS >(),
                      recvbuf, recvcounts, displs, internal::getMpiType< TR >(),
                      root, comm );
#else
  static_assert( std::is_same< TS, TR >::value,
                 "MpiWrapper::gather() for serial run requires send and receive buffers are of the same type" );
  std::size_t const sendBufferSize = sendcount * sizeof(TS);
  std::size_t const recvBufferSize = recvcounts[0] * sizeof(TR);
  GEOS_ERROR_IF_NE_MSG( sendBufferSize, recvBufferSize, "size of send buffer and receive buffer are not equal" );
  memcpy( recvbuf, sendbuf, sendBufferSize );
  return 0;
#endif
}
 
 
template< typename TS, typename TR >
int MpiWrapper::scatter( TS const * const sendbuf,
                         int sendcount,
                         TR * const recvbuf,
                         int recvcount,
                         int MPI_PARAM( root ),
                         MPI_Comm MPI_PARAM( comm ))
{
#ifdef GEOS_USE_MPI
  return MPI_Scatter( sendbuf, sendcount, internal::getMpiType< TS >(),
                      recvbuf, recvcount, internal::getMpiType< TR >(),
                      root, comm );
#else
  static_assert( std::is_same< TS, TR >::value,
                 "MpiWrapper::scatter() for serial run requires send and receive buffers are of the same type" );
  std::size_t const sendBufferSize = sendcount * sizeof(TS);
  std::size_t const recvBufferSize = recvcount * sizeof(TR);
  GEOS_ERROR_IF_NE_MSG( sendBufferSize, recvBufferSize, "size of send buffer and receive buffer are not equal" );
  memcpy( recvbuf, sendbuf, sendBufferSize );
  return 0;
#endif
}
 
template< typename TS, typename TR >
int MpiWrapper::scatterv( TS const * const sendbuf,
                          const int * sendcounts,
                          const int * MPI_PARAM( displs ),
                          TR * const recvbuf,
                          int recvcount,
                          int MPI_PARAM( root ),
                          MPI_Comm MPI_PARAM( comm ))
{
#ifdef GEOS_USE_MPI
  return MPI_Scatterv( sendbuf, sendcounts, displs, internal::getMpiType< TS >(),
                       recvbuf, recvcount, internal::getMpiType< TR >(),
                       root, comm );
#else
  static_assert( std::is_same< TS, TR >::value,
                 "MpiWrapper::scatterv() for serial run requires send and receive buffers are of the same type" );
  std::size_t const sendBufferSize = sendcounts * sizeof(TS);
  std::size_t const recvBufferSize = recvcount * sizeof(TR);
  GEOS_ERROR_IF_NE_MSG( sendBufferSize, recvBufferSize, "size of send buffer and receive buffer are not equal" );
  memcpy( recvbuf, sendbuf, sendBufferSize );
  return 0;
#endif
}
 
 
 
template< typename T >
int MpiWrapper::iRecv( T * const buf,
                       int count,
                       int MPI_PARAM( source ),
                       int tag,
                       MPI_Comm MPI_PARAM( comm ),
                       MPI_Request * MPI_PARAM( request ) )
{
#ifdef GEOS_USE_MPI
  GEOS_ERROR_IF( (*request)!=MPI_REQUEST_NULL,
                 "Attempting to use an MPI_Request that is still in use." );
  return MPI_Irecv( buf, count, internal::getMpiType< T >(), source, tag, comm, request );
#else
  std::map< int, std::pair< int, void * > > & pointerMap = getTagToPointersMap();
  std::map< int, std::pair< int, void * > >::iterator iPointer = pointerMap.find( tag );
 
  if( iPointer==pointerMap.end() )
  {
    pointerMap.insert( {tag, {1, buf} } );
  }
  else
  {
    GEOS_ERROR_IF( iPointer->second.first != 0,
                   "Tag does is assigned, but pointer was not set by iSend." );
    memcpy( buf, iPointer->second.second, count*sizeof(T) );
    pointerMap.erase( iPointer );
  }
  return 0;
#endif
}
 
template< typename T >
int MpiWrapper::recv( array1d< T > & buf,
                      int MPI_PARAM( source ),
                      int tag,
                      MPI_Comm MPI_PARAM( comm ),
                      MPI_Status * MPI_PARAM( request ) )
{
#ifdef GEOS_USE_MPI
  MPI_Status status;
  int count;
  MPI_Probe( source, tag, comm, &status );
  MPI_Get_count( &status, MPI_CHAR, &count );
 
  GEOS_ASSERT_EQ( count % sizeof( T ), 0 );
  buf.resize( count / sizeof( T ) );
 
  return MPI_Recv( reinterpret_cast< char * >( buf.data() ),
                   count,
                   MPI_CHAR,
                   source,
                   tag,
                   comm,
                   request );
#else
  GEOS_ERROR( "Not implemented!" );
  return MPI_SUCCESS;
#endif
}
 
template< typename T >
int MpiWrapper::iSend( arrayView1d< T > const & buf,
                       int MPI_PARAM( dest ),
                       int tag,
                       MPI_Comm MPI_PARAM( comm ),
                       MPI_Request * MPI_PARAM( request ) )
{
#ifdef GEOS_USE_MPI
  GEOS_ERROR_IF( (*request)!=MPI_REQUEST_NULL,
                 "Attempting to use an MPI_Request that is still in use." );
  return MPI_Isend( reinterpret_cast< void const * >( buf.data() ),
                    buf.size() * sizeof( T ),
                    MPI_CHAR,
                    dest,
                    tag,
                    comm,
                    request );
#else
  GEOS_ERROR( "Not implemented." );
  return MPI_SUCCESS;
#endif
}
 
template< typename T >
int MpiWrapper::send( T const * const buf,
                      int count,
                      int dest,
                      int tag,
                      MPI_Comm comm )
{
#ifdef GEOS_USE_MPI
  return MPI_Send( buf, count, internal::getMpiType< T >(), dest, tag, comm );
#else
  GEOS_ERROR( "Not implemented without MPI" );
#endif
}
 
template< typename T >
int MpiWrapper::iSend( T const * const buf,
                       int count,
                       int MPI_PARAM( dest ),
                       int tag,
                       MPI_Comm MPI_PARAM( comm ),
                       MPI_Request * MPI_PARAM( request ) )
{
#ifdef GEOS_USE_MPI
  GEOS_ERROR_IF( (*request)!=MPI_REQUEST_NULL,
                 "Attempting to use an MPI_Request that is still in use." );
  return MPI_Isend( buf, count, internal::getMpiType< T >(), dest, tag, comm, request );
#else
  std::map< int, std::pair< int, void * > > & pointerMap = getTagToPointersMap();
  std::map< int, std::pair< int, void * > >::iterator iPointer = pointerMap.find( tag );
 
  if( iPointer==pointerMap.end() )
  {
    pointerMap.insert( {tag, {0, const_cast< T * >(buf)}
                       } );
  }
  else
  {
    GEOS_ERROR_IF( iPointer->second.first != 1,
                   "Tag does is assigned, but pointer was not set by iRecv." );
    memcpy( iPointer->second.second, buf, count*sizeof(T) );
    pointerMap.erase( iPointer );
  }
  return 0;
#endif
}
 
template< typename U, typename T >
U MpiWrapper::prefixSum( T const value, MPI_Comm comm )
{
  U localResult;
 
#ifdef GEOS_USE_MPI
  U const convertedValue = value;
  int const error = MPI_Exscan( &convertedValue, &localResult, 1, internal::getMpiType< U >(), MPI_SUM, comm );
  MPI_CHECK_ERROR( error );
#endif
  if( commRank() == 0 )
  {
    localResult = 0;
  }
 
  return localResult;
}
 
 
template< typename T >
T MpiWrapper::allReduce( T const & value, Reduction const op, MPI_Comm const comm )
{
  T result;
  allReduce( &value, &result, 1, getMpiOp( op ), comm );
  return result;
}
 
template< typename T >
void MpiWrapper::allReduce( Span< T const > const src, Span< T > const dst, Reduction const op, MPI_Comm const comm )
{
  GEOS_ASSERT_EQ( src.size(), dst.size() );
  allReduce( src.data(), dst.data(), LvArray::integerConversion< int >( src.size() ), getMpiOp( op ), comm );
}
 
template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
void MpiWrapper::allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, int const count, Reduction const op, MPI_Comm const comm )
{
  static_assert( std::is_trivially_copyable< typename get_value_type< SRC_CONTAINER_TYPE >::type >::value,
                 "The type in the source container must be trivially copyable." );
  static_assert( std::is_trivially_copyable< typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
                 "The type in the destination container must be trivially copyable." );
  static_assert( std::is_same< typename get_value_type< SRC_CONTAINER_TYPE >::type,
                               typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
                 "Source and destination containers must have the same value type." );
  GEOS_ASSERT_GE( src.size(), count );
  GEOS_ASSERT_GE( dst.size(), count );
  allReduce( src.data(), dst.data(), count, getMpiOp( op ), comm );
}
 
template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
void MpiWrapper::allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, Reduction const op, MPI_Comm const comm )
{
  static_assert( std::is_trivially_copyable< typename get_value_type< SRC_CONTAINER_TYPE >::type >::value,
                 "The type in the source container must be trivially copyable." );
  static_assert( std::is_trivially_copyable< typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
                 "The type in the destination container must be trivially copyable." );
  static_assert( std::is_same< typename get_value_type< SRC_CONTAINER_TYPE >::type,
                               typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
                 "Source and destination containers must have the same value type." );
  GEOS_ASSERT_EQ( src.size(), dst.size() );
  allReduce( src.data(), dst.data(), LvArray::integerConversion< int >( src.size() ), getMpiOp( op ), comm );
}
 
template< typename T >
T MpiWrapper::sum( T const & value, MPI_Comm comm )
{
  return MpiWrapper::allReduce( value, Reduction::Sum, comm );
}
 
template< typename T >
void MpiWrapper::sum( Span< T const > src, Span< T > dst, MPI_Comm comm )
{
  MpiWrapper::allReduce( src, dst, Reduction::Sum, comm );
}
 
template< typename T >
T MpiWrapper::min( T const & value, MPI_Comm comm )
{
  return MpiWrapper::allReduce( value, Reduction::Min, comm );
}
 
template< typename T >
void MpiWrapper::min( Span< T const > src, Span< T > dst, MPI_Comm comm )
{
  MpiWrapper::allReduce( src, dst, Reduction::Min, comm );
}
 
template< typename T >
T MpiWrapper::max( T const & value, MPI_Comm comm )
{
  return MpiWrapper::allReduce( value, Reduction::Max, comm );
}
 
template< typename T >
void MpiWrapper::max( Span< T const > src, Span< T > dst, MPI_Comm comm )
{
  MpiWrapper::allReduce( src, dst, Reduction::Max, comm );
}
 
 
template< typename T >
T MpiWrapper::reduce( T const & value, Reduction const op, int root, MPI_Comm const comm )
{
  T result;
  reduce( &value, &result, 1, getMpiOp( op ), root, comm );
  return result;
}
 
template< typename T >
void MpiWrapper::reduce( Span< T const > const src, Span< T > const dst, Reduction const op, int root, MPI_Comm const comm )
{
  GEOS_ASSERT_EQ( src.size(), dst.size() );
  reduce( src.data(), dst.data(), LvArray::integerConversion< int >( src.size() ), getMpiOp( op ), root, comm );
}
 
template< typename FIRST, typename SECOND, MpiWrapper::PairReduction const OP >
MpiWrapper::PairType< FIRST, SECOND >
MpiWrapper::allReduce( PairType< FIRST, SECOND > const & localPair, MPI_Comm comm )
{
#ifdef GEOS_USE_MPI
  auto const type = internal::getMpiPairType< FIRST, SECOND >();
  auto const mpiOp = internal::getMpiPairReductionOp< FIRST, SECOND, OP >();
  PairType< FIRST, SECOND > pair{ localPair.first, localPair.second };
  MPI_Allreduce( MPI_IN_PLACE, &pair, 1, type, mpiOp, comm );
  return pair;
#else
  return localPair;
#endif
}
 
template< typename FIRST, typename SECOND, typename CONTAINER, MpiWrapper::PairReduction const OP >
MpiWrapper::PairType< FIRST, SECOND >
MpiWrapper::allReduce( CONTAINER const & pairs, MPI_Comm const comm )
{
  using PAIR_T = PairType< FIRST, SECOND >;
  std::function< PAIR_T( PAIR_T, PAIR_T ) > const getMin = []( PAIR_T const & a, PAIR_T const & b ) {
    return ( std::tie( a.first, a.second ) < std::tie( b.first, b.second ) ) ? a : b;
  };
  std::function< PAIR_T( PAIR_T, PAIR_T ) > const getMax = []( PAIR_T const & a, PAIR_T const & b ) {
    return ( std::tie( a.first, a.second ) > std::tie( b.first, b.second ) ) ? a : b;
  };
  PAIR_T const defaultPair{
    OP == PairReduction::Min ? std::numeric_limits< FIRST >::max() : std::numeric_limits< FIRST >::lowest(),
    OP == PairReduction::Min ? std::numeric_limits< SECOND >::max() : std::numeric_limits< SECOND >::lowest()
  };
  // based on the operation, pair will be the minimum / maximum element (or defaultPair if pairs is empty)
  PAIR_T pair = std::accumulate( pairs.begin(), pairs.end(), defaultPair,
                                 OP == PairReduction::Min ? getMin : getMax );
  return allReduce< FIRST, SECOND, OP >( pair, comm );
}
 
template< typename FIRST, typename SECOND >
MpiWrapper::PairType< FIRST, SECOND > MpiWrapper::min( MpiWrapper::PairType< FIRST, SECOND > const & pair, MPI_Comm comm )
{ return allReduce< FIRST, SECOND, PairReduction::Min >( pair, comm ); }
 
template< typename FIRST, typename SECOND, typename CONTAINER >
MpiWrapper::PairType< FIRST, SECOND > MpiWrapper::min( CONTAINER const & pairs, MPI_Comm comm )
{ return allReduce< FIRST, SECOND, CONTAINER, PairReduction::Min >( pairs, comm ); }
 
template< typename FIRST, typename SECOND >
MpiWrapper::PairType< FIRST, SECOND > MpiWrapper::max( MpiWrapper::PairType< FIRST, SECOND > const & pair, MPI_Comm comm )
{ return allReduce< FIRST, SECOND, PairReduction::Max >( pair, comm ); }
 
template< typename FIRST, typename SECOND, typename CONTAINER >
MpiWrapper::PairType< FIRST, SECOND > MpiWrapper::max( CONTAINER const & pairs, MPI_Comm comm )
{ return allReduce< FIRST, SECOND, CONTAINER, PairReduction::Max >( pairs, comm ); }
 
} /* namespace geos */
 
#endif /* GEOS_COMMON_MPIWRAPPER_HPP_ */