Merge arrays from multiple ranks using MPI

Question

I have code that asks 3 processes to produce 2D vectors that look like this (increments by 0.05 on the second column), with zeroed values (not shown) for the rows that belong to other processes:

Proc 1
100,0,1,1,1
100,0.05,4,1,5
100,0.10,10,5,11
100,...,...,...,...
100,0.30,12,11,15

Proc 2
100,0.35,40,35,55
100,0.40,45,39,55
100,...,...,...,...
100,0.65,180,90,215

Proc 3
100,0.70,130,67,145
100,...,...,...,...
100,1.0,1,1,1

and I am attempting to use MPI_Allreduce to produce a 2D vector of the same size:

100,0,1,1,1
100,0.05,2,2,3
100,0.1,3,2,4
100,0.15,3,2,6
100,0.2,4,2,6
100,0.25,14,10,20
100,0.30,15,11,21
100,0.35,10,6,19
100,0.4,13,9,21
100,0.45,16,12,25
100,0.5,33,17,55
100,0.55,70,33,155
100,0.6,80,30,190
100,0.65,110,45,200
100,0.7,145,134,161
100,0.75,131,127,138
100,0.8,123,120,129
100,0.85,117,114,122
100,0.9,111,111,113
100,0.95,110,108,112
100,1.0,1,1,1

So far I have tried the following:

    // Function that calculates the average number of steps, incrementing over arraysize and probability
    // Output vector has 21 rows and columns represent [Arraysize, Probability, Average, Min Value, Max Value]
    std::vector<std::vector<double>> forest_fire_average_steps(int arraySize, int numberOfRuns, int rank, int numProcs)
    {
        int count = 21 / numProcs;
        int start = rank * count;
        int end = start + count;

        // init probability
        double p;

        // Create empty vector to store results
        std::vector<std::vector<double>> stepsResults(21, std::vector<double>(5, 0));

        if (rank == 0) p = 0;
        if (rank == 1) p = 0.35;
        if (rank == 2) p = 0.7;
        
        if (rank < 3)
        {
            // Iterates over a range of probability values, from 0 to 1 in 0.05 increments.
            for (int i = start; i < end; ++i)
            {
                // Adds the array size and probability value to the first two columns of the row.
                stepsResults[i][0] = arraySize;
                stepsResults[i][1] = p;

                // Runs the forest fire model 'numberOfRuns' times with a defined array size and p
                // and stores the results.
                std::vector<double> runSteps;
                for (int j = 0; j < numberOfRuns; ++j)
                {
                    runSteps.push_back(forest_fire(arraySize, p).stepCount);
                }

                // Sums the runSteps vector.
                // Could have used std::reduce which is more efficient but the HPC G++ compiler was out of date.
                int sum = 0;

                for (auto& n : runSteps)
                {
                    sum += n;
                }

                // Calculates and stores the average of the results.
                double averageSteps = sum / runSteps.size();
                stepsResults[i][2] = averageSteps;
                stepsResults[i][3] = *min_element(runSteps.begin(), runSteps.end());
                stepsResults[i][4] = *max_element(runSteps.begin(), runSteps.end());

                // Increments probability.
                p += 0.05;
            }
        }

        std::vector<std::vector<double>> finalResults(21, std::vector<double>(5, 0));
        for (unsigned int i=0;i<21;++i){
            int ierr = MPI_Allreduce(stepsResults[i].data(), finalResults[i].data(), 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        }

        MPI_Finalize();

        if (rank == 3)
        {
            for (int i = 0; i < 21; i++)
            {
                    for (int j = 0; j < finalResults[i].size(); j++)
                    {
                        if (j < 4)
                        {
                            std::cout << finalResults[i][j] << ",";
                        } else
                        {
                            std::cout << finalResults[i][j];
                        }
                        
                    }

                std::cout << std::endl;
            }
        }

        return finalResults;
    }

I'm not entirely sure why but the output looks like this, where it has clearly skipped some rows (0.6,0.65, etc.) and added extra rows of zeroes:

100,0,1,1,1
100,0.05,2,2,3
100,0.1,3,2,4
100,0.15,3,2,6
100,0.2,4,2,6
100,0.35,10,6,19
100,0.4,13,9,21
100,0.45,16,12,25
100,0.5,33,17,55
100,0.55,70,33,155
100,0.7,145,134,161
100,0.75,131,127,138
100,0.8,123,120,129
100,0.85,117,114,122
100,0.9,111,111,113
0,0,0,0,0
0,0,0,0,0
0,0,0,0,0
0,0,0,0,0
0,0,0,0,0
0,0,0,0,0

Is there something wrong with my implementation that is causing this?

Answer 1

I'm surprised this program does not abort with a segfault. Your problem is in your definition of a matrix as a vector<vector<double>>. MPI wants contiguous buffers, and a vector-of-vectors is a bunch of small arrays, randomly scattered in memory. Write a class for your object and let it store a single vector<double> in which you index with i+j*N or so.