Step 2 — Build the Sequential Reference Solver

We now have the mathematical ingredients of the solver: a rectangular grid, the fictitious-domain coefficient, a finite-difference operator, and the PCG iteration. The next task is to connect them into one program that can take a grid size, compute a solution, and save the result.

This first implementation is deliberately simple. Every operation runs in a well-defined order on a single process, which makes it much easier to inspect the numerical method and establish a reference result for the later parallel versions.

The execution path we want is:

read M and N
    ↓
construct PoissonSolver
    ↓
initialize grid data
    ↓
run PCG
    ↓
check convergence
    ↓
save the solution

1. Start with a small executable

The top-level program does not need to contain the numerical method itself. Its main responsibility is to define the concrete problem, create the solver, start the computation, and write the result.

For this assignment, the physical region and source function can be defined as

bool region(double x, double y) {
    return (x > 0 && y > 0 && y < 3 && (y + 3 * x < 9));
}

double f_func(double x, double y) {
    return 1.0;
}

The executable then reads the grid dimensions from the command line:

int M = atoi(argv[1]);
int N = atoi(argv[2]);

and creates the solver over the rectangular computational domain \([0,3]\times[0,3]\):

PoissonSolver solver(
    M, N,
    0.0, 3.0,
    0.0, 3.0,
    region,
    f_func,
    timer
);

solver.solve();

This separation is useful because task.cpp describes which problem is being solved, while PoissonSolver contains how the numerical problem is solved.

A different geometry can therefore reuse the same solver by changing region(). A different source term can reuse it by changing f_func().

2. Organize the solver state

The PCG algorithm needs several grid functions at the same time. A convenient way to manage them is to keep them inside a solver class.

The reference implementation stores

std::vector<std::vector<double>> u;
std::vector<std::vector<double>> f;
std::vector<std::vector<double>> r;
std::vector<std::vector<double>> z;
std::vector<std::vector<double>> k;
std::vector<std::vector<double>> p;
std::vector<std::vector<double>> A_p;
std::vector<std::vector<double>> M_inv;

Each array has a direct numerical meaning:

All of these arrays have size

because the grid includes both the interior nodes and the outer boundary.

The constructor also computes the grid spacing

and the fictitious-domain parameter

In code:

hx = (x_max - x_min) / M;
hy = (y_max - y_min) / N;

eps = std::max(hx, hy) * std::max(hx, hy);

hx2 = hx * hx;
hy2 = hy * hy;

The arrays can initially be filled with zero, except for the inverse preconditioner, whose entries are computed after k is available (M_inv is allocated with ones and then overwritten by initialize_M_inv()).

3. Initialize the numerical problem

Before the PCG iteration starts, we need to construct the discrete data that represents the PDE.

The right-hand side is initialized by evaluating each grid point:

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

For this problem, points inside the trapezoid receive

while points in the fictitious region receive

The coefficient field is initialized in the same way:

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

After these two steps, the geometric problem has become two rectangular arrays:

f[i][j]    source term
k[i][j]    domain coefficient

The rest of the solver no longer needs to reason explicitly about the shape of the trapezoid.

4. Implement the discrete operator

The most important numerical kernel in the sequential solver is the operation

We implement it as a function instead of constructing the complete matrix \(A\).

For every interior grid point, first estimate the coefficient on the four edges surrounding the node:

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]);

Then apply the variable-coefficient finite-difference stencil:

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

For one point \((i,j)\), this operation reads only

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

so the operator has a local five-point stencil structure.

This function will later become one of the main targets for parallelization, so it is worth keeping it isolated and easy to test.

5. Build the diagonal preconditioner

The PCG method also needs an inexpensive approximation to \(A^{-1}\).

We use the inverse of the diagonal of the discrete operator. At every interior point,

The preconditioner stores

The implementation is therefore:

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]);

        double diag =
            (kx_plus + kx_minus) / hx2
          + (ky_plus + ky_minus) / hy2;

        M_inv[i][j] = 1.0 / diag;
    }
}

Applying it during PCG is then just

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 preconditioner is created once and reused throughout the iteration.

