/*
 * 03-matrix-multiply.c  —  LAB 1: Matrix Multiplication with OpenMP
 * Session 10 · 2311CSC501J Parallel Processing · Unit IV — OpenMP in Depth
 *
 * WHAT THIS SHOWS  (this is Lab 1 for the course)
 *   Matrix multiplication C = A x B is the textbook "embarrassingly parallel"
 *   problem. Every output cell C[i][j] is an independent dot product of row i
 *   of A with column j of B — no two output cells depend on each other. That
 *   means we can compute them all at the same time.
 *
 *   We compute the product TWICE:
 *       1. serial     — one thread, the plain triple loop.
 *       2. parallel   — the SAME triple loop with ONE line added:
 *                          #pragma omp parallel for
 *                       on the OUTER (row) loop, so each thread owns a band of
 *                       rows and fills them independently.
 *
 *   We time both with omp_get_wtime() (a high-resolution wall-clock timer that
 *   ships with OpenMP), check that the two results agree, and print the speedup
 *   Speedup = T(serial) / T(parallel).
 *
 * WHY IT IS NOT A PERFECT Nx SPEEDUP  (discuss this in class)
 *   - Memory bandwidth: all threads hammer the same RAM; the bus saturates.
 *   - Cache behaviour: reading B column-by-column is cache-unfriendly.
 *   - Fork/join + scheduling overhead, and any serial parts (allocation, init).
 *   - Turbo/thermal effects: one busy core clocks higher than eight busy cores.
 *   This is Amdahl's Law and coordination cost (Session 08) made concrete.
 *
 * DATA SCOPING — the #1 source of OpenMP bugs
 *   In the parallel loop, i, j, k are loop indices. i is made private
 *   automatically (it is the loop variable of the omp for). j and k MUST be
 *   private too — if two threads shared one k, they would corrupt each other's
 *   inner loop. We declare j and k INSIDE the loop body (C99), which makes them
 *   private by construction. A, B, C, and n are shared (read-only, or each
 *   thread writes disjoint cells of C, so there is no race).
 *
 * HOW TO COMPILE & RUN
 *       gcc -O2 -fopenmp 03-matrix-multiply.c -o matmul
 *       ./matmul
 *       OMP_NUM_THREADS=1 ./matmul     # force serial team -> speedup ~1x
 *       OMP_NUM_THREADS=4 ./matmul     # 4 threads
 *       OMP_NUM_THREADS=8 ./matmul     # 8 threads
 *
 *   -O2 matters here: without optimization the triple loop is very slow.
 *   macOS (Apple clang): clang -O2 -Xpreprocessor -fopenmp -lomp 03-matrix-multiply.c -o matmul
 *   Or paste into OnlineGDB (language C) and press Run — OpenMP is already on.
 *   (On OnlineGDB, keep N at 512 or drop to 256 so it finishes quickly.)
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

#define N 512    /* matrix dimension (N x N). 512 or 1024 are good class sizes. */

/* Allocate one contiguous N*N block of doubles and index it as m[i*N + j].
 * Contiguous memory is cache-friendlier than an array of row pointers. */
static double *alloc_matrix(void) {
    double *m = (double *)malloc((size_t)N * N * sizeof(double));
    if (!m) {
        fprintf(stderr, "Out of memory allocating a %dx%d matrix.\n", N, N);
        exit(EXIT_FAILURE);
    }
    return m;
}

/* Fill A and B with simple, repeatable values (no rand, so results are stable). */
static void init_matrices(double *A, double *B) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            A[i * N + j] = (double)((i + j) % 10) + 1.0;   /* 1..10 */
            B[i * N + j] = (double)((i * j) % 7)  + 1.0;   /* 1..7  */
        }
    }
}

/* SERIAL matrix multiply: the classic triple loop, one thread. */
static void multiply_serial(const double *A, const double *B, double *C) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            double sum = 0.0;                 /* accumulate the dot product */
            for (int k = 0; k < N; k++) {
                sum += A[i * N + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

/* PARALLEL matrix multiply: identical maths, ONE pragma added.
 * schedule(static) because every row is the same amount of work (balanced).
 * default(none) forces us to state the scope of every variable — good hygiene
 * that catches accidental sharing bugs at compile time. */
static void multiply_parallel(const double *A, const double *B, double *C) {
    #pragma omp parallel for schedule(static) default(none) shared(A, B, C)
    for (int i = 0; i < N; i++) {          /* i is private (the omp for index)   */
        for (int j = 0; j < N; j++) {      /* j declared here -> private          */
            double sum = 0.0;              /* sum is local -> private             */
            for (int k = 0; k < N; k++) {  /* k declared here -> private          */
                sum += A[i * N + k] * B[k * N + j];
            }
            C[i * N + j] = sum;            /* each thread writes disjoint rows    */
        }
    }
}

/* Compare two result matrices; return the largest absolute difference.
 * Floating-point addition can reorder, so we allow a tiny tolerance. */
static double max_difference(const double *X, const double *Y) {
    double worst = 0.0;
    for (int i = 0; i < N * N; i++) {
        double d = fabs(X[i] - Y[i]);
        if (d > worst) worst = d;
    }
    return worst;
}

int main(void) {
    double *A  = alloc_matrix();
    double *B  = alloc_matrix();
    double *Cs = alloc_matrix();   /* result from the serial run   */
    double *Cp = alloc_matrix();   /* result from the parallel run */

    init_matrices(A, B);

    /* Report the team size OpenMP will use (set via OMP_NUM_THREADS). */
    int threads = 0;
    #pragma omp parallel
    {
        #pragma omp single
        threads = omp_get_num_threads();
    }

    printf("Matrix multiply  C = A x B   (%d x %d, doubles)\n", N, N);
    printf("OpenMP team size: %d thread(s)\n\n", threads);

    /* --- Time the serial version --- */
    double t0 = omp_get_wtime();
    multiply_serial(A, B, Cs);
    double t1 = omp_get_wtime();
    double serial_time = t1 - t0;
    printf("Serial   : %.4f s\n", serial_time);

    /* --- Time the parallel version --- */
    double t2 = omp_get_wtime();
    multiply_parallel(A, B, Cp);
    double t3 = omp_get_wtime();
    double parallel_time = t3 - t2;
    printf("Parallel : %.4f s\n", parallel_time);

    /* --- Correctness: the two results must match --- */
    double diff = max_difference(Cs, Cp);
    printf("\nMax difference between serial and parallel result: %.2e\n", diff);
    printf("%s\n", (diff < 1e-6) ? "CORRECT — parallel result matches serial."
                                 : "MISMATCH — a scoping/race bug crept in!");

    /* --- Speedup --- */
    if (parallel_time > 0.0) {
        double speedup = serial_time / parallel_time;
        printf("\nSpeedup  = T(serial) / T(parallel) = %.2fx on %d thread(s)\n",
               speedup, threads);
        printf("Efficiency = speedup / threads      = %.0f%%\n",
               100.0 * speedup / threads);
        printf("\nNot a perfect %dx? Memory bandwidth, cache misses, and fork/join\n", threads);
        printf("overhead eat into it — exactly the coordination cost from Session 08.\n");
    }

    free(A); free(B); free(Cs); free(Cp);
    return 0;
}