Week 4 · MPI II · Companion to lecture 04b

One File,
Three Questions

Does each rank know its own work? How do the answers come back? Is the answer right? The file prints 41.6667 at every p, and at three of them it is wrong.

Concept map

Monday gave you the vocabulary. Today is one file: Code/fa2026-mpi, branch lecture04, src/main.cpp. Every method is a question first, then the code. Click a node to jump to it.

todaylater
Three questions, in order: does every rank know its own slice, how do p partial answers become one, and is the answer right. It prints 41.6667 every time. That is not the same as being right.

The problem: area under a curve, by trapezoids

Cut [a, b] into n trapezoids, h = (b − a) / n, and

area ≈ h · [ f(x₀)/2 + f(x₁) + … + f(xₙ₋₁) + f(xₙ)/2 ]

// Today: f(x) = x², a = 0, b = 5, n = 1,024,000,000
// Exact: 125/3 = 41.666…

The exact answer is known, the work splits trivially, and n is a dial from "too small to bother" to "worth a cluster". An answer you can check by hand is a gift. Know it before you run anything.

Question 1 · Its own slice

You are rank 4 of 10. Every rank runs this same program. Nobody sends you anything. You know a = 0, b = 5, n = 1,024,000,000, and your own rank. Where does your strip start, where does it end, how many trapezoids are in it?

Rank 4 of 10. Where does your strip start?

h = (b − a)/n first, because everything else is measured in h. The strip is n/p trapezoids wide, so rank 4 starts 4 × (n/p) trapezoids in: myA = a + myRank·(n/p)·h = 2.0. It ends one strip later, myB = 2.5. Count = n/p = 102,400,000.

Nobody had to be told. Rank and p are enough. That is SPMD.

Derive it for any rank


  


  
In the file: firstMethod()


double a, b, h, n, result;        // n is a double
a = A;  b = B;  n = N;            // 0.0, 5.0, 1024000000
h = (b - a) / n;

// What am I working on?
myA = a + myRank * (n / numProcs) * h;
myB = myA + (n / numProcs) * h;
myResult = Trap(myA, myB, n / numProcs, h);

  
Check your answer against the file in order: the constants, then the three lines that matter. No rank needs a message to start: rank, p and four constants are everything. If something about double n bothers you, write it down and hold it. It comes back in "Is it right?".

  







  
Question 2 · Ten partial sums. How do we add them up?

  
Your toolbox is Week 3's: MPI_Send and MPI_Recv. Nothing else exists yet. Count the messages before you write the loop. That number is Monday's C(p).


  

    
How many messages, who sends, and what does rank 0 do with the strip it computed itself?

    

      
      
      
      
    

    
    

      
 p − 1 messages. Everyone but rank 0 sends one double; rank 0 receives p − 1 of them. And rank 0 seeds the total with its own myResult before the loop starts — the most commonly dropped line. No MPI_Probe here, unlike Week 3's strings: rank 0 knows exactly what is coming, one double from each rank.

      
If you said log p: excellent, but do it with Send and Recv only. The tree is Week 5.

    

  


  
In the file: the merge, by hand


