Step 5 — Combine MPI and OpenMP

The MPI solver from Step 4 divides the global grid across processes. Each process owns one rectangular subdomain, exchanges halo values with its neighbors, and participates in global reductions.

The next stage adds a second level of parallelism inside each MPI process.

The decomposition remains unchanged:

global grid
    ↓
MPI subdomains
    ↓
one subdomain per process

Inside each subdomain, OpenMP threads now cooperate on the local numerical work:

MPI rank
    ↓
local subdomain
    ↓
OpenMP threads

This gives us a hybrid execution model:

The numerical method does not need another redesign. The main task is to combine the communication structure from Step 4 with the shared-memory loop parallelism from Step 3.

1. Separate the responsibilities of MPI and OpenMP

The hybrid solver is easiest to understand when the two models have clearly separated responsibilities.

MPI manages work that crosses process boundaries:

MPI
├── divide the global grid
├── identify neighboring subdomains
├── exchange halo values
├── combine distributed scalar products
└── gather the final solution

OpenMP manages work that stays inside one process:

OpenMP
├── initialize local grid values
├── apply the local stencil
├── update PCG vectors
├── apply the local preconditioner
├── compute local reductions
└── pack and unpack local buffers

The result is a two-level hierarchy:

global problem
│
├── MPI rank 0
│   ├── thread 0
│   ├── thread 1
│   ├── thread 2
│   └── thread 3
│
├── MPI rank 1
│   ├── thread 0
│   ├── thread 1
│   ├── thread 2
│   └── thread 3
│
└── ...

Each MPI rank still owns one local subdomain. OpenMP changes how the work inside that subdomain is executed.

2. Reuse the MPI solver

A useful property of the current implementation is that the MPI-only and MPI+OpenMP versions use the same solver code.

MPIPoissonSolver already inherits the local numerical operations from PoissonSolver, including loops such as

#pragma omp parallel for schedule(static) collapse(2)
for (int i = 1; i < M; ++i) {
    for (int j = 1; j < N; ++j) {
        ...
    }
}

These directives appear in the local stencil, vector updates, initialization, preconditioning, and local reductions.

When the MPI program is compiled without OpenMP support, those loops execute sequentially inside each process.

When the same sources are compiled with

-fopenmp

the directives become active and every MPI process creates an OpenMP team for the local loops.

This means the transition from MPI to hybrid execution is mostly a build and runtime configuration change:

same domain decomposition
same halo exchange
same MPI_Allreduce calls
same PCG control flow
        ↓
OpenMP enabled for local work

3. Parallelize the local stencil inside each rank

Consider the distributed apply_A() operation.

MPI has already ensured that the halo of p is valid:

exchange_halo(p);

After that exchange, each rank has all values needed to evaluate its local five-point stencil.

The local operation can therefore use the same OpenMP loop from Step 3:

#pragma omp parallel for schedule(static) collapse(2)
for (int i = 1; i < M; ++i) {
    for (int j = 1; j < N; ++j) {
        double kx_plus  = 0.5 * (k[i][j] + k[i+1][j]);
        double kx_minus = 0.5 * (k[i][j] + k[i-1][j]);
        double ky_plus  = 0.5 * (k[i][j] + k[i][j+1]);
        double ky_minus = 0.5 * (k[i][j] + k[i][j-1]);

        out[i][j] = -(
            (
                kx_plus  * (v[i+1][j] - v[i][j])
              - kx_minus * (v[i][j] - v[i-1][j])
            ) / hx2
            +
            (
                ky_plus  * (v[i][j+1] - v[i][j])
              - ky_minus * (v[i][j] - v[i][j-1])
            ) / hy2
        );
    }
}

The sequence is now:

MPI halo exchange
        ↓
halo values are ready
        ↓
OpenMP parallel stencil

MPI provides the remote values. OpenMP distributes the local stencil work among threads.

4. Keep pointwise PCG updates local

The vector updates do not require communication between ranks because every process updates only the entries that belong to its own subdomain.

For example,

is evaluated locally as

#pragma omp parallel for schedule(static) collapse(2)
for (int i = 1; i < M; ++i)
    for (int j = 1; j < N; ++j)
        u[i][j] += alpha * p[i][j];