6. Initialize the PCG iteration

The solution array starts from zero,

The initial residual is

We can compute it using the same apply_A() operation:

std::vector<std::vector<double>> A_u(
    M + 1,
    std::vector<double>(N + 1, 0.0)
);

apply_A(u, A_u);

for (int i = 1; i < M; ++i)
    for (int j = 1; j < N; ++j)
        r[i][j] = f[i][j] - A_u[i][j];

Then apply the diagonal preconditioner:

apply_preconditioner();

which gives

The first search direction is

In code:

p = z;

Finally compute the initial scalar product

which will be reused when calculating the first PCG step.

7. Implement the scalar reductions

PCG repeatedly needs scalar products over the complete grid.

For example,

The sequential implementation is simply:

double compute_rz() {
    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];

    return rz;
}

Likewise,

becomes

double compute_p_Ap() {
    double p_Ap = 0.0;

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

    return p_Ap;
}

The residual norm can be computed as

double compute_l2_norm() {
    double sum = 0.0;

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

    return std::sqrt(sum);
}

These loops look simple, but they are structurally different from the stencil and vector-update loops: every grid point contributes to one shared scalar.

That distinction becomes important when we introduce parallel execution.

8. Assemble the PCG loop

We can now put all the operations together.

Initialization proceeds as

initialize_f();
initialize_k();
initialize_M_inv();

initialize_r();
apply_preconditioner();
initialize_p();

double rz_prev = compute_rz();
double r_norm = compute_l2_norm();

Then the iterative part becomes:

for (int iter = 0;
     r_norm > tolerance && iter < max_iter;
     ++iter) {

    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 corresponding vector updates are small:

u[i][j] += alpha * p[i][j];

r[i][j] -= alpha * A_p[i][j];

and

p[i][j] = z[i][j] + beta * p[i][j];

The complete iteration can therefore be summarized as

p
│
├── apply A ───────────────► Ap
│
├── (p, Ap)
│
├── alpha
│
├── update u
│
├── update r
│
├── ||r||
│
├── apply preconditioner ──► z
│
├── (r, z)
│
├── beta
│
└── update p
        │
        └──────── next iteration

At this stage, every mathematical equation introduced in Step 1 has a corresponding operation in the program.

9. Build the sequential executable

The repository uses the same solver source for both the sequential and OpenMP stages.

The sequential executable can be compiled with

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

Notice that this command does not contain

-fopenmp

so the compiler produces the sequential baseline.

The current source file already contains OpenMP directives for later stages. When the same file is compiled without OpenMP support, those directives do not activate parallel execution. This allows the assignment to use the same numerical implementation for the sequential and shared-memory versions.

For a cleaner teaching implementation, you can also begin with the loops without any OpenMP directives and add them explicitly in Step 3.

10. Run the required baseline cases

The first assignment asks us to run the sequential solver on

With the executable above:

./task_seq 10 10
./task_seq 20 20
./task_seq 40 40

The repository automates exactly this sequence in the first stage of run.sh.

Each run should report whether PCG converged and should produce a solution file such as

solution/solution_M_10_N_10.csv
solution/solution_M_20_N_20.csv
solution/solution_M_40_N_40.csv

The top-level program saves the complete \((M+1)\times(N+1)\) solution grid as CSV after solver.solve() completes.

11. Validate the sequential baseline

Before using this solver as the foundation for later versions, check three things.

First, convergence. The residual norm should decrease until the stopping criterion is reached. The reference implementation periodically writes the residual to debug.log, which provides a simple way to inspect the iteration.

Second, the numerical field. Load the CSV result and inspect the solution as a 2D heat map or a 3D surface. Refining the grid should preserve the main shape of the solution while giving a finer numerical representation.

Third, repeatability. Running the same sequential executable with the same grid should give the same result up to the expected floating-point behavior.

Keep the execution time as well. It will later provide the baseline

used to calculate parallel speedup,

Once these checks are complete, we have a reference implementation with a known numerical result and a known execution path.

The next step is to examine these loops one by one and decide which can be executed concurrently. That will lead directly to the OpenMP version.