Write a parallel program for matrix multiplication.

Concept

In sequential matrix multiplication C = A × B (all n×n), each output element C[i][j] = Σₖ A[i][k]·B[k][j] depends only on row i of A and column j of B, never on any other output element. This makes the outer (i, j) index space embarrassingly parallel: the n² independent output elements can be computed by different processors/threads simultaneously, using either the shared-memory model (OpenMP, threads reading a common A and B) or the message-passing model (MPI, each process holding a partition of A and B). The only inherently sequential work is the inner k-loop's accumulation for a single output element — but that accumulation itself can be done independently and concurrently by whichever thread/process owns that element, so it never forces threads to wait on each other.

Application — shared-memory (OpenMP) program

A standard shared-memory solution parallelises the outer row loop and lets every thread compute complete rows of C on its own:

#include <stdio.h>
#include <omp.h>
#define N 500

double A[N][N], B[N][N], C[N][N];

int main() {
    int i, j, k;

    /* ... initialise A and B ... */

    #pragma omp parallel for private(i, j, k) schedule(static)
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            double sum = 0.0;
            for (k = 0; k < N; k++) {
                sum += A[i][k] * B[k][j];
            }
            C[i][j] = sum;
        }
    }

    return 0;
}

1. The compiler directive #pragma omp parallel for forks a team of threads and splits the iterations of the outer i-loop among them — each thread owns a contiguous band of rows.

2. Each thread runs the full (j, k) nest for its own rows. private(i, j, k) gives every thread its own copies of the loop counters i, j, k (which are declared once, outside the parallel region, so without this clause they would be shared and clash); sum needs no explicit privatisation because it is declared inside the loop body itself — each thread executing that block automatically gets its own instance on its own stack. Because no thread ever writes into another thread's row of C, there is no data race and therefore no lock or barrier is needed inside the loop.

3. A and B are only read, never written, during the parallel region, so they can be safely shared by all threads without synchronisation.

4. schedule(static) divides the N row-indices into equal contiguous chunks up front (low overhead, ideal here since every row costs the same N·N multiply-adds).

5. With p threads/processors the compute cost falls from O(n³) sequential time to O(n³/p), plus a small, one-time thread-fork overhead.

Cross-check — worked trace and message-passing (MPI) variant

Trace on a tiny 2×2 case with 2 threads, A = [[1, 2], [3, 4]] and B = [[5, 6], [7, 8]], schedule(static) assigning row i = 0 to thread T0 and row i = 1 to thread T1:

• T0 computes row 0 independently: C[0][0] = 1·5 + 2·7 = 19, and C[0][1] = 1·6 + 2·8 = 22.

• T1 computes row 1 independently and concurrently: C[1][0] = 3·5 + 4·7 = 43, and C[1][1] = 3·6 + 4·8 = 50.

• Neither thread reads or writes the other's row, so the two rows can genuinely run in parallel; the assembled result C = [[19, 22], [43, 50]] is exactly what the single-threaded triple loop would produce, confirming correctness.

On a distributed-memory cluster the same row-partition idea is realised with explicit communication instead of shared arrays: the root process scatters (roughly) N/p rows of A to every process (MPI_Scatter), broadcasts the whole of B to every process (MPI_Bcast), each process computes its own block of rows of C completely independently using the identical triple loop above, and the root gathers the row-blocks back (MPI_Gather). If N is not exactly divisible by p, MPI_Scatterv/MPI_Gatherv (which take explicit per-process counts and displacements) must be used instead so that every row is assigned to exactly one process and none is dropped. The result is identical to the sequential C, since each process still computes exactly the same C[i][j] = Σₖ A[i][k]·B[k][j] sum — parallelism only changes which processor performs which independent sum, never the value being summed, so correctness follows directly from the same argument as the Concept section.