The Trapezoidal Rule and Distributed I/O

Having covered the fundamentals of MPI initialization and point-to-point message passing, we can now apply these techniques to a concrete numerical computation: numerical integration via the trapezoidal rule.

Along the way, we explore how to partition mathematical workloads, distinguish between local and global data structures, and handle standard input/output (I/O) in a distributed-memory cluster.

3.2 The Trapezoidal Rule for Numerical Integration

The trapezoidal rule approximates the definite integral of a function:

\int_a^b f(x) \, dx

Geometrically, this represents the area under the curve $y = f(x)$ between the vertical lines $x = a$ , $x = b$ , and the $x$ -axis.

flowchart LR
  subgraph TrapConcept["Trapezoidal Approximation"]
      direction TB
      A["Interval [a, b]"] --> B["Divide into n equal subintervals of width h = (b - a) / n"]
      B --> C["Approximate curve on [x_i, x_{i+1}] with secant line"]
      C --> D["Area of trapezoid i = (h / 2) * [ f(x_i) + f(x_{i+1}) ]"]
      D --> E["Sum all trapezoid areas across interval"]
  end

If we divide the interval $[a, b]$ into $n$ equal subintervals of width:

h = \frac{b - a}{n}

The endpoints of the subintervals are:

x_0 = a, \quad x_1 = a + h, \quad x_2 = a + 2h, \quad \dots, \quad x_n = b

Summing the areas of the $n$ trapezoids yields the composite trapezoidal formula:

\text{Area} \approx h \left[ \frac{f(x_0)}{2} + f(x_1) + f(x_2) + \dots + f(x_{n-1}) + \frac{f(x_n)}{2} \right]

Serial Algorithm

In pseudocode, the sequential computation proceeds as follows:

/* Serial Trapezoidal Rule */
h = (b - a) / n;
approx = (f(a) + f(b)) / 2.0;
for (i = 1; i <= n - 1; i++) {
    x_i = a + i * h;
    approx += f(x_i);
}
approx = h * approx;

3.3 Parallelizing the Trapezoidal Rule

To convert this serial computation into a parallel program, we apply Foster’s four-stage design methodology:

Partitioning: Decompose the overall computation into the finest-grained tasks possible. Here, computing the area of each individual trapezoid represents a basic task, along with an accumulator task that sums the partial areas.
Communication: Establish data flow channels. Each trapezoid evaluation task transmits its calculated area to the summation task.
Aggregation: Group fine-grained tasks into composite tasks. Because the number of trapezoids $n$ is typically orders of magnitude larger than the number of available cores ( $p$ ), calculating a single trapezoid per process would introduce prohibitive communication overhead. We therefore aggregate contiguous groups of $n / p$ trapezoids into blocks.
Mapping: Assign each aggregated block to an MPI process rank.

flowchart TD
  subgraph Foster["Figure 3.5: Task Partitioning & Communication"]
      direction TB
      T0["Compute area trap 0"]
      T1["Compute area trap 1"]
      TDOTS["..."]
      TN["Compute area trap n-1"]
      
      SUM["Add areas into total"]
      
      T0 --> SUM
      T1 --> SUM
      TDOTS --> SUM
      TN --> SUM
  end

Subinterval Calculation per Process

Assuming the number of processes (comm_sz) evenly divides $n$ , each process is assigned:

\text{local\_n} = \frac{n}{\text{comm\_sz}}

Each process calculates its private subinterval endpoints based on its unique my_rank:

\text{local\_a} = a + \text{my\_rank} \times \text{local\_n} \times h

\text{local\_b} = \text{local\_a} + \text{local\_n} \times h

Local vs. Global Variables

In distributed-memory systems, variables exist in distinct address spaces:

Local variables: Variables whose contents differ across processes or have significance only to the process hosting them (e.g., local_a, local_b, local_n, local_int).
Global variables: Conceptual variables that represent identical values across all processes (e.g., the global bounds $a$ and $b$ , and the total trapezoid count $n$ ).

3.4 First MPI Implementation of the Trapezoidal Rule

Below is the complete C source code for Program 3.2 and Program 3.3:

/* Program 3.2: First version of MPI trapezoidal rule */
#include <stdio.h>
#include <mpi.h>

/* Prototype for integration function */
double Trap(double left_endpt, double right_endpt, int trap_count, double base_len);

/* Target mathematical function: f(x) = x^2 */
double f(double x) {
    return x * x;
}

