Step 8 — Verify Correctness and Measure Performance

The solver is now available in five execution modes:

Sequential
OpenMP
MPI
MPI + OpenMP
MPI + CUDA

The final step is to turn those implementations into a defensible numerical and performance study.

There are two separate questions to answer:

Correctness
    Does every implementation solve the same discrete problem?

Performance
    How does the cost change as we add threads, processes, and GPUs?

Keeping these questions separate makes the results easier to interpret. A fast run is useful only after the numerical result has been checked, and a correct solver still needs a careful timing methodology before its scaling can be explained.

1. Start from one numerical reference

The sequential solver from Step 2 is the reference implementation for the parallel versions.

The first assignment asks us to run the sequential program on

10 × 10
20 × 20
40 × 40

grids.

These runs serve two purposes.

First, they confirm that the PCG solver converges on progressively finer grids.

Second, they produce a reference solution against which the OpenMP, MPI, and hybrid implementations can be compared.

The parallel versions should preserve

the same physical domain
the same grid spacing
the same fictitious-domain coefficient
the same finite-difference stencil
the same diagonal preconditioner
the same PCG recurrence

while changing only the execution model.

For the \(40\times40\) validation cases, the comparison chain is therefore

Sequential
    ↓
OpenMP
    ↓
MPI
    ↓
MPI + OpenMP

The MPI+CUDA solver can then be checked against the same numerical formulation on the larger grids used for the GPU experiments.

2. Check convergence before comparing solutions

Every solver reports whether PCG reached the requested tolerance.

The CPU implementation prints messages of the form

[OK] Converged in ... iterations, residual = ...

and records the residual norm periodically in debug.log.

The distributed solver computes its residual norm from

using a local reduction followed by MPI_Allreduce.

The CUDA solver follows the same structure, with the rank-local contribution computed on the GPU before the MPI reduction.

A useful validation record for each run is therefore

grid size
execution mode
process count
thread count
GPU count
iteration count
final residual norm

If two implementations stop at very different iteration counts, that is a signal to inspect the numerical path before drawing conclusions from runtime.

3. Compare floating-point results with a tolerance

Parallel reductions can change the order in which floating-point additions are performed.

For example, the sequential scalar product

is accumulated in one order, while OpenMP may combine several thread-local partial sums and MPI combines rank-local values through a collective reduction.

CUDA adds another reduction stage on the device.

The execution order therefore changes across implementations even when the mathematical expression is the same.

A practical comparison should measure a numerical difference such as

or

and check that it stays within a suitable tolerance.

For two arrays a and b, a simple Python check can be written as

import numpy as np

a = np.loadtxt("reference.csv", delimiter=",")
b = np.loadtxt("candidate.csv", delimiter=",")

diff = a - b

print("L2 difference :", np.linalg.norm(diff))
print("Linf difference:", np.max(np.abs(diff)))

This gives a more useful correctness test than comparing CSV files as text.

4. Preserve each solution before the next run

The current executables generate filenames from the grid size:

solution/solution_M_<M>_N_<N>.csv

The backend and resource configuration are not included in that filename.

As a result, several runs of the same grid can write to the same path.

When collecting validation data, preserve the output immediately after each configuration. For example:

mkdir -p validation

./task_seq 40 40
cp solution/solution_M_40_N_40.csv \
   validation/solution_seq_40_40.csv

export OMP_NUM_THREADS=4
./task_omp 40 40
cp solution/solution_M_40_N_40.csv \
   validation/solution_omp_40_40_t4.csv

mpiexec -np 4 ./task_mpi 40 40
cp solution/solution_M_40_N_40.csv \
   validation/solution_mpi_40_40_p4.csv

The same idea can be used for hybrid and GPU configurations.

This gives every numerical comparison an explicit source file instead of depending on whichever run happened most recently.

5. Visualize the solution field

A scalar error norm is useful for automated checking, while a plot provides a quick way to inspect the overall numerical field.

The repository includes plot_solution.py, which reads a CSV solution and produces both

2D heatmap
3D surface

representations.

The script currently loads

solution/solution_M_800_N_1200.csv

and saves

solution/solution_M_800_N_1200_2D.png
solution/solution_M_800_N_1200_3D.png

at 300 DPI.

The two views answer slightly different questions.

The 2D heatmap makes the spatial structure and domain shape easy to inspect.

The 3D surface makes the magnitude and smoothness of \(u(x,y)\) easier to see.

