Step 4 — Decompose the Domain with MPI

OpenMP kept the complete grid inside one process and divided loop iterations among threads. MPI changes the ownership of the data itself.

The global \(M\times N\) grid is now divided into rectangular subdomains, and each MPI process stores and updates only one part of the solution. The PCG algorithm remains the same, but two operations now require communication:

1. the finite-difference stencil needs values from neighboring subdomains;
2. scalar products such as \((r,z)\) and \((p,Ap)\) need contributions from every process.

These two requirements lead to the main communication patterns of the MPI solver:

neighboring stencil values
        ↓
    halo exchange

distributed scalar products
        ↓
    MPI_Allreduce

The goal of this step is to preserve the numerical method from the sequential solver while changing the data layout from one global grid to many local grids.

1. From one global grid to local subdomains

Suppose the global computational grid is distributed across four MPI processes. A two-dimensional decomposition may look like this:

global grid

+-------------------+-------------------+
|                   |                   |
|      rank 0       |      rank 1       |
|                   |                   |
+-------------------+-------------------+
|                   |                   |
|      rank 2       |      rank 3       |
|                   |                   |
+-------------------+-------------------+

Each process now owns only the values inside its rectangle.

For a global field \(u\), instead of allocating the entire array on every process, rank \(q\) stores a local field

covering its own part of the global domain together with the boundary data needed by the stencil.

This changes the meaning of the solver dimensions. In the sequential program,

M, N

describe the complete grid. In the MPI solver, each process constructs its own MPIPoissonSolver using the dimensions and physical coordinates of its local subdomain.

The global grid spacing is still

so every subdomain belongs to the same finite-difference discretization.

2. Choose a two-dimensional process grid

If the program runs with \(P\) MPI processes, we first choose a process-grid shape

For example,

P = 4

1 × 4      2 × 2      4 × 1

are all valid factorizations mathematically, but they produce very different local subdomain shapes.

The assignment requires a two-dimensional rectangular decomposition and asks that the ratio between the local numbers of grid nodes in the two directions remain in the range

It also requires neighboring subdomains to differ by at most one grid interval when an exact division is impossible.

The repository handles this in DomainDecomposer.

For every factorization of \(P\), it estimates the local aspect ratio:

double local_M = static_cast<double>(M) / tx;
double local_N = static_cast<double>(N) / ty;
double ratio = local_M / local_N;

Candidates outside the required range are rejected:

if (ratio < 0.5 || ratio > 2.0) {
    continue;
}

Among the remaining candidates, the implementation selects the decomposition whose local blocks are closest to square:

double score = std::fabs(ratio - 1.0);

Once \(p_x\) and \(p_y\) are known, each global direction is divided independently.

A one-dimensional division can be written conceptually as

where \(b\) is the number of blocks. Some blocks receive \(q\) intervals and the remaining blocks receive \(q+1\). This keeps the difference between block sizes to at most one.

The result is a balanced two-dimensional partition of the original grid.

3. Build an MPI Cartesian topology

After choosing

we create a two-dimensional Cartesian communicator:

int dims[2] = {px, py};
int periods[2] = {0, 0};

MPI_Comm cart_comm;

MPI_Cart_create(
    MPI_COMM_WORLD,
    2,
    dims,
    periods,
    1,
    &cart_comm
);

periods = {0, 0} means the computational domain is not periodic. Processes on the outer edge therefore have no neighbor beyond the physical boundary.

Each process has Cartesian coordinates

which identify its position in the process grid.

Conceptually:

coords

(0,0)  (0,1)  (0,2)
(1,0)  (1,1)  (1,2)

These coordinates determine which part of the global grid belongs to each process.

The Cartesian communicator also gives us a direct way to find the four neighbors needed by the five-point stencil:

MPI_Cart_shift(cart_comm, 0, 1, &left,   &right);
MPI_Cart_shift(cart_comm, 1, 1, &bottom, &top);

For a process in the middle of the grid, all four neighbors are valid ranks.

For a process at a physical boundary, MPI returns

MPI_PROC_NULL

for the missing neighbor.

This representation maps naturally to the finite-difference stencil.

4. Add one layer of halo cells

A local process can evaluate most stencil points using values already stored inside its own subdomain.

The difficulty appears near an MPI boundary.

Suppose two ranks meet along the \(x\) direction:

rank A                         rank B

owned cells                    owned cells
...  a  b  c  |               d  e  f  ...
              boundary

