Step 3 — Parallelize the Solver with OpenMP

The sequential solver gives us a numerical reference and, just as importantly, a clear execution structure. We can now look at the expensive loops inside the PCG method and ask which iterations can execute concurrently.

OpenMP is the first parallel model introduced in the assignment. All threads still operate on the same process and share the same arrays in memory:

The grid is therefore stored exactly as it was in the sequential program. There is no domain decomposition, no halo region, and no explicit data exchange between workers.

The main change is how iterations of selected loops are scheduled across CPU threads.

For this solver, the computational work falls into three useful categories:

Understanding these three patterns is enough to parallelize most of the sequential solver.

1. Start with independent grid-point operations

Several loops perform one operation at every grid point without modifying neighboring output values.

For example, the sequential initialization of the right-hand side is

for (int i = 0; i <= M; ++i) {
    for (int j = 0; j <= N; ++j) {
        double x = x_min + i * hx;
        double y = y_min + j * hy;

        f[i][j] =
            inside_region(x, y)
            ? source_func(x, y)
            : 0.0;
    }
}

Each iteration writes to a different element of f. Evaluating f[i][j] does not depend on the value being written by another iteration.

This makes the loop directly suitable for OpenMP:

#pragma omp parallel for schedule(static) collapse(2)
for (int i = 0; i <= M; ++i) {
    for (int j = 0; j <= N; ++j) {
        double x = x_min + i * hx;
        double y = y_min + j * hy;

        f[i][j] =
            inside_region(x, y)
            ? source_func(x, y)
            : 0.0;
    }
}

The same transformation can be applied when initializing the fictitious-domain coefficient:

#pragma omp parallel for schedule(static) collapse(2)
for (int i = 0; i <= M; ++i) {
    for (int j = 0; j <= N; ++j) {
        double x = x_min + i * hx;
        double y = y_min + j * hy;

        k[i][j] =
            inside_region(x, y)
            ? 1.0
            : 1.0 / eps;
    }
}

The reference solver uses this pattern in both initialize_f() and initialize_k(), with schedule(static) and collapse(2). The same OpenMP structure appears throughout the local grid operations in the solver.

2. Why use collapse(2)?

Our grid traversal consists of two nested loops:

for (int i = ...)
    for (int j = ...)
        ...

Without collapse(2), OpenMP distributes iterations of the outer i loop. With

collapse(2)

the two-dimensional iteration space is treated as one larger set of independent iterations.

Conceptually, an \(M\times N\) traversal becomes approximately

individual pieces of work that OpenMP can distribute among threads.

This is useful for a regular rectangular grid because each \((i,j)\) pair performs roughly the same amount of computation.

The reference implementation also uses

schedule(static)

which assigns portions of this iteration space to threads ahead of time.

For these regular numerical loops, the workload at each grid point is very similar, so static scheduling provides a simple distribution with little runtime scheduling overhead.

3. Parallelize the PCG vector updates

The PCG algorithm repeatedly updates complete grid functions.

For the solution,

The sequential implementation is

for (int i = 1; i < M; ++i)
    for (int j = 1; j < N; ++j)
        u[i][j] += alpha * p[i][j];

Every iteration writes only to its own u[i][j], so the loop becomes

#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 residual update follows the same pattern,

so we can write

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

The search-direction update,

is also pointwise:

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

The diagonal preconditioner has the same structure:

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

The current SM25 implementation parallelizes update_u, update_r, update_p, and apply_preconditioner in exactly this way.

These are the easiest operations to parallelize because no iteration needs to modify data produced by another iteration.

4. Parallelize the finite-difference stencil

apply_A() requires a little more thought because each output point reads several neighboring input points:

        v[i][j+1]
             |
v[i-1][j] -- v[i][j] -- v[i+1][j]
             |
        v[i][j-1]

Reading neighboring values does not create a race condition as long as the input array v remains unchanged during the operation.

Each iteration reads from v and writes to one unique location:

out[i][j]

so different grid points can still be evaluated concurrently.

The OpenMP version is therefore:

#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
        );
    }
}

This is one of the central OpenMP kernels in the solver. The numerical stencil has not changed. OpenMP only changes which CPU thread evaluates each grid point.

The reference implementation uses this parallel structure directly in apply_A().

5. Reductions need special treatment

Inner products are different from the previous loops.

Consider

The sequential implementation uses one accumulator:

double rz = 0.0;

for (int i = 1; i < M; ++i)
    for (int j = 1; j < N; ++j)
        rz += r[i][j] * z[i][j];

Simply adding parallel for here would allow multiple threads to update rz at the same time.

That introduces a race condition.

OpenMP provides the reduction clause for this pattern:

double rz = 0.0;

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

Conceptually, each thread accumulates a private partial result,

and OpenMP combines those partial values at the end:

The same pattern is required for

and for the residual norm

For example:

double sum = 0.0;

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

double r_norm = std::sqrt(sum);

The reference implementation uses OpenMP reductions in compute_rz(), compute_p_Ap(), and compute_l2_norm().

This distinction will become even more important in the MPI version. With OpenMP, all threads still share one address space, so OpenMP can combine the thread-local values inside one process. Once the grid is distributed across multiple MPI processes, these reductions will also need communication between processes.