For the final report, the first assignment explicitly asks for a figure of the solution on the largest reported grid. The repository’s \(800\times1200\) visualization is a natural output for that part of the report.

6. Separate correctness grids from performance grids

The \(40\times40\) problem is useful for implementation checks because it runs quickly and makes it easy to repeat many configurations.

It is too small to tell the whole performance story.

Parallel execution introduces overhead from

OpenMP thread management
MPI halo exchange
MPI global reductions
GPU kernel launches
Host ↔ Device transfers

and those costs can dominate when each worker receives only a small amount of numerical work.

The first assignment therefore distinguishes the small correctness runs from larger reporting grids such as

400 × 600
800 × 1200

The repository follows the same pattern for CPU experiments.

Its Polus workflow also contains a matched comparison stage for

800 × 1200
1600 × 2400
3200 × 4800

using MPI, MPI+OpenMP, one-GPU MPI+CUDA, and two-GPU MPI+CUDA configurations.

This gives us two useful experiment groups:

small grids
    implementation correctness

large grids
    performance and scaling

7. Use the original sequential program as the performance baseline

The coursework defines speedup relative to the original sequential implementation.

Let

be the sequential runtime for a fixed grid, and let

be the runtime of one parallel configuration on the same grid.

Then

is the speedup.

For a CPU configuration using \(p\) parallel workers, efficiency can be reported as

For OpenMP,

For pure MPI,

For a hybrid CPU configuration, the resource description should preserve both

and

If efficiency is normalized by the total CPU execution width, use

and state that convention explicitly in the report.

For GPU configurations, the assignment still asks for speedup relative to the original sequential CPU program. GPU efficiency is less naturally expressed with the same CPU-worker formula, so the resource configuration and speedup should be reported explicitly rather than forcing unlike hardware into one ambiguous efficiency number.

8. Build experiment tables before drawing plots

A useful OpenMP table has the form

The MPI table can use the same structure with process count in the second column.

For hybrid execution, keep both dimensions visible:

For CUDA, keep GPU count visible:

Once the tables are complete, speedup and efficiency plots become mechanical transformations of recorded measurements.

9. Measure more than total runtime

The current CPU executable already divides its runtime into

Initialization
Laplace operator
Update operations
Reduction
Finalization
Total runtime

The MPI executable adds

MPI halo exchange
MPI Allreduce

and the CUDA executable further separates

CUDA Laplace kernel
CUDA update kernels
CUDA reduction kernels
Host → Device copies
Device → Host copies
MPI halo exchange
MPI Allreduce

These categories are especially important for the GPU assignment, which asks for detailed timing of parallel loops, initialization and finalization, memory transfers, and communication.

A useful decomposition is

The exact components depend on the backend, but keeping the categories separate lets us explain why total runtime changes.

10. Be precise about parallel wall-clock time

The current MPI and MPI+CUDA programs print the timing summary from rank 0.

That is useful for inspecting where rank 0 spends time, and the repository labels the MPI summary accordingly.

For a strict parallel wall-clock measurement, the critical path is determined by the slowest participating rank.

There are two practical ways to make the final timing robust:

use the scheduler/job wall-clock measurement

or

reduce per-rank elapsed times with MPI_MAX

For example, an extended implementation could combine a local elapsed time as

double local_time = timer.get("laplace");
double global_time = 0.0;

MPI_Reduce(
    &local_time,
    &global_time,
    1,
    MPI_DOUBLE,
    MPI_MAX,
    0,
    cart_comm
);

This reports the longest rank time to rank 0.

The current repository does not perform this MPI_MAX aggregation in its printed solver summaries, so the timing source used in the final tables should be stated clearly.

11. Explain scaling through work and overhead

For a fixed global grid, adding CPU workers reduces the local computational work available to each worker.

At the same time, some overheads remain.

For OpenMP, the main effects include

parallel-region overhead
thread scheduling
shared-memory bandwidth
reduction synchronization

For MPI, the local grid becomes smaller as the process count grows, while each PCG iteration still requires

halo communication
global reductions

A useful conceptual model is

For the hybrid solver,

depends on both the number of MPI subdomains and the number of threads working inside each subdomain.

Two configurations with the same nominal CPU width can therefore behave differently:

4 MPI × 1 thread
2 MPI × 2 threads
1 MPI × 4 threads

because they create different numbers of MPI boundaries and different amounts of shared-memory parallelism.

This is exactly why process count and thread count should remain separate in the result tables.

12. Explain CUDA performance through kernels, transfers, and MPI