The same applies to

and

So these operations follow a simple hierarchy:

one MPI rank
      ↓
one local array
      ↓
OpenMP parallel loop

No MPI communication is needed until a later operation depends on data from another subdomain.

5. Combine OpenMP reduction with MPI_Allreduce

The scalar products are where the two parallel models interact most clearly.

Take

In the MPI-only version, each process computes a local sum and then uses MPI_Allreduce:

local grid
    ↓
local sum
    ↓
MPI_Allreduce
    ↓
global sum

In the hybrid version, the local sum itself is parallelized with OpenMP:

double local_rz = 0.0;

#pragma omp parallel for reduction(+:local_rz) \
    collapse(2) schedule(static)
for (int i = 1; i < M; ++i) {
    for (int j = 1; j < N; ++j) {
        local_rz += r[i][j] * z[i][j];
    }
}

After the OpenMP region finishes, one local scalar remains:

double global_rz = 0.0;

MPI_Allreduce(
    &local_rz,
    &global_rz,
    1,
    MPI_DOUBLE,
    MPI_SUM,
    cart_comm
);

The complete reduction is therefore hierarchical:

thread 0 ─┐
thread 1 ─┼─ OpenMP reduction ─► local scalar ─┐
thread 2 ─┤                                  │
thread 3 ─┘                                  │
                                             ├─ MPI_Allreduce
other MPI ranks ─────────────────────────────┘
                                             ↓
                                        global scalar

The same pattern is used for

and

This is one of the defining structures of the hybrid solver:

6. Keep MPI communication outside the OpenMP loops

The current solver uses OpenMP to parallelize local memory operations while the MPI communication calls remain outside those parallel loops.

For example, halo exchange follows this pattern:

OpenMP pack boundary buffer
        ↓
MPI_Irecv / MPI_Isend
        ↓
MPI_Waitall
        ↓
OpenMP unpack received buffer

The boundary packing can use multiple threads:

#pragma omp parallel for schedule(static)
for (int j = 0; j <= N; ++j)
    send_left[j] = field[1][j];

Then the process posts the MPI communication:

MPI_Irecv(recv_left.data(), N + 1, MPI_DOUBLE, left_rank, 0, cart_comm, &reqs[rq++]);
MPI_Isend(send_left.data(), N + 1, MPI_DOUBLE, left_rank, 1, cart_comm, &reqs[rq++]);

with the same request bookkeeping (reqs and rq) introduced in Step 4.

After MPI_Waitall, OpenMP can again help copy the received values into the halo:

#pragma omp parallel for schedule(static)
for (int j = 0; j <= N; ++j)
    field[0][j] = recv_left[j];

This keeps the execution model simple: local array work is threaded, while the communication calls are issued in the normal process-level control flow.

The same principle applies to global reductions. Threads first produce the rank-local scalar, and MPI_Allreduce is called after the OpenMP reduction has completed.

7. Understand the hybrid PCG iteration

Putting the two levels together gives the following iteration:

MPI halo exchange of p
        ↓
OpenMP parallel apply_A
        ↓
OpenMP local reduction for (p, Ap)
        ↓
MPI_Allreduce
        ↓
compute global alpha
        ↓
OpenMP update u
OpenMP update r
        ↓
OpenMP local residual reduction
        ↓
MPI_Allreduce
        ↓
global residual norm
        ↓
OpenMP apply preconditioner
        ↓
OpenMP local reduction for (r, z)
        ↓
MPI_Allreduce
        ↓
compute global beta
        ↓
OpenMP update p

The algorithm still follows the same PCG recurrence used in every previous step.

What changes is the execution hierarchy:

MPI handles global dependencies
OpenMP handles local parallel work

This lets us reuse both the distributed-memory structure from Step 4 and the shared-memory kernels from Step 3.

8. Build the hybrid executable

The MPI-only executable was compiled with

mpicxx -O3 -o task_mpi \
    task_mpi.cpp \
    src/conjugate_gradient.cpp \
    src/mpi_conjugate_gradient.cpp \
    -Iinclude \
    -lm \
    -std=c++11

To activate the OpenMP regions, add -fopenmp:

mpicxx -O3 -o task_mpi_omp \
    task_mpi.cpp \
    src/conjugate_gradient.cpp \
    src/mpi_conjugate_gradient.cpp \
    -Iinclude \
    -lm \
    -std=c++11 \
    -fopenmp

The number of MPI processes is selected with

mpirun -np <processes>

and the number of OpenMP threads inside each process is controlled with

export OMP_NUM_THREADS=<threads>

For example:

export OMP_NUM_THREADS=4
mpirun -np 2 ./task_mpi_omp 40 40

This configuration creates two MPI processes, each using four OpenMP threads.

Conceptually:

2 MPI ranks × 4 OpenMP threads
                =
8 OpenMP worker threads across the two ranks

The important distinction is that these eight threads do not all share the same global grid. Four belong to one MPI subdomain and four belong to the other.

9. Run the required hybrid checks

The assignment asks us to test the \(40\times40\) problem with four OpenMP threads per MPI process using

1 MPI process × 4 OpenMP threads
2 MPI processes × 4 OpenMP threads

The repository automates these cases with

export OMP_NUM_THREADS=4

mpirun -np 1 ./task_mpi_omp 40 40
mpirun -np 2 ./task_mpi_omp 40 40

The first configuration is useful because it exercises the hybrid binary with only one MPI rank. The second introduces communication between two subdomains while retaining four local threads per rank.

For both runs, compare the numerical solution with the sequential, OpenMP, and MPI reference results.

The solver should still represent the same discretized Poisson problem and converge under the same numerical settings. Use a numerical tolerance rather than requiring bit-identical values: the OpenMP reductions can combine floating-point sums in a thread-dependent order, so the hybrid result need not match the other implementations in the last bits.

10. Explore larger hybrid configurations

Once the \(40\times40\) cases are working, larger grids can be used to explore different combinations of processes and threads.

The repository includes examples such as

export OMP_NUM_THREADS=1
mpirun -np 2 ./task_mpi_omp 400 600

export OMP_NUM_THREADS=2
mpirun -np 2 ./task_mpi_omp 400 600

export OMP_NUM_THREADS=4
mpirun -np 2 ./task_mpi_omp 400 600

export OMP_NUM_THREADS=8
mpirun -np 2 ./task_mpi_omp 400 600

and analogous runs on the larger \(800\times1200\) grid.

These experiments let us vary two quantities independently:

The total amount of CPU parallelism is approximately

while the number of distributed subdomains is still only \(P\).

This distinction matters because changing \(P\) changes the communication pattern, while changing \(T\) changes how much shared-memory parallelism is used inside each subdomain.

11. Think about processes and threads separately

Two configurations can use the same total number of CPU workers and still behave differently.

For example:

4 MPI processes × 1 thread
2 MPI processes × 2 threads
1 MPI process   × 4 threads

all use four CPU workers in total.

Their execution structures are different:

4 × 1
many subdomains
more MPI boundaries
little shared-memory parallelism

2 × 2
fewer subdomains
two threads per local domain

1 × 4
one subdomain
no inter-rank halo exchange
four shared-memory threads

This is the main reason to study the hybrid version separately from pure MPI and pure OpenMP.

The best balance depends on the grid size, machine topology, communication cost, and the amount of local work available to each thread. We will compare these effects later using measured runtimes instead of assuming that one configuration is always preferable.

12. What did hybrid parallelism add?

The numerical solver has now passed through three CPU execution models.

Sequential:

one process
one thread
one complete grid

OpenMP:

one process
many threads
one complete grid

MPI:

many processes
one local subdomain per process
communication between subdomains

MPI+OpenMP:

many processes
one local subdomain per process
many threads inside each process

The hybrid version preserves the two communication patterns introduced by MPI:

halo exchange
MPI_Allreduce

and combines them with the local loop parallelism introduced by OpenMP:

parallel stencil
parallel vector updates
parallel local reductions

The resulting solver has two levels of parallelism while keeping the same finite-difference operator and PCG recurrence.

The next step changes the local execution model again. MPI will continue to manage the distributed grid, while the local stencil, vector updates, and reductions will move from CPU threads to CUDA kernels running on the GPU.