To evaluate the stencil at c, rank A needs the value d.

To evaluate the stencil at d, rank B needs the value c.

We therefore extend the local arrays with halo cells (also called ghost cells):

rank A                              rank B

owned             halo        halo             owned
...  a  b  c  |   d'          c'   |  d  e  f  ...

Here d' is a local copy of rank B’s boundary value, while c' is a local copy of rank A’s value.

With this extra layer, the local stencil code can keep the same indexing pattern used by the sequential solver:

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

The MPI-specific work happens before the stencil is evaluated: neighboring processes synchronize the halo values.

In task_mpi.cpp, the owned subdomain bounds are first computed from the Cartesian coordinates. If a neighbor exists in one direction, the local range is extended by one grid position in that direction before constructing MPIPoissonSolver.

This gives the local solver enough storage for the halo layer while preserving the physical boundary of the original problem.

5. Exchange halo values

MPIPoissonSolver stores send and receive buffers for the four directions:

send_left       recv_left
send_right      recv_right
send_bottom     recv_bottom
send_top        recv_top

Before applying a stencil to a distributed field, each process copies its owned boundary values into these buffers.

For example, the left boundary is packed with

for (int j = 0; j <= N; ++j)
    send_left[j] = field[1][j];

The right boundary uses

field[M-1][j]

and the bottom and top boundaries use the corresponding rows.

The communication itself uses nonblocking point-to-point MPI operations. For the left neighbor, for example:

MPI_Request reqs[8];
int rq = 0;

if (left_rank != MPI_PROC_NULL) {
    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++]);
}

The right, bottom, and top neighbors use the same pattern with their own send and receive buffers.

After posting all transfers, the process waits for them to complete:

MPI_Waitall(rq, reqs, MPI_STATUSES_IGNORE);

The received data is then copied into the halo cells:

field[0][j] = recv_left[j];
field[M][j] = recv_right[j];

field[i][0] = recv_bottom[i];
field[i][N] = recv_top[i];

This gives every process the neighboring values required by its local five-point stencil.

The communication pattern can be summarized as

        top
         ↑
         │
left ← local rank → right
         │
         ↓
       bottom

The amount of local computation grows with the area of the subdomain, while halo communication depends on its boundary. Keeping subdomains reasonably compact therefore also gives a useful communication shape.

6. Decide when halo exchange is required

The halo does not need to be exchanged before every operation.

Pointwise updates such as

operate only on local entries.

The finite-difference operator needs neighboring values, so its input field must have valid halos.

During initialization, the MPI solver exchanges the coefficient field before constructing the diagonal preconditioner:

initialize_k();
exchange_halo(k);
initialize_M_inv();

It also exchanges the initial solution before computing the first residual:

exchange_halo(u);
initialize_r();

Inside the PCG loop, the search direction must be synchronized before computing

exchange_halo(p);
apply_A(p, A_p);

This gives a simple rule:

Exchange a distributed field before an operation reads that field across a subdomain boundary.

The local numerical kernels remain very close to the sequential implementation once the halo contains the correct neighboring values.

7. Turn local reductions into global reductions

The stencil introduces nearest-neighbor communication. PCG introduces a second communication pattern through its scalar products.

Consider

After domain decomposition, each process owns only part of this sum:

The complete value is

Each rank therefore begins by computing a local scalar:

double local_rz = 0.0;

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

Then MPI combines the contributions from all processes:

double global_rz = 0.0;

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

MPI_Allreduce is especially convenient here because every process needs the final scalar in order to continue the PCG iteration.

The same pattern is used for

and

So the three distributed reductions become

local (r,z)
      ↓
MPI_Allreduce
      ↓
global (r,z)

local (p,Ap)
      ↓
MPI_Allreduce
      ↓
global (p,Ap)

local ||r||²
      ↓
MPI_Allreduce
      ↓
global ||r||²

After the global scalar is available, every rank computes the same

and

then independently updates its local portion of the PCG vectors.

8. The distributed PCG iteration

With halo exchange and global reductions in place, the PCG iteration becomes

exchange halo of p
        ↓
local apply_A(p)
        ↓
local (p, Ap)
        ↓
MPI_Allreduce
        ↓
global alpha
        ↓
local update u
local update r
        ↓
local ||r||²
        ↓
MPI_Allreduce
        ↓
global residual norm
        ↓
local preconditioner
        ↓
local (r, z)
        ↓
