Reconstructing Matrix scattered over many MPI ranks

Question

I have a program that does some matrix computations using MPI (MPICH). Each rank has a slice of the matrix and does the computations on their slice to get a new slice of the matrix. Sometimes I need to synchronize the matrix on the manager for snapshotting, output, etc. so I need all ranks to send their current matrix slice back to the manager.

I've come up with a solution to do this, and it seems to work OK; but I wanted to know whether this is a good approach or if I'm doing something stupid since I've never worked with MPI before. One thing I'm particularly curious about is whether I'm deleting my requests vector correctly. Thanks in advance!

Here's the code:

#define __WORLD MPI_COMM_WORLD

// reconstructs the matrix on the manager node
void reconstructMatrix(double** current, const int COLS, const int ROWS,
                       const int MY_RANK, const int MPI_SIZE)
{
    auto requests = new std::vector<MPI_Request*>();
    requests->reserve(COLS);

    // Manager receives data from all ranks
    if (MY_RANK == 0)
    {
        // number of columns assigned to all ranks
        const auto SLICE_SIZE = COLS / MPI_SIZE;

        // receive from rank 1 to rank n-1
        for (auto src = 1; src < MPI_SIZE; src++)
        {
            for (auto i = (SLICE_SIZE * src); (i < COLS) && (i < SLICE_SIZE * (src + 1)); i++)
            {
                auto temp = new MPI_Request();
                MPI_Irecv(current[i], ROWS, MPI_DOUBLE, src, 1, __WORLD, temp);
                requests->push_back(temp);
            }
        }
    }
    else // Workers send data to manager  
    {
        // index 0 and COLS + 1 are ghost columns
        for (auto i = 1; i <= COLS; i++)
        {
            auto temp = new MPI_Request();
            MPI_Isend(current[i], ROWS, MPI_DOUBLE, 0, 1, __WORLD, temp);
            requests->push_back(temp);
        }
    }
    // wait for sync to finish
    waitAll(requests);

    requests->clear();
    delete requests;
    requests = nullptr;
}

// waits for all MPI requests and deletes the vector
void waitAll(std::vector<MPI_Request*>*& requests)
{
    for (auto& req : *requests)
    {
        MPI_Status stat;
        MPI_Wait(req, &stat);
        if (stat.MPI_ERROR != MPI_SUCCESS)
        {
            std::cerr << "MPI Request failed!\n\tCode: " << stat.MPI_ERROR
                      << "\n\tFrom: " << stat.MPI_SOURCE;
        }
        delete req;
    }
}

Answer 1

Note: I've written a matrix multiplication benchmark program in C and MPI:https://github.com/thefangbear/matrix-mpi. It uses a row-major format to store the two matrices in 1D arrays for better cache coherency. It also supports both synchronous and asynchronous send/receive modes, tunable via a preprocessor macro, so you can take a look.

With that being said, let's turn to a few improvable parts in your code:

1. You don't have to use an array of pointers to represent the matrix format. You can use an 1D array with row-major/column-major format. If you use a "2D array" (array of pointers), each row/colum slice of your matrix are not guaranteed to be stored in consecutive locations.

2. Isend/Irecv lead to degraded performance. You can verify this via benchmarking. MPI sends/receives extra messages when asynchronous mode is used, which makes basic operations more expensive. And since you have blocking anyways somewhere in your code, ISend/IRecv do not make a huge difference.

3. It is common for worker nodes to send back processed slices of data to be merged and outputted in the master node. You can technically just invoke printf in a loop on each worker node to output its own slice of data (since mpirun redirects all stdout streams on each node to your terminal screen), but MPI does not guarantee the order of output.

4. If you want a more elegant semantics, you should checkout MPI_Scatterv/MPI_Gatherv scatter/gather functions: https://www.mpich.org/static/docs/v3.1/www3/MPI_Scatterv.html. (But I don't think you'll get significant performance gains by using them)