int main(void) {
    int my_rank, comm_sz, n = 1024, local_n;
    double a = 0.0, b = 3.0, h, local_a, local_b;
    double local_int, total_int;
    int source;

    MPI_Init(NULL, NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

    h = (b - a) / n;                 /* h is identical for all processes */
    local_n = n / comm_sz;           /* Number of trapezoids per process */

    local_a = a + my_rank * local_n * h;
    local_b = local_a + local_n * h;
    local_int = Trap(local_a, local_b, local_n, h);

    if (my_rank != 0) {
        /* Workers send their local integral to Process 0 */
        MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    } else {
        /* Process 0 initializes total with its own local portion */
        total_int = local_int;
        for (source = 1; source < comm_sz; source++) {
            MPI_Recv(&local_int, 1, MPI_DOUBLE, source, 0,
                     MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            total_int += local_int;
        }
    }

    if (my_rank == 0) {
        printf("With n = %d trapezoids, our estimate\n", n);
        printf("of the integral from %f to %f = %.15e\n", a, b, total_int);
    }

    MPI_Finalize();
    return 0;
} /* main */

/* Program 3.3: Trap function */
double Trap(double left_endpt, double right_endpt, int trap_count, double base_len) {
    double estimate, x;
    int i;

    estimate = (f(left_endpt) + f(right_endpt)) / 2.0;
    for (i = 1; i <= trap_count - 1; i++) {
        x = left_endpt + i * base_len;
        estimate += f(x);
    }
    estimate = estimate * base_len;

    return estimate;
} /* Trap */

3.5 Dealing with Input and Output in MPI

In the initial implementation, values for $a$ , $b$ , and $n$ were hardcoded. In realistic applications, parameters must be supplied interactively or read from disk.

3.5.1 Console Output (`stdout` and `stderr`)

Most MPI implementations allow all processes in MPI_COMM_WORLD to call printf() and write to standard output. However, MPI does not enforce scheduling or serialization for shared output devices.

Consider Program 3.4:

#include <stdio.h>
#include <mpi.h>

int main(void) {
    int my_rank, comm_sz;
    MPI_Init(NULL, NULL);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

    printf("Proc %d of %d > Does anyone have a toothpick?\n", my_rank, comm_sz);

    MPI_Finalize();
    return 0;
}

Running this with 6 processes reveals non-deterministic output:

Run 1:                               Run 2:
Proc 0 of 6 > Does anyone have...     Proc 0 of 6 > Does anyone have...
Proc 1 of 6 > Does anyone have...     Proc 1 of 6 > Does anyone have...
Proc 2 of 6 > Does anyone have...     Proc 4 of 6 > Does anyone have...
Proc 5 of 6 > Does anyone have...     Proc 2 of 6 > Does anyone have...
Proc 3 of 6 > Does anyone have...     Proc 5 of 6 > Does anyone have...
Proc 4 of 6 > Does anyone have...     Proc 3 of 6 > Does anyone have...

Processes compete concurrently for access to stdout. Their lines interleave unpredictably depending on operating system scheduling and network transmission times.

3.5.2 Console Input (`stdin`)

Unlike standard output, standard input (stdin) is restricted:

The majority of MPI implementations connect stdin only to Process 0.
If multiple processes attempted to read from stdin simultaneously, the system would have no way to determine which process should consume which keystrokes or lines.

Therefore, interactive input follows a coordinator-distribution pattern:

Process 0 prompts the user and calls scanf().
Process 0 transmits the input parameters to each worker process.

Below is Program 3.5, implementing Get_input:

void Get_input(int my_rank, int comm_sz, double* a_p, double* b_p, int* n_p) {
    int dest;

    if (my_rank == 0) {
        printf("Enter a, b, and n\n");
        scanf("%lf %lf %d", a_p, b_p, n_p);

        /* Send values sequentially to every other process */
        for (dest = 1; dest < comm_sz; dest++) {
            MPI_Send(a_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
            MPI_Send(b_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
            MPI_Send(n_p, 1, MPI_INT,    dest, 0, MPI_COMM_WORLD);
        }
    } else {
        /* Worker processes receive parameters from Process 0 */
        MPI_Recv(a_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(b_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(n_p, 1, MPI_INT,    0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
} /* Get_input */

While functional, this approach requires $3 \times (p - 1)$ separate point-to-point messages. In the next section, we examine collective communications, which dramatically reduce communication overhead for global operations.