6. What stays sequential?

Adding OpenMP does not mean every line in solve() should execute simultaneously.

The PCG iteration still has an algorithmic order:

apply A to p
      ↓
compute (p, Ap)
      ↓
compute alpha
      ↓
update u and r
      ↓
compute residual norm
      ↓
apply preconditioner
      ↓
compute (r, z)
      ↓
compute beta
      ↓
update p

For example,

cannot be calculated until both scalar products needed for that iteration are available.

Likewise, the next search direction cannot be constructed until the new residual has been computed and preconditioned.

OpenMP therefore parallelizes the work inside these computational stages. The overall PCG dependency graph remains unchanged.

This gives the OpenMP implementation a useful structure:

PCG iteration
│
├── parallel stencil
│
├── parallel reduction
│
├── scalar calculation
│
├── parallel vector update
│
├── parallel reduction
│
├── parallel preconditioner
│
├── parallel reduction
│
├── scalar calculation
│
└── parallel vector update

7. Build the OpenMP executable

The sequential and OpenMP versions can use the same source files.

The sequential build from Step 2 was

g++ -O3 -o task_seq     task.cpp     src/conjugate_gradient.cpp     -Iinclude     -lm     -std=c++11

To activate the OpenMP directives, add

-fopenmp

when compiling:

g++ -O3 -o task_omp     task.cpp     src/conjugate_gradient.cpp     -Iinclude     -lm     -std=c++11     -fopenmp

The number of threads can then be controlled through

export OMP_NUM_THREADS=4

followed by a normal program invocation:

./task_omp 40 40

For example,

export OMP_NUM_THREADS=1
./task_omp 40 40

export OMP_NUM_THREADS=4
./task_omp 40 40

export OMP_NUM_THREADS=16
./task_omp 40 40

The assignment explicitly asks us to test the \(40\times40\) problem with 1, 4, and 16 OpenMP threads and compare the results with the sequential program.

8. Verify numerical correctness first

The first OpenMP experiment should answer a numerical question:

Does changing the number of threads change the problem we are solving?

Run the same \(40\times40\) grid with

sequential
OpenMP × 1 thread
OpenMP × 4 threads
OpenMP × 16 threads

and compare the resulting solution fields.

The main shape of the numerical solution and the convergence behavior should remain consistent across all configurations.

Small differences in scalar reductions can occur because floating-point addition is order dependent. A parallel reduction may combine partial sums in a different order from the sequential loop.

For that reason, validation should use an appropriate numerical tolerance rather than requiring every floating-point value in the CSV files to have an identical bit pattern.

The OpenMP implementation should also converge under the same numerical settings used by the sequential reference.

9. Measure where the time goes

The executable already separates runtime into several categories:

Initialization time
Laplace operator time
Update operator time
Reduction time
Finalize time
Total runtime

This breakdown is useful because OpenMP does not affect every category in the same way.

The Laplace/operator time contains the stencil work from apply_A(), and also from initialize_M_inv(), which reads the same 5-point stencil of k.

The update time contains pointwise operations such as vector updates, coefficient initialization, and preconditioning.

The reduction time contains scalar products and the residual norm.

Initialization and finalization also include work outside the main iterative kernels, such as constructing solver state and saving the final CSV.

Looking only at total runtime can hide these differences. Comparing the timing categories helps reveal which parts benefit from additional CPU threads and which parts begin to limit overall speedup.

10. Compute speedup and efficiency

With the sequential runtime

and an OpenMP runtime using \(p\) threads,

the speedup is

Parallel efficiency can then be written as

For example, the experiment table can be organized as

The assignment first uses the \(40\times40\) case to verify the OpenMP implementation. Its reporting material also includes substantially larger grids, where parallel performance is easier to observe because each thread has more numerical work to perform.

The repository follows the same idea. After the \(40\times40\) OpenMP test, run.sh includes larger experiments on \(400\times600\) and \(800\times1200\) grids.

Small grids are especially useful for correctness checks, while larger grids give a more meaningful view of performance because OpenMP thread-management and synchronization costs become smaller relative to the amount of useful computation.

11. What did OpenMP change?

It is useful to compare the program before moving on to MPI.

The sequential version stores one complete grid and executes every loop on one CPU thread.

The OpenMP version still stores one complete grid in one process. Multiple threads now cooperate on the local loops:

                  shared arrays
        u, f, k, r, z, p, Ap, M_inv
                       │
          ┌────────────┼────────────┐
          │            │            │
       thread 0     thread 1      ...
          │            │
          └──── parallel loops ─────┘

No grid values need to be explicitly sent between threads.

This shared-memory model makes OpenMP a relatively small change to the sequential implementation. Most numerical functions keep the same data structures, formulas, and control flow; the expensive grid loops gain OpenMP directives and reductions gain explicit reduction clauses.

The next stage changes the structure more substantially. With MPI, one process will no longer own the entire computational grid. Each process will own a subdomain, neighboring stencil values will sometimes live in another process, and the scalar reductions used by PCG will need to combine values across the distributed program.

That is the problem we solve in Step 4.