Module 3

MPI: Collective Communication

From messages to collective operations

In Module 2, every process computed part of a sum. Process 0 then received the partial sums one at a time and added them together.

That pattern is common enough that MPI provides a function for it. Today we will combine results with reduction, share input with broadcast, and use those ideas to estimate an area.

This builds on Chapter 3’s trapezoidal-rule example (§3.2), input/output discussion (§3.3), and collective operations (§3.4). We will keep the cyclic distribution from our sum example to make the transition easier. The optional timing section uses §3.6.1.

Create the programs as we go. Complete versions are available here:

A quick check before we start

Run the previous sum program with two and four processes. Both should produce 55.

mpirun -n 2 ./mpi_sum
mpirun -n 4 ./mpi_sum

Which variable holds one process’s contribution? Which process ends up with the complete answer? Does process 0 do any of the arithmetic itself?

Combine the answers with MPI_Reduce

Copy mpi_sum.c to a new file:

cp mpi_sum.c mpi_collective_sum.c

Keep the initialization, n, and the loop that computes local_sum. Replace the entire send/receive section with:

int total = 0;
MPI_Reduce(&local_sum, &total, 1, MPI_INT, MPI_SUM,
           0, MPI_COMM_WORLD);

if (my_rank == 0) {
    printf("Sum from 1 through %d: %d\n", n, total);
}

Keep MPI_Finalize() and return 0; at the end. Add MPI_Barrier(MPI_COMM_WORLD); immediately before finalization on every process. This is the workaround we found for PicoCluster’s intermittent shutdown hang, not a general requirement of MPI. The complete examples below include it, including the invalid-input exit path.

mpicc -g -Wall -o mpi_collective_sum mpi_collective_sum.c
mpirun -n 4 ./mpi_collective_sum

The result should still be 55. We have replaced several explicit messages with one collective operation: an operation involving every process in the communicator.

Argument Meaning
&local_sum This process’s contribution
&total Where the result is stored on the root process
1 Each process contributes one element
MPI_INT The elements are C integers
MPI_SUM Add the contributions
0 Put the result on rank 0, the root
MPI_COMM_WORLD The participating communicator

The count is 1, not comm_sz: each process contributes one integer. Process 0 contributes its own local_sum too.

Every process calls MPI_Reduce, but only rank 0 receives the combined result. The other processes have a variable named total, but this call does not fill it with the answer there.

MPI can choose an efficient communication algorithm, such as combining contributions through a tree. We specify the operation; we do not specify the route each message takes. A reduction is not a promise that this tiny calculation will run faster.

Everyone must participate

Do not put the reduction inside if (my_rank == 0). Only the printing belongs there.

All processes in this communicator must call matching collective operations in the same order. For these examples, they must agree on the root, count, data type, and reduction operation. Collectives do not use the message tags we supplied to MPI_Send and MPI_Recv.

Discuss: What happens if only process 0 calls the reduction? The program is incorrect and may hang: the operation is missing the other contributions.

Read once, then broadcast

So far, changing n has required editing and recompiling. Let us ask for it when the program runs.

Have only process 0 read from the terminal. Then share its value with everyone using MPI_Bcast:

MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);

On rank 0, n supplies the value. On every other rank, the call fills in that process’s own n. Rank 0 calls the broadcast too; we do not write separate receive calls.

Operation Direction Result
Broadcast One process to everyone Everyone has a copy of the root’s value
Reduction Everyone to one process The root has a combined result

Broadcast does not create shared memory. Each process still owns its own variables.

Complete the input-and-sum program

Replace the fixed n with the input and broadcast below. The rest is the same cyclic sum and reduction we just wrote.

#include <stdio.h>
#include <mpi.h>

