sys/mpiWorld.h

Go to the documentation of this file.

/**
 *  @file      mpiWorld.h
 *  @brief     The MPI implementation that splits a grid among the ranks.
 *  @author    Mikica Kocic
 *  @copyright GNU General Public License (GPLv3).
 */

#include <mpi.h>

#ifdef _STANDALONEMPI
    #include <cstdio>
    #include <cstdlib>
    #include <cmath>
    #define MPIWorld_sanityCheck main
#endif

/** Implements the message passing interface (MPI) for high-performance computing (HPC).
 *
 *  In the MPI environment, the calculations on the grid are split among the ranks
 *  so the edges are exchanged between each two neighbors in the following way:
 *  <pre>
 *        .---- left -----.-------.---- right ----.
 *        | ghost | edge  |       | edge  | ghost |
 *        +-------+-----------------------+-------+
 *   ...  |       |     grid chunk i      |       |
 *        +-------+-----------------------+-------+------------+-------+
 *                                |       |   grid chunk i+1   |       |  ...
 *                                +-------+--------------------+-------+
 *                                :     left      :
 *                                : ghost | edge  :
 *                                `-------v-------´
 *                         exchanged between two neighbors              </pre>
 *
 *  The overlapping edges of one rank become other's ghosts (of the same size).
 *  The ghost of the left-most rank is determined from the inner boundary conditions and
 *  the ghost of the right-most rank is determined from the outer boundary conditions.
 *
 *  When compiling, use `mpic++` with `-D_USEMPI`. To run, use `mpiexec`.
 *
 *  @warning The maximal slicing is not compliant with MPI since the boundary value
 *  problem for the slicing differential equation requires access to the whole grid!
 */
class MPIWorld
{
    int worldSize;         //!< The size of the MPI world
    int myRank;            //!< Our rank in the MPI world
    int leftRank;          //!< The rank of the left neighbor
    int rightRank;         //!< The rank of the right neighbor
    MPI_Request waitLeft;  //!< Outstanding MPI_Isend request to the left rank
    MPI_Request waitRight; //!< Outstanding MPI_Isend request to the right rank

public:

    /** Initializes the MPI environment and finds our rank in the MPI world.
     */
    MPIWorld ()
    {
        // Initialize the MPI environment
        //
        MPI_Init( NULL, NULL );

        // Find out the rank and the size of the MPI
        //
        MPI_Comm_rank( MPI_COMM_WORLD, &myRank );
        MPI_Comm_size( MPI_COMM_WORLD, &worldSize );

        leftRank  = isFirstInRank() ? -1 : myRank - 1;
        rightRank = isLastInRank()  ? -1 : myRank + 1;

        waitLeft = waitRight = MPI_REQUEST_NULL;

        // std::printf( "my rank %d, left %d, right %d\n", myRank, leftRank, rightRank );
    }

    ~MPIWorld ()
    {
         MPI_Finalize ();
    }

    int size () const {
        return worldSize;
    }

    int rank () const {
        return myRank;
    }

    bool isFirstInRank () const {
        return myRank == 0;
    }

    bool isLastInRank () const {
        return myRank == worldSize - 1;
    }

    /** Send our edges (they become other's ghosts) and
     *  receive other's edges (they become our ghosts).
     */
    void exchangeBoundaries
    (
        void* left_ghost,  //!< Where to put the right edge of the left neighbor
        void* left_edge,   //!< Pointer to the left edge of our data
        void* right_edge,  //!< Pointer to the right edge of our data
        void* right_ghost, //!< Where to put the left edge of the right neighbor
        int   edge_size    //!< The edge (or ghost) size
        )
    {
        waitLeft = waitRight = MPI_REQUEST_NULL;

        MPI_Datatype data_type = sizeof(Real) == sizeof(double)
                               ? MPI_DOUBLE : MPI_LONG_DOUBLE;

        // Send our edges (they become other's ghosts)

        if ( leftRank >= 0 ) {
            MPI_Isend( left_edge, edge_size, data_type, leftRank, 0,
                       MPI_COMM_WORLD, &waitLeft );
        }
        if ( rightRank >= 0 ) {
            MPI_Isend( right_edge, edge_size, data_type, rightRank, 0,
                       MPI_COMM_WORLD, &waitRight );
        }

        // Receive other's edges (they become our ghosts)
        //
        if ( leftRank >= 0 ) {
            MPI_Recv( left_ghost, edge_size, data_type, leftRank, 0,
                      MPI_COMM_WORLD, MPI_STATUS_IGNORE );
        }
        if ( rightRank >= 0 ) {
            MPI_Recv( right_ghost, edge_size, data_type, rightRank, 0,
                      MPI_COMM_WORLD, MPI_STATUS_IGNORE );
        }
    }

    /** Wait the other ranks to have received our edges.
     */
    void waitExchange ()
    {
        if ( waitLeft != NULL ) {
            MPI_Status status;
            MPI_Wait( &waitLeft, &status );
        }
        if ( waitRight != NULL ) {
            MPI_Status status;
            MPI_Wait( &waitRight, &status );
        }
    }
};

/** A sanity check test-bed with benchmarking.
 *  To test, replace `main` by MPIWorld_sanityCheck then invoke
 *  `mpiexec -n 4 bim-solver 10 10 10 1`.
 */
int MPIWorld_sanityCheck( int argc, char** argv )
{
    const Int  nGhost = argc >= 2 ? atol ( argv[1] ) : 1;
    const Int  mLen   = argc >= 3 ? atol ( argv[2] ) : 1;
    const Int  nTotal = argc >= 4 ? atol ( argv[2] ) : 80000;
    const bool debug  = argc >= 5;

    MPIWorld mpi;

    Int   nLen        = nTotal / mpi.size(); // nLen should be a multiple of mpi.size()
    Real* gridRow     = new Real[ nLen + 2 * nGhost ];
    Real* myData      = gridRow + nGhost;
    Real* left_ghost  = gridRow;
    Real* left_edge   = gridRow + nGhost;
    Real* right_edge  = gridRow + nGhost + nLen - nGhost;
    Real* right_ghost = gridRow + nGhost + nLen;

    for ( Int m = 0; m <  mLen; ++m )
    {
        left_ghost  [0] = -1;
        left_edge   [0] = ( mpi.rank () + 1 ) * 1000 + m * 100 + 1;
        myData      [0] = ( mpi.rank () + 1 ) * 1000 + m * 100 + 2;
        right_edge  [0] = ( mpi.rank () + 1 ) * 1000 + m * 100 + 3;
        right_ghost [0] = -1;

        if ( debug ) {
            std::printf( "PRE %*s %5.0lf|%5.0lf%5.0lf%5.0lf|%5.0lf\n", mpi.rank() * 16,
                "", *left_ghost, *left_edge, *myData, *right_edge, *right_ghost );
        }

        mpi.exchangeBoundaries(
            left_ghost, left_edge, right_edge, right_ghost, nGhost );

        // Do a sample calculation (equally divided among the ranks)
        // The calculation depends on the left ghost.
        for ( Int n = nGhost; n < nGhost + nLen; ++n ) {
            gridRow[n] = gridRow[n-1] + cos(n) * exp( n * sin(n) );
        }

        mpi.waitExchange ();

        if ( debug ) {
            printf( "POST%*s %5.0lf|%5.0lf%5.0lf%5.0lf|%5.0lf\n", mpi.rank() * 16,
                   "", *left_ghost, *left_edge, *myData, *right_edge, *right_ghost );
        }
    }

    return 0;
}