The GPU version introduces another layer.

Inside each rank, the main PCG operations execute on the GPU:

stencil
vector updates
preconditioner
local reductions

The current halo implementation stages boundary values through host memory:

GPU boundary
    ↓
Device → Host
    ↓
MPI exchange
    ↓
Host → Device
    ↓
GPU halo

So a useful model for one distributed GPU iteration is

This explains why GPU performance should be evaluated from total runtime and the timing breakdown together.

The assignment specifically asks the MPI+GPU implementation to be compared with the MPI/OpenMP implementation and asks unexpected performance behavior to be explained.

If a GPU kernel is fast while total runtime improves only modestly, the transfer and communication categories provide the first place to look.

If a larger grid improves GPU speedup, the likely reason should be tested against the timing data: larger problems provide more local numerical work relative to fixed launch and transfer overheads.

The report should use the measured timing categories to support such explanations instead of treating them as assumptions.

13. Use the hardware numbers as context, not as measured performance

The GPU assignment provides approximate hardware figures for the Polus CPUs and Tesla P100 GPUs.

It lists double-precision peak throughput of roughly

POWER8 CPU    0.3 TFLOP/s
Tesla P100    4.7 TFLOP/s

and memory bandwidth of roughly

POWER8 CPU    230 GB/s
Tesla P100    700 GB/s

These numbers help explain why the GPU can offer much more raw arithmetic throughput and memory bandwidth.

They are hardware reference values.

The Poisson-PCG solver also contains communication, reductions, halo copies, and synchronization, so application speedup should come from measured program runtime.

A useful report can therefore present the hardware figures as context and the solver timings as the actual evidence.

14. Build the final comparison around matched configurations

The repository’s Polus stage 9 is designed for matched large-grid comparisons.

For each of

800 × 1200
1600 × 2400
3200 × 4800

it contains LSF configurations for

MPI
MPI + OpenMP
MPI + CUDA with 1 GPU
MPI + CUDA with 2 GPUs

This is a good structure for the final cross-backend table because every row uses the same numerical grid.

A comparison sheet could therefore look like

Repeat the same structure for the larger grids.

The values should come from the actual Polus logs rather than being filled from expectations.

15. Assemble the report as an argument

A strong final report should connect the entire development path.

A practical structure is

1. Problem setup
   Poisson equation
   fictitious domain
   grid and boundary conditions

2. Numerical method
   finite-difference operator
   diagonal preconditioner
   PCG iteration

3. Sequential reference
   implementation
   convergence
   baseline timing

4. OpenMP
   parallel loops
   reductions
   correctness and scaling

5. MPI
   2D decomposition
   halo exchange
   MPI_Allreduce
   correctness and scaling

6. MPI + OpenMP
   hybrid execution
   process/thread configurations

7. MPI + CUDA
   GPU kernels
   device memory
   GPU reductions
   staged halo exchange

8. Results
   solution visualization
   timing tables
   speedup and efficiency plots
   operator timing breakdown

9. Discussion
   scaling limits
   communication cost
   transfer cost
   unexpected behavior

The first assignment explicitly asks the report to include the numerical method, OpenMP/MPI/hybrid implementations, experiment results, the largest-grid solution figure, and speedup plots.

The GPU extension adds the requirement to compare MPI+GPU against the earlier CPU implementations, explain correctness assessment, report detailed timing components, and discuss unexpected performance behavior.

The report therefore becomes the final connection between the mathematics, the implementation, and the measured execution behavior.

16. What have we built?

The complete project now follows one numerical problem through a sequence of execution models:

Poisson equation
      ↓
finite-difference system
      ↓
PCG solver
      ↓
sequential reference
      ↓
OpenMP
      ↓
MPI
      ↓
MPI + OpenMP
      ↓
MPI + CUDA
      ↓
IBM Polus experiments
      ↓
correctness and performance analysis

Across these steps, the numerical core remains recognizable.

The changes happen around data ownership and execution:

OpenMP
    shares one grid among CPU threads

MPI
    distributes the grid across processes

MPI + OpenMP
    combines distributed subdomains with local CPU threads

MPI + CUDA
    keeps the MPI decomposition and executes local PCG kernels on GPUs

The final evaluation closes the loop by checking that every implementation still solves the same problem and by measuring the cost of each parallel execution model.

At that point, the SM25 repository is more than a collection of parallel versions of the same solver. It provides a complete path from the mathematical problem to a reproducible supercomputing experiment.