int main(void) {
    int my_rank, comm_sz;
    MPI_Init(NULL, NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

    int n = 0;
    if (my_rank == 0) {
        printf("How many integers (1-1000000000)? ");
        fflush(stdout);
        if (scanf("%d", &n) != 1 || n < 1 || n > 1000000000) {
            n = 0;
        }
    }
    MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if (n == 0) {
        if (my_rank == 0) {
            printf("Please enter an integer from 1 through 1000000000.\n");
        }
        MPI_Barrier(MPI_COMM_WORLD); /* PicoCluster shutdown workaround. */
        MPI_Finalize();
        return 1;
    }

    long long local_sum = 0;
    for (int i = my_rank + 1; i <= n; i += comm_sz) {
        local_sum += i;
    }

    long long total = 0;
    MPI_Reduce(&local_sum, &total, 1, MPI_LONG_LONG_INT, MPI_SUM,
               0, MPI_COMM_WORLD);

    if (my_rank == 0) {
        printf("Sum from 1 through %d: %lld\n", n, total);
    }
    MPI_Barrier(MPI_COMM_WORLD); /* PicoCluster shutdown workaround. */
    MPI_Finalize();
    return 0;
}

For a noticeable workload, enter 1000000000 (one billion). The expected sum is 500000000500000000. Start with the small checks below, then compare this larger input with one, two, and four processes.

The loop bound n still fits in an int, but the answer does not. This complete version therefore uses long long for local_sum and total, MPI_LONG_LONG_INT for their reduction, and %lld to print the result. Keep those three changes together. The input limit keeps the loop and sum within their data types on PicoCluster. scanf returns the number of successfully read values; we expect one. fflush(stdout) makes the prompt appear before we wait for input.

If input fails or is outside the allowed range, rank 0 sets n to zero. We broadcast that zero before deciding to stop, so everyone takes the same exit path. If rank 0 simply returned while the others waited in a collective, the program could get stuck.

Compile and run:

mpicc -g -Wall -o mpi_collective_sum mpi_collective_sum.c
mpirun -n 4 ./mpi_collective_sum

Enter 10. The result should be:

How many integers (1-1000000000)? 10
Sum from 1 through 10: 55

Try 100 (5050), then 3 with four processes (6). The process with no terms still contributes zero and participates in the reduction. Try 0 to check that the whole program exits cleanly.

Pause and trace: Each process calls broadcast, then reduction. Only rank 0 reads and prints. Why does broadcasting n have to happen before the loop?

A more useful calculation: area under a curve

The book uses the trapezoidal rule to introduce parallel numerical integration. We will estimate the area under \(f(x)=x^2\) from 0 to 1. The exact answer is \(1/3\), so we have something to check.

Divide the interval into n equal pieces of width:

\[h = \frac{b-a}{n}.\]

For each piece, approximate the area by a trapezoid:

\[\text{area of one trapezoid} = h\frac{f(\text{left})+f(\text{right})}{2}.\]

That is width × average height. Add those small areas to estimate the total. You do not need to derive an integration formula to follow the code.

For two trapezoids, the intervals are \([0,0.5]\) and \([0.5,1]\). Their estimated areas are 0.0625 and 0.3125, giving 0.375. More trapezoids bring this example closer to \(1/3\).

Give each process some trapezoids

Number the trapezoids from 0 through n−1. With eight trapezoids and three processes:

Rank Trapezoid indices
0 0, 3, 6
1 1, 4, 7
2 2, 5

This is the same cyclic distribution as before, but our index starts at zero. Every trapezoid is assigned exactly once, even when n is not divisible by the process count.

The book’s first version assigns a consecutive block to each process and assumes an evenly divisible number of trapezoids. Our version assigns individual trapezoids cyclically. It is easier to connect to Tuesday’s loop, although neighboring trapezoids repeat some endpoint evaluations. We will keep that same algorithm when comparing process counts.

Write mpi_trap.c

#include <stdio.h>
#include <mpi.h>

double f(double x) {
    return x * x;
}

int main(void) {
    int my_rank, comm_sz;
    MPI_Init(NULL, NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

    const double a = 0.0;
    const double b = 1.0;
    const int n = 500000000;
    const double h = (b - a) / n;
    double local_area = 0.0;

    for (int i = my_rank; i < n; i += comm_sz) {
        double left = a + i * h;
        double right = a + (i + 1) * h;
        local_area += h * (f(left) + f(right)) / 2.0;
    }

    double total_area = 0.0;
    MPI_Reduce(&local_area, &total_area, 1, MPI_DOUBLE, MPI_SUM,
               0, MPI_COMM_WORLD);

    if (my_rank == 0) {
        printf("Trapezoids: %d\n", n);
        printf("Processes: %d\n", comm_sz);
        printf("Estimated area: %.12f\n", total_area);
        printf("Exact area:     %.12f\n", 1.0 / 3.0);
    }
    MPI_Barrier(MPI_COMM_WORLD); /* PicoCluster shutdown workaround. */
    MPI_Finalize();
    return 0;
}

f is a small C function: give it an x, and it returns x * x. To change the curve later, we can change that function.

Both area examples now use 500 million trapezoids to make the arithmetic take noticeable time. This deliberately does more work than is needed for a good estimate of this simple curve; we are studying runtime as well as correctness. Each process keeps only a few numbers, so increasing n does not allocate an array of that size. Actual runtime depends on the board and compiler.

All processes use the same fixed interval and number of trapezoids, so no broadcast is needed here. Each process computes its own local_area. The reduction adds doubles this time, so the MPI type is MPI_DOUBLE.

mpicc -g -Wall -O2 -o mpi_trap mpi_trap.c
mpirun -n 1 ./mpi_trap
mpirun -n 2 ./mpi_trap
mpirun -n 4 ./mpi_trap

The estimated area should be close to 0.333333333333. Changing the process count should not change the mathematical problem.

Try setting n to 2 and running with one, two, and four processes. All should give 0.375. With four processes, two have no trapezoids and contribute zero.

Restore a larger n afterward. Tiny differences in the last digits across process counts can occur because floating-point addition is rounded, and regrouping additions can change those roundoff errors. Increasing n also does not guarantee that every printed digit improves indefinitely.

Optional: how long did the parallel work take?

Once the calculation works, make a copy named mpi_trap_timed.c. MPI has its own wall-clock timer:

double start = MPI_Wtime();
/* Work to measure goes here. */
double elapsed = MPI_Wtime() - start;

MPI_Wtime() returns a double measured in seconds. We subtract two readings on the same process; clocks on different nodes need not share the same starting point. No extra timer header or feature macro is needed beyond mpi.h.

Before the computation loop, add:

MPI_Barrier(MPI_COMM_WORLD);
double start = MPI_Wtime();

The barrier waits until all processes have entered it before any returns. This keeps a process from starting the timed computation while another has not yet reached the starting point. It does not make their clocks equal or guarantee exactly simultaneous starts.

Immediately after the area reduction, add:

double elapsed = MPI_Wtime() - start;
double max_elapsed = 0.0;
MPI_Reduce(&elapsed, &max_elapsed, 1, MPI_DOUBLE, MPI_MAX,
           0, MPI_COMM_WORLD);

Inside the existing rank-0 printing block, add:

printf("Computation and reduction: %.6f seconds\n", max_elapsed);

Each process measures its local elapsed time. The second reduction selects the largest duration, rather than reporting whichever process happened to finish fastest. This follows the book’s timing approach.

The measured region includes the calculation and the reduction needed to produce the answer. It excludes MPI startup, the initial barrier, printing, the second reduction used to collect timings, and the final shutdown barrier. It is a useful measure of that region, not the total runtime of the launched job.

A broadcast or reduction is not a substitute for a barrier: completing a collective does not generally mean every other process has finished that collective.

The complete timed version includes these changes.

mpicc -g -Wall -O2 -o mpi_trap_timed mpi_trap_timed.c
mpirun -n 1 ./mpi_trap_timed
mpirun -n 2 ./mpi_trap_timed
mpirun -n 4 ./mpi_trap_timed

Start with the supplied n = 500000000. If the one-process run is still too short, try 1000000000; if it is inconveniently slow, try 100000000. Recompile both area programs after changing n, and use the same value for every process count. Aim for a one-process computation lasting a few seconds on your board rather than a particular universal iteration count. Keep the compiler options the same. -O2 enables compiler optimizations; it does not create MPI processes for us.

Repeat each configuration at least three times and record all runs. Compare medians, and note which nodes were used and whether the cluster was busy. The book reports minimum times to approximate a quiet system; here we use medians to describe a typical run on our shared cluster. Keep the same summary method for every process count. More processes can be slower when communication and waiting outweigh the arithmetic saved.

Processes Run 1 (s) Run 2 (s) Run 3 (s) Median (s)
1
2
4

For a first comparison, divide the one-process time by the multiple-process time. For example, 0.8 seconds divided by 0.5 seconds gives a speedup of 1.6. Label this as speedup relative to the one-process MPI version; a tuned serial implementation would be a separate baseline.

Check your understanding

  • Why must rank 0 call MPI_Reduce even though it receives the answer?
  • Why does only rank 0 call scanf, while everyone calls MPI_Bcast?
  • Does broadcast let one process change another process’s variable afterward?
  • What should a process contribute if it has no trapezoids?
  • Why use MPI_SUM for areas but MPI_MAX for elapsed times?
  • Why might a four-process run take longer than a two-process run?

If there is time

  1. Change the curve to 2.0 * x. The area from 0 to 1 is exactly 1 mathematically; the trapezoidal rule is exact for a straight line apart from floating-point rounding.
  2. Move the input-and-broadcast pattern into the area program so rank 0 reads n. Require a positive, bounded value and make sure everyone takes the same path on invalid input.
  3. Look up MPI_Allreduce in §3.4.4. It combines contributions and gives the result to every process. How would that differ from the reduction used here?

Before moving on

Be able to point to the input, each process’s local work, and the collective operation that produces the answer. Next we can distribute arrays with scatter and gather, and compare cyclic and block distributions in more detail.

Back to top