if (myRank != 0) {
    MPI_Send(&myResult, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
} else {
    result = myResult;                       // rank 0 is a worker too
    for (int i = 1; i < numProcs; i++) {
        MPI_Recv(&myResult, 1, MPI_DOUBLE, i, 0,
                 MPI_COMM_WORLD, &myStatus);
        result += myResult;
    }
}

  
Week 3's receive loop with a double instead of a string. Nine receives at p = 10, in rank order: Monday's C(p) = p − 1, in code.


  

    
Step the loop

    
The loop order is i = 1 … p−1. Rank 0 waits for rank 1 even if rank 9 was ready first. Each worker's bar is how long its strip took (an illustration, not a measurement); rank 0's row shows which receive it is blocked on.

    

      
      
      
      
    

    
    

  

  
Hold that thought — rank 0 waiting on rank 1 while rank 9's message sits in the buffer — until Week 5. That ordering is the prototype for R03.

  







  
Question 3 · That loop says one thing. Design the call.

  
In words: add one double from every rank, and leave the total on rank 0. You are writing the library function that does it. What must the caller tell you? List the arguments, then count them. There are seven.


  

    
Name the seven

    
Each card is one thing the caller has to say. Decide what it is, then click to see the name the library uses. The seventh usually needs a nudge, and it is the one that matters.

    

    
  


  
In the file: secondMethod()

  

    


// firstMethod: the merge by hand
if (myRank != 0) {
    MPI_Send(...);
} else {
    result = myResult;
    for (i = 1; i < numProcs; i++) {
        MPI_Recv(...);
        result += myResult;
    }
}

    

    


// secondMethod: the same merge, one line
MPI_Reduce(&myResult, &result, 1,
           MPI_DOUBLE, MPI_SUM, 0,
           MPI_COMM_WORLD);

// Same answer. The merge is now the
// library's problem, and it is free
// to use a tree.

    

  

  
The loop is not wrong; it is hand-built. MPI_Reduce says what you want instead of how to get it, and the library picks the algorithm — typically a tree, which is Monday's ⌈log₂ p⌉.


  

    
A student keeps the if (myRank != 0) from firstMethod and puts the MPI_Reduce inside it, so that rank 0 does not call it. What happens?

    

      
      
      
      
    

    
    

      
 Every rank in the communicator calls every collective. No if around it. The root is not exempt: the root argument is how the answer arrives, not who gets to skip the call. The standard calls the wrapped version erroneous; what you observe is garbage at best, a hang at worst.

    

  

  







  
Rank 0 read the input. Now what?

  
Monday's rule: one rank reads, because the launcher decides who gets stdin. So on ranks 1…p−1, a, b, h and n hold garbage, and nobody else can start. Run the call you just designed backwards.


  

    
Reduce takes seven. This one takes five. Which two go away?

    
Click the two arguments that have no job when one rank's values are copied to everyone, then check.

    

    

    

  


  
In the file: thirdMethod()


if (myRank == 0) {
    a = A;  b = B;  n = N;
    h = (b - a) / n;
}
MPI_Bcast(&a, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&b, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&h, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&n, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  

    

      
One to all

      
Reduce, reversed. On rank 0 the buffer is the input; on everyone else it is the output. One root, everyone else receives, and every rank calls it. Root still means root; it just means source now instead of destination.

    

    

      
Before and after

      
Before the broadcasts, a, b, h and n are uninitialized on ranks 1…p−1. After them, every rank holds rank 0's values. Four broadcasts instead of one packed array costs about 24 µs at p = 10 (the deck's model number), against 0.1 s of compute. Right habit, invisible here. Week 6 is the version an array cannot handle.

    

  

  







  
Question 3 · Is it right?

  
Same program, same n, eight values of p. This file's own arithmetic, printed to 12 decimals. At the precision this program normally prints, every one of these reads 41.6667.


  

    
Which answers are wrong?

    

      
pprinted

      
        
141.666666666672

        
241.666666666672

        
341.666666603371

        
441.666666666675

        
641.666666460962

        
741.666666567016

        
841.666666666675

        
1041.666666666668

        
exact41.666666666667

      
    

    

      
      
      
      
    

    
    

      
 p = 3 is wrong from the 8th decimal; 6 and 7 from the 7th. The others agree with 125/3 to the 11th decimal, which is rounding. These are exactly what lecture04 prints; on Cardinal the last digit or two may differ, because MPI_Reduce may add in a different order, but the wrong digits will not move.

    

  


  

    
Second question. What do 1, 2, 4, 8 and 10 have that 3, 6 and 7 do not?

    

      
      
      
      
    

    
    

      
 n/p is a whole number only when p divides 2¹⁶ · 5⁶. At p = 3, 6 and 7 it is not, and the code has a double n and a Trap() that takes whole numbers. The picture below is the bug with numbers small enough to draw.

    

  


  
The bug: 10 trapezoids, 3 ranks

  

    

      
      
    

    
    

  

  
Precisely: the code computes the strip edges with the full fraction, so each strip really is 3⅓ trapezoids wide. Trap() takes a long, so each rank is told 3. It spaces its points h apart from the left edge and uses the true right edge for the last one, so that last trapezoid is 1⅓ h wide and is counted as h. The effect equals dropping a ⅓ sliver per strip: one whole trapezoid, never counted.


  

    
The real file, at any p

    

      
      
    

    

    
Lost trapezoids = n − p · ⌊n/p⌋ = n mod p: the arithmetic the file does. The result line is a closed form of Trap()'s loop on f(x) = x², computed here; it reproduces the wrong digits in the table above and differs from the file only in the trailing rounding digits, where the file's billion-term loop accumulates its own rounding.

  

  
At p = 6 each rank loses ⅔, six ranks lose four trapezoids of 1,024,000,000, and the answer is off by 2 × 10⁻⁷ — forty thousand times the rounding error at p = 1 — printed as 41.6667. A test where everything divides evenly is not a test.

  







  
10 trapezoids, 3 ranks. Nobody gets skipped.

  
New rule: you may only count in whole trapezoids. No fractions, anywhere. Two integers do all of it: how many I own, and how far in I start.


  

    
How many does each rank get: 3, 3, 4 — or 4, 3, 3? Does it matter?

    

      
      
      
      
    

    
    

      
 n / p is the base share, n % p is what is left over, and the leftovers go one apiece to the first n % p ranks. The integer most people miss is the second one: first, the global index of my first trapezoid — myRank · base plus one for each lower rank that got an extra. Then myA = a + first · h and the strip is myN trapezoids wide. Rank 2 computes all of that from its own rank alone.

    

  


  
In the file: fourthMethod()


long n = N;
long base  = n / numProcs;
long extra = n % numProcs;
long myN   = base + (myRank < extra ? 1 : 0);
long first = myRank * base + (myRank < extra ? myRank : extra);

myA = a + first * h;
myB = myA + myN * h;
myResult = Trap(myA, myB, myN, h);

  

    

      
Checked

      
p = 5, 6, 7, 9 all within 2 × 10⁻¹¹ of 125/3. Every one of the 1,024,000,000 trapezoids owned by exactly one rank. This is why every lab checklist says correct at p = 7, not just at powers of two.

    

    

      
The honest alternative

      
Require that p divides n. That is a legitimate design, and it is Pacheco's. But then the program must refuse and MPI_Abort, not truncate silently. Refusing is honest; a quiet wrong answer is not.

    

  


  

    
Explore the split

    

      
      
      
      
    

    

  

  







  
main(): three methods, one stopwatch

  
Every rank times itself. Only rank 0 reports. Nothing fences the region it reports. Count what is missing before you read on.


MPI_Init(NULL, NULL);

timer.startTimer();  firstMethod();   timer.stopTimer();  t1 = ...;
timer.startTimer();  secondMethod();  timer.stopTimer();  t2 = ...;
timer.startTimer();  thirdMethod();   timer.stopTimer();  t3 = ...;

if (myRank == 0) std::cout << "Times: " << t1 << ...;


  

    
Two rules from Week 3 are missing here. One from Monday is too. Which three?

    

      
      
      
      
    

    
    

      
 From Week 3: no barrier before any startTimer, and none before any stopTimer. From Monday: no MAX across ranks — rank 0 reports only its own clock, and without the closing barrier that is not the slowest rank's time. Bonus if you spotted it: std::cout sits inside each timed method.

    

  


  

    
Where the fences go

    
On the branch, main() continues past these lines with fourthMethod timed properly: barrier, MPI_Wtime, Reduce with MPI_MAX. Toggle to see the shape, then scroll to it in the real file.

    

    

  

  







  
Predict, then diagnose: the model says a tie

  
Monday's model: 0.1 s of compute per rank at p = 10 against tens of microseconds of merge. All three methods should land within noise of each other.


  

    
What the model predicts

    

      
    

    

    
Model, not a measurement: compute = (n/p) · t_f with t_f = 1 ns; one double costs α + 8β with α = 2 µs, β = 1/(20 GB/s); the hand merge is p − 1 of those in sequence, Reduce and each Bcast ⌈log₂ p⌉ levels. Real messages differ; the ratios are what matter.

  


  

    
You run it and t₁ (firstMethod) is milliseconds slower than t₂ and t₃. Which explanation does the lecture put on the table?

    

      
      
      
      
    

    
    

      
 If Reduce "wins" by milliseconds, the merge is not why: the model puts all three within tens of microseconds, below the noise of a 0.1 s run. What t₁ carries: no barrier after MPI_Init, so t₁ absorbs every first-call cost (lazy connection setup on the first messages is the likelier one; MPI_Init normally ends with a launcher-level fence, so raw launch skew is small). And rank 0 starts t₂ late: the workers finish firstMethod the moment their small Send is buffered and go straight into secondMethod's compute, while rank 0 is still receiving and printing.

      
These are hypotheses. Don't argue it; test it. ./Default/MPI --twice runs firstMethod again. If the second call is fast, it was a first-call effect. This page does not tell you the answer: that is your run, below.

    

  

  
A gap of milliseconds is not Send/Recv versus Reduce; it is something else, and naming it is the skill. That framing is the Diagnosis section of a reflection in miniature.

  







  
Run it on Cardinal

  
Branch lecture04. It runs all four methods and prints 12 digits; nothing to edit. Predict first: p = 7 is the bug at twelve digits, p = 10 --twice asks whether t₁ is a first-call cost. About 20 minutes.

  
    
    
  1. 
      
    Clone once, check out the lecture branch
    
      
    ssh <your-username>@cardinal.osc.edu
git clone -o class git@gitlab.com:cs4540-spring-2019/cs-4170-5170-fall-2026/mpi.git mpi
cd mpi
git fetch class
git checkout -b lecture04 class/lecture04
    
    
    2. 
    
  2. 
      
    Load the modules, build
    
      
    Modules do not persist between logins: load them every time. The makefile's OSC target links -lboost_mpi even though this file uses no Boost, so the boost module must be loaded or the link fails.
    
      
    module load mvapich/4.1 boost/1.83.0
cd Default && make OSC && cd ..
    
    
    3. 
    
  3. 
      
    Get ten nodes of your own
    
      
    One process per node, so every message crosses the network. The login node is shared; time nothing there.
    
      
    salloc --nodes=10 --ntasks-per-node=1 -c 2 --time=00:15:00 -A PCS0313
module load mvapich/4.1 boost/1.83.0
    
    
    4. 
  


  
Experiments

  

    
A · p = 7: the bug, at twelve digits

    
Write down the digits you expect from firstMethod before you run it. The explorer in "Is it right?" will tell you if you were right.

    
srun -c 2 -n 7 ./Default/MPI

    
What you should see
Methods 1–3 print 41.666666567016 or within a digit or two of it — wrong from the 7th decimal, and still 41.6667 at four. fourthMethod lands within 2 × 10⁻¹¹ of 41.666666666667. Then try -n 6 (four trapezoids of a billion lost) and -n 8 (nothing lost: 8 divides 2¹⁶ · 5⁶).

  

  

    
B · p = 10, --twice: is t₁ a first-call cost?

    
The model says a tie. Before you run, write the explanation you would give if t₁ comes out milliseconds slower.

    
srun -c 2 -n 10 ./Default/MPI --twice

    
What you should see
firstMethod runs twice. Compare the two: if the second call is fast and the first is not, the gap was a first-call effect, not Send/Recv versus Reduce. If both are slow, the hypothesis is wrong and you need another one. The lecture does not tell you which; your run does.

  

  

    
C · One run is an anecdote

    
Five runs, and the minimum. That is what the reflection wants.

    
for i in 1 2 3 4 5; do
    srun -c 2 -n 10 ./Default/MPI --twice | grep -i time
done

    
What you should see
Five "Times" lines. The minimum of each column is your number; the spread between runs is the noise the model said a tens-of-microseconds difference would hide in. Then exit to release the allocation.

  

  
Before Friday (R02): commit your R02 prediction and build the skeleton. R02 — Monte Carlo π — is due Sunday, Sep 20. It divides a fixed 200,000,000 samples by p, and its sweep includes p = 6. Compute π from the samples you actually drew: today's bug is tomorrow's TODO 1. Monday: Reduce and Bcast were two collectives; Week 5 does the rest, and R03 predicts the ratio between a hand-rolled loop and a collective — today's firstMethod versus secondMethod is its prototype.

  







  
Check yourself

  
Six questions, instant feedback, nothing recorded.


  
1. firstMethod, rank 4 of 10, a = 0, b = 5, n = 1,024,000,000. What is myA, and who told rank 4?

    

    
h = 5/1,024,000,000; n/p = 102,400,000; four strips of (n/p)·h = 0.5 each. Rank and p are enough: SPMD.


  
2. In the hand merge, what is the line most people drop, and why does it matter?

    

    
p − 1 messages, and rank 0 seeds the total with its own strip. No Probe needed: one double from each rank, and rank 0 knows it.


  
3. Why does MPI_Reduce need a communicator argument at all, when the program has only one?

    

    
The seventh argument is the one that needs a nudge and the one that matters. Every rank in the communicator calls every collective; the root is how the answer arrives, not who may skip the call.


  
4. Bcast takes five arguments to Reduce's seven. Which two are gone?

    

    
Buffer, count, datatype, root, comm. Root still means root; it now means source instead of destination.


  
5. p = 1, 2, 4, 8, 10 print the right answer and p = 3, 6, 7 do not. The reason:

    

    
At p = 6: four trapezoids of a billion, off by 2 × 10⁻⁷, still printed as 41.6667. A test where everything divides evenly is not a test.


  
6. fourthMethod, n = 10, p = 3. What does rank 2 own?

    

    
Ranks own trapezoids 0–3, 4–6, 7–9. Every trapezoid owned by exactly one rank; every count a whole number.


  

  

  







  
  
← → keys also work