MPI_Allreduce
        ↓
global beta
        ↓
local update p

The code follows the same sequence:

for (iter = 0; r_norm > tolerance && iter < max_iter; ++iter) {
    exchange_halo(p);
    apply_A(p, A_p);

    double p_Ap = compute_p_Ap();
    double alpha = rz_prev / p_Ap;

    update_u(alpha);
    update_r(alpha);

    r_norm = compute_l2_norm();

    apply_preconditioner();

    double rz = compute_rz();
    double beta = rz / rz_prev;

    update_p(beta);
    rz_prev = rz;
}

The numerical PCG recurrence is unchanged. The distributed version adds communication at the points where the recurrence depends on data owned by other processes.

This is the central design of the MPI solver:

local numerical kernels
        +
nearest-neighbor halo exchange
        +
global scalar reductions

9. Reconstruct the global solution

During the solve, each process stores only its local solution. For plotting and comparison with the sequential result, we eventually want one global \((M+1)\times(N+1)\) array again.

The repository performs this reconstruction on rank 0.

Each process first describes the owned part of its local domain:

struct Domain2D {
    int x_min;
    int x_max;
    int y_min;
    int y_max;
    int size;
};

These descriptions are collected with

MPI_Gather(&domain, sizeof(Domain2D), MPI_BYTE,
           domains.data(), sizeof(Domain2D), MPI_BYTE, 0, cart_comm);

so rank 0 knows where every block belongs in the global grid.

Each process then flattens its owned local solution into a one-dimensional buffer. Because subdomains can differ slightly in size, the buffers do not necessarily contain the same number of elements.

The values are therefore collected using

MPI_Gatherv(sendbuf.data(), sendbuf.size(), MPI_DOUBLE,
            recv_buf.data(), recv_counts.data(), displs.data(),
            MPI_DOUBLE, 0, cart_comm);

with one receive count and displacement for every rank.

Rank 0 reconstructs

by placing each received block back at its global coordinates, and finally writes the same CSV format used by the sequential and OpenMP versions.

This makes direct numerical comparison between implementations much easier.

10. Build and run the MPI version

The MPI-only executable is compiled with mpicxx:

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

There is no -fopenmp flag in this build. The MPI source already contains OpenMP directives that will be used in the hybrid version later; without OpenMP enabled, the local loops execute sequentially inside each MPI process.

Run the program with mpirun.

For example:

mpirun -np 1 ./task_mpi 40 40
mpirun -np 2 ./task_mpi 40 40
mpirun -np 4 ./task_mpi 40 40

These are the process counts required for the initial MPI correctness test in the assignment.

For every run, inspect the reported process-grid decomposition and verify that the program converges to the same numerical problem as the sequential reference.

11. Measure communication separately

MPI introduces costs that did not exist in the sequential solver.

The current implementation records separate timing categories for

Laplace operator
Update operations
Local reductions
MPI halo exchange
MPI Allreduce
Finalization
Total runtime

This distinction will be useful later when we analyze scaling.

The total iteration time now contains several qualitatively different pieces:

Increasing the number of processes reduces the amount of local grid work per rank, while communication remains part of every PCG iteration.

For the \(40\times40\) cases in this step, the immediate goal is numerical correctness across 1, 2, and 4 MPI processes. Larger-grid timing and scaling results can be analyzed together with the other implementations in the final evaluation step.

A practical debugging hint for this stage: a halo-exchange bug usually produces correct-looking results at -np 1 (no exchange ever happens) and only fails once the subdomains actually communicate at -np 2 or -np 4. If the single-process run converges but the multi-process runs do not, the halo exchange is the first place to look.

12. What changed when we moved to MPI?

The transition from OpenMP to MPI changes the ownership model of the solver.

With OpenMP:

one process
    ↓
one complete grid
    ↓
many threads

With MPI:

global grid
    ↓
many subdomains
    ↓
one subdomain per process

The finite-difference method and PCG equations stay the same. Their data dependencies now cross process boundaries.

The stencil dependency becomes a halo exchange:

neighbor data
    ↓
MPI_Isend / MPI_Irecv

The scalar dependency becomes a global reduction:

local scalar
    ↓
MPI_Allreduce
    ↓
global scalar

Once these two communication patterns are implemented correctly, the sequential PCG solver becomes a distributed-memory solver.

The next step keeps this MPI decomposition and adds shared-memory parallelism inside every process, producing the hybrid MPI+OpenMP version.