One last MPI example
We have sent individual messages, shared input with broadcast, reduced partial answers, distributed arrays, and exchanged neighboring values. Today we will put a few of those ideas together in a Monte Carlo simulation: using random samples to estimate an answer.
This is the dartboard problem from Chapter 3, programming assignment 3.2 (pp. 155–156), in the second-edition textbook. It also appeared in our older MPI project. We will use it as a class example, then turn to the Exam 1 study guide. Project 2 remains the prime-number census.
Why darts can estimate pi
Imagine a square running from -1 to 1 in both directions, with a circle of radius 1 centered at the origin. The square has area 4 and the circle has area \(\pi\). A uniformly placed dart has about a \(\pi/4\) chance of landing inside the circle. So:
\[ \pi \approx 4\,\frac{\text{darts inside the circle}}{\text{total darts}}. \]
A point \((x,y)\) is inside when \(x^2+y^2\leq1\). We can test that without a square root.
Predict: If 790 of 1,000 darts land inside, what is our estimate? Does a larger number of darts guarantee a smaller error in every individual run?
Decide what moves between processes
| Stage | What every rank does | MPI operation |
|---|---|---|
| Share input | Receive the total number of tosses chosen by rank 0 | MPI_Bcast |
| Divide work | Compute its own number of tosses | Local arithmetic |
| Sample | Generate points and count its own hits | Local loop |
| Combine | Contribute its hit count to rank 0’s total | MPI_Reduce with MPI_SUM |
We do not need to scatter every dart. Each rank can create its own samples. We do not need to gather their coordinates either; only the total hit count matters.
Do all the tosses, including the remainder
long long local_tosses = total_tosses / size;
if (rank < total_tosses % size) local_tosses++;First give every rank the quotient. Then give one extra toss to each of the first few ranks until the remainder is assigned. For ten tosses on four ranks, the counts are 3, 3, 2, 2. For two tosses on four ranks, they are 1, 1, 0, 0. Even a rank with no tosses still participates in the broadcast and reduction.
Unlike scattering an existing array, this distribution needs no MPI_Scatterv: each rank can calculate how much work to generate.
Generate points and count hits
srand(12345u + (unsigned int) rank);
long long local_hits = 0;
for (long long toss = 0; toss < local_tosses; toss++) {
double x = 2.0 * rand() / (RAND_MAX + 1.0) - 1.0;
double y = 2.0 * rand() / (RAND_MAX + 1.0) - 1.0;
if (x * x + y * y <= 1.0) local_hits++;
}rand() produces a pseudorandom integer from 0 through RAND_MAX. Dividing by RAND_MAX + 1.0 produces a value in \([0,1)\); multiplying by 2 and subtracting 1 maps it to \([-1,1)\). We call srand once before sampling, not inside the loop.
The seed starts the sequence. A different seed on each rank avoids having all ranks simply repeat the same sequence. These basic C generators are sufficient for this demonstration; different seeds alone do not prove statistically independent streams. Scientific simulations need more careful random-number choices.
Our fixed seeds make a configuration repeatable with the same C library. Changing the process count changes how the samples are generated, so the estimate can change. That is different from Project 2, where the exact prime count must stay the same.
Reduce counts, then estimate
MPI_Reduce(&local_hits, &total_hits, 1, MPI_LONG_LONG_INT,
MPI_SUM, 0, MPI_COMM_WORLD);Every rank contributes one count, not one value per toss. MPI_LONG_LONG_INT matches C’s long long; the 1 is the number of values contributed by each rank. Rank 0 uses 4.0 * total_hits / total_tosses to calculate a floating-point estimate. Reducing the counts is simpler and safer than averaging local estimates when the ranks have different numbers of tosses.
Complete program
The downloadable mpi_pi.c accepts an optional toss count, with ten million as the default. We use atoll from <stdlib.h> to turn the argument string into a long long, much as we used atoi for an int in Module 5. For this classroom example, enter a whole number from 1 through 1,000,000,000; this simple conversion is not a complete input validator. Rank 0 reads the argument, then broadcasts the value so every rank uses the same toss count. If the count is outside our range, rank 0 prints a message with printf and all ranks finish MPI and exit.
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
int main(int argc, char *argv[]) {
int rank, size;
MPI_Init(NULL, NULL);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
long long total_tosses = 10000000;
if (rank == 0) {
if (argc > 2) {
total_tosses = 0;
} else if (argc == 2) {
total_tosses = atoll(argv[1]);
}
}
MPI_Bcast(&total_tosses, 1, MPI_LONG_LONG_INT, 0, MPI_COMM_WORLD);
if (total_tosses < 1 || total_tosses > 1000000000LL) {
if (rank == 0) printf("Use: ./mpi_pi [tosses from 1 to 1000000000]\n");
MPI_Finalize();
return 1;
}
long long local_tosses = total_tosses / size;
if (rank < total_tosses % size) local_tosses++;
srand(12345u + (unsigned int) rank);
long long local_hits = 0, total_hits = 0;
MPI_Barrier(MPI_COMM_WORLD);
double start = MPI_Wtime();
for (long long toss = 0; toss < local_tosses; toss++) {
double x = 2.0 * rand() / (RAND_MAX + 1.0) - 1.0;
double y = 2.0 * rand() / (RAND_MAX + 1.0) - 1.0;
if (x * x + y * y <= 1.0) local_hits++;
}
MPI_Reduce(&local_hits, &total_hits, 1, MPI_LONG_LONG_INT,
MPI_SUM, 0, MPI_COMM_WORLD);
double elapsed = MPI_Wtime() - start;
double max_elapsed = 0.0;
MPI_Reduce(&elapsed, &max_elapsed, 1, MPI_DOUBLE,
MPI_MAX, 0, MPI_COMM_WORLD);
if (rank == 0) {
double estimate = 4.0 * total_hits / total_tosses;
double error = estimate - 3.141592653589793;
if (error < 0) error = -error;
printf("Tosses: %lld; processes: %d\n", total_tosses, size);
printf("Hits: %lld\n", total_hits);
printf("Pi estimate: %.8f; absolute error: %.8f\n", estimate, error);
printf("Timed calculation: %.6f seconds (maximum rank time)\n", max_elapsed);
}
MPI_Finalize();
return 0;
}mpicc -g -Wall -O2 -o mpi_pi mpi_pi.c
mpirun -n 1 ./mpi_pi 10000000
mpirun -n 2 ./mpi_pi 10000000
mpirun -n 4 ./mpi_pi 10000000Start with 1000 tosses if you want a quick check; then increase the useful work. The output reports the total tosses, hits, estimate, error, and elapsed time. There is no single expected hit count. The estimate should be in \([0,4]\), and the hit count must be between 0 and the toss count.
The program starts its timer after input, work division, and seeding. It times the sampling loop and the reduction of hit counts. Each rank subtracts its own MPI_Wtime readings; a second reduction reports the maximum of those elapsed durations. Printing, startup, and the reduction of timing values are outside the measured region. The initial barrier does not synchronize hardware clocks or guarantee simultaneous starts.
Repeat each configuration at least three times and keep all times. With fixed seeds, repeating the same configuration measures timing variability, not new random samples. For a separate accuracy experiment, change the seed and rebuild. More tosses generally improve the estimate statistically, but an individual larger run can have a larger error. More processes divide the work; they do not automatically improve accuracy when the total number of tosses stays fixed.
Keep MPI_Finalize in the source. The intermittent multi-node finalization problem is separate from this algorithm. The classroom fallback is /usr/bin/mpirun -hosts pc0 -genv HWLOC_COMPONENTS -gl -n 4 ./mpi_pi 10000000. Those are four processes on one board; label that setup when comparing times.
Review together
- Why broadcast the toss count instead of letting every rank prompt for input?
- Why reduce hit counts instead of gathering every dart?
- What would happen if every rank used the same seed and the same number of tosses?
- Why might four processes be slower for a very small toss count?
- Which results should be identical across process counts in the prime project? Which may differ in this simulation?
For next time
Work through the Exam 1 study guide, especially the tracing and debugging exercises. The paper exam is Thursday, October 8. Continue Project 2; it remains due Tuesday, October 6, at 11:59 p.m.