A matrix-vector multiply parallelized using boost::MPI

Question

This program performs a matrix-vector multiply using MPI to split the computation up. It is not extremely robust (for example, it doesn't handle the case where the number of MPI processes does not divide evenly into the dimensions of the matrix), but I mainly wanted to get some practice using MPI. Is my use of the functions correct, and what are some ways I can improve this code?

#include <iostream>
#include <string>
#include <boost/mpi.hpp>
#include <numeric>

void matVecMult( double ** mat, double * vec, double * result, int beg, int end, unsigned N) {
    for (auto i = beg; i < end; i++) {
        result[i] = std::inner_product(mat[i], mat[i] + N, vec, 0.0);
    }
}

int main(int argc, char ** argv) {
    boost::mpi::environment env;
    boost::mpi::communicator world;

    int dim = (argc > 1) ? std::stoi(argv[1]) : 4;
    auto N = static_cast<unsigned>(dim);

    auto **matrix = new double *[N]();
    auto  *vector = new double[N]();
    auto  *result = new double[N]();

    for(unsigned i = 0; i < N; i++){
        matrix[i] = new double[N]();
        for(unsigned j = 0; j < N; j++){
            if(i==j) {matrix[i][j] = i+j;}
        }
        vector[i] = i + i;
    }

    int beg = world.rank() * (N/world.size());
    int end = beg + (N/world.size());

    matVecMult(matrix, vector, result, beg, end, N);

    world.isend(0, world.rank(), result + beg, (N/world.size()));
    if(world.rank() == 0){
        for(int i = 0; i < world.size(); i++) {
            world.recv(i,i, result+(i*N/world.size()), (N/world.size()));
        }
    }

    for(unsigned i = 0; i < N; i++){
        delete[] matrix[i];
    }

    delete[] matrix;
    delete[] vector;
    delete[] result;

    return 0;
}

Answer 1

An obvious way to improve the code is to use standard containers to manage memory instead of raw pointers.

For this code, I would choose std::vector<double> for vector and result, and probably std::vector<std::vector<double>> for matrix (though note that this isn't the most cache-friendly choice for a 2-d matrix). Remember, we can refer to a std::vector's elements as a plain array using the data() member if needed (here, the iterators should be sufficient).

That allows us to eliminate all the delete[] operations, as destructors will take care of that for us - including the non-local exits where new[] throws std::bad_alloc.

---

This loop looks wasteful:

    for(unsigned j = 0; j < N; j++){
        if(i==j) {matrix[i][j] = i+j;}
    }

As the body is conditional on i==j, that's simply equivalent to

        matrix[i][i] = i+i;

---

I don't have experience with MPI, so not commenting on its use. It's probably worth making N be const, and possibly also defining a useful variable:

auto const chunk_size = N/world.size();

Is the variable env needed? If constructing a boost::mpi::environment has some useful side-effect, it may be worth an explanatory comment, as it currently looks like an unused variable (if this is something obvious and expected in MPI, then disregard this comment).