The Trapezoid,
Three Ways, in Threads

Concept map

Monday gave you the words. Today is one file that uses all of them. Click a node to jump to it.

Same problem. What changed?

f(x) = x², [0, 5], n = 1.024·10⁸ today (the Cardinal job uses Week 4's 1.024·10⁹). Exact answer 125/3 = 41.666…

TrapRule1: split by hand, critical inside

As written in 2018. main opens #pragma omp parallel and makes one call. Every thread runs the whole function.

void parallelTrap(double a, double b, int n, double& result){
    double h, x, myRetValue;
    int myId       = omp_get_thread_num();
    int numThreads = omp_get_num_threads();
    int myN, myA, myB;

    h   = (b - a) / n;
    myN = n / numThreads;
    myA = a + myId * myN * h;
    myB = myA + myN * h;

    myRetValue = (f(myA) + f(myB)) / 2.0;
    for (int i = 1; i < myN; i++) {
        x = myA + i * h;
        myRetValue += f(x);
    }
    myRetValue = myRetValue * h;

    #pragma omp critical
    result += myRetValue;
}

[0, 5], n = 1.024e8                      exact 41.667
t=1   41.666666666671155   CHECK ok
t=2   32.916666507978192   CHECK WRONG
t=3   32.407406631586333   CHECK WRONG
t=4   29.479166538494464   CHECK WRONG
t=5   41.666666666666224   CHECK ok

The old input: [0, 16], n = 16          serial 1368
t=1, 2, 4, 8   1368        CHECK ok
t=3, 5         1127.5      CHECK WRONG

The bug, explored

The same arithmetic TrapRule1 does, computed here. Move t and watch where each thread thinks its strip starts.

TrapRule2: the caller merges

The worker now returns its strip. It no longer touches result, so it can be tested alone. Main owns the merge.

// parallelTrap() now RETURNS its strip (same int split)
#pragma omp parallel num_threads(numThreads)
{
    //Incorrect First option
    // #pragma omp critical
    // result += parallelTrap(a, b, n);

    double localResult;
    localResult += parallelTrap(a, b, n);

    #pragma omp critical
    result += localResult;
}

Also note num_threads(numThreads) with numThreads = 4 hard-coded in main: Monday's backup bug. The port drops it so a sweep over t can work.

$ make all          (g++ 13 -O2)
warning: 'localResult' is used uninitialized
t=1..7   same digits as the TrapRule1 port    (it happened to be 0)

$ clang++ -O2       (same source, [0,16], n=16)
t=1   trap2   nan
t=2   trap2   nan
t=4   trap2   nan

trap2_CriticalCall   t=1..7   CHECK ok

What the critical holds

A critical section is correct for any block of code. The question is only what you put inside it.

The fix is Week 4's fix

fourthMethod with thread ids. Counts in long, coordinates in double, the first n % t threads take one extra.

// myStrip(): used by trap1, trap2, trap2_Reduction, trap2_CriticalCall
int  tid = omp_get_thread_num();
int  t   = omp_get_num_threads();

double h = (b - a) / n;

long base  = n / t;
long extra = n % t;

long myN   = base;
long first = tid * base + extra;

if (tid < extra) {                 // the first n % t threads take one more
    myN   = base + 1;
    first = tid * base + tid;
}

double myA = a + first * h;
double myB = a + (first + myN) * h;
// … the serial trapezoid over [myA, myB] with myN steps, times h …

// trap2: TrapRule2's shape, fixed
#pragma omp parallel
{
    double local = myStrip(a, b, n);  // declared AND set

    #pragma omp critical
    result += local;                  // t updates, not n
}

// trap2_Reduction: the merge as a clause
#pragma omp parallel reduction(+:result)
result += myStrip(a, b, n);
// reduction works on a bare region too:
// a private result per thread, added at the join

TrapRule3: the serial loop, plus one line

Put this next to TrapRule1. Everything the programmer computed by hand is now the runtime's job.

double h        = (b - a) / n;
double retValue = (f(a) + f(b)) / 2.0;

#pragma omp parallel for reduction(+:retValue)
for (long i = 1; i <= n - 1; i++) {
    retValue += f(a + i * h);
}
retValue = h * retValue;

// The 2019 file still calls this serialTrap().
// The name is the only thing left over.

Right to 10⁻¹³, and different in the last digits

Every correct version, t = 1, 2, 3, 4, 7. All CHECK ok against 125/3.

           trap3                 trap1 (manual split)
t=1   41.666666666671155    41.666666666671155
t=2   41.666666666671055    41.666666666671048
t=3   41.666666666667766    41.666666666668483
t=4   41.666666666669542    41.666666666669542
t=7   41.666666666668064    41.666666666668419
exact 41.666666666666667

trap3, ten runs at t=4:  ...542 ×8   ...549 ×2
trap3, ten runs at t=7:  ...064 ×9   ...071 ×1

Run it on Cardinal

Same branch as Monday: lecture09. If you built it then, skip to the experiments. 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/openmp.git openmp
   cd openmp
   git fetch class
   git checkout -b lecture09 class/fa2026-openmp/lecture09

   copy

2. Build, then get cores of your own

   One warning is expected ('h' may be used uninitialized — Monday's bug). The login node is shared; anything you time goes in an allocation.

   module load gcc/12.3.0
   cd Default && make OSC && cd ..
   salloc -N 1 -n 1 -c 8 -t 00:15:00 -A PCS0313

   copy

   Every run: OMP_NUM_THREADS=t srun -c 8 ./Default/OMP <method> [n [a b]]. The trapezoid methods: trap_Serial trap2_CriticalCall trap1 trap2 trap2_Reduction trap3. Default n = 1.024·10⁸.

Experiments

for m in trap_Serial trap1 trap2 trap2_Reduction trap3; do
    OMP_NUM_THREADS=3 srun -c 8 ./Default/OMP $m | grep -E '^trap|CHECK'
done

for t in 1 8; do
    OMP_NUM_THREADS=$t srun -c 8 ./Default/OMP trap2_CriticalCall | grep -E 'CHECK|ELAPSED'
    OMP_NUM_THREADS=$t srun -c 8 ./Default/OMP trap1              | grep -E 'CHECK|ELAPSED'
done

cd Default && make OSC && cd ..
OMP_NUM_THREADS=2 srun -c 8 ./Default/OMP trap1
OMP_NUM_THREADS=5 srun -c 8 ./Default/OMP trap1
OMP_NUM_THREADS=3 srun -c 8 ./Default/OMP trap1 16 0 16
OMP_NUM_THREADS=4 srun -c 8 ./Default/OMP trap1 16 0 16
git checkout src/main.cpp && cd Default && make OSC && cd ..    # put it back