ayodeji-ibrahim: Implemented blocked and parallel matrix multiplication

Co-authored-by: Muhammad Zahid <muhammad.zahid@example.com>

Author: Deji10

README.md

@@ -235,3 +235,58 @@ git push origin student-name
     - Use small test cases to debug your blocked and parallel implementations.
 
 Good luck, and enjoy optimizing your matrix multiplication!
+
+---
+
+## Performance Results
+
+### Environment
+- **Platform**: GitHub Codespaces (Linux x86_64)
+- **Compiler**: g++ with `-O3 -fopenmp`
+- **Threads**: 4 (OpenMP)
+- **Block Size**: 64
+
+### Results Table
+
+| Case | Dimensions (m × n × p) | Naive (s) | Blocked (s) | Parallel (s) | Blocked Speedup | Parallel Speedup |
+|------|------------------------|-----------|-------------|--------------|-----------------|------------------|
+| 0    | 64 × 64 × 64           | 0.000433  | 0.000355    | 0.001822     | 1.22×           | 0.24×            |
+| 1    | 128 × 64 × 128         | 0.000982  | 0.000858    | 0.000854     | 1.14×           | 1.15×            |
+| 2    | 100 × 128 × 56         | 0.000747  | 0.000543    | 0.000621     | 1.38×           | 1.20×            |
+| 3    | 128 × 64 × 128         | 0.001224  | 0.002144    | 0.001733     | 0.57×           | 0.71×            |
+| 4    | 32 × 128 × 32          | 0.000233  | 0.000143    | 0.001070     | 1.63×           | 0.22×            |
+| 5    | 200 × 100 × 256        | 0.007649  | 0.006110    | 0.007922     | 1.25×           | 0.97×            |
+| 6    | 256 × 256 × 256        | 0.025777  | 0.025024    | 0.018829     | 1.03×           | 1.37×            |
+| 7    | 256 × 300 × 256        | 0.031356  | 0.028475    | 0.031300     | 1.10×           | 1.00×            |
+| 8    | 64 × 128 × 64          | 0.000560  | 0.000402    | 0.000884     | 1.39×           | 0.63×            |
+| 9    | 256 × 256 × 257        | 0.020814  | 0.028982    | 0.022235     | 0.72×           | 0.94×            |
+
+All implementations validated against the reference output (`output.raw`) with tolerance `1e-2`. All cases marked **OK**.
+
+### Analysis
+
+**Correctness**: Every implementation produces identical results to the reference output for all 10 test cases.
+
+**Cache Optimization (Blocked)**:
+- Blocked multiplication gives a **modest improvement (1.0× – 1.6×)** in most cases (0, 1, 2, 4, 5, 6, 7, 8).
+- For tiny matrices (case 4: 32×128×32), blocking gives the largest gain (1.63×) because the entire problem fits inside the cache and blocking reduces redundant memory traffic.
+- Cases 3 and 9 show a *slowdown* with blocking. These results reflect measurement noise in the Codespaces environment — when execution time is in the millisecond range, normal variance can flip the ordering between runs.
+
+**Parallel (OpenMP)**:
+- Parallelization helps **only when the matrices are large enough** to amortize OpenMP thread creation overhead.
+- For tiny matrices (cases 0, 4, 8), parallel is dramatically *slower* (0.22× – 0.63×) because thread setup costs exceed the actual compute work.
+- For larger matrices (case 6: 256×256×256), parallel achieves the expected speedup (~1.37×) using 4 threads.
+- A more impressive parallel speedup would be visible at matrix dimensions above ~512×512, where the problem size exceeds OpenMP's fixed overhead.
+
+**Block Size Choice**: 64 was chosen because modern CPU cache lines are 64 bytes (8 doubles per line). This block size keeps the working set within L1 cache (typically 32–64 KB) without excessive loop overhead.
+
+**Thread Count**: 4 threads were used to match the Codespaces default core allocation.
+
+### Challenges
+
+1. **Small Test Cases**: The provided test cases are too small to fully showcase OpenMP parallelism. On these sizes, the OpenMP setup cost is comparable to (or exceeds) the actual work. Larger matrices (1024×1024 and above) would demonstrate parallel speedups closer to the theoretical 4× on 4 threads.
+
+2. **Measurement Noise**: With timings in the sub-millisecond range, run-to-run variance is significant. Reported numbers are from a single run; averaging across multiple runs would give more stable speedup figures.
+
+3. **Text-format I/O**: The `.raw` files are space-separated text, not binary doubles. Reading was done using `ifstream >> double` with the first two integers as dimensions.
+

main.cpp

@@ -1,114 +1,138 @@
 #include <iostream>
 #include <fstream>
+#include <cstdlib>
+#include <cmath>
 #include <string>
+#include <algorithm>
 #include <omp.h>
-#include <cmath>
 
-void naive_matmul(float *C, float *A, float *B, uint32_t m, uint32_t n, uint32_t p) {
-    //TODO : Implement naive matrix multiplication
-}
+using namespace std;
 
-void blocked_matmul(float *C, float *A, float *B, uint32_t m, uint32_t n, uint32_t p, uint32_t block_size) {
-    // TODO: Implement blocked matrix multiplication
-    // A is m x n, B is n x p, C is m x p
-    // Use block_size to divide matrices into submatrices
+void naive_matmul(double *A, double *B, double *C, int m, int n, int p) {
+    for (int i = 0; i < m; i++) {
+        for (int j = 0; j < p; j++) {
+            C[i*p + j] = 0.0;
+            for (int k = 0; k < n; k++) {
+                C[i*p + j] += A[i*n + k] * B[k*p + j];
+            }
+        }
+    }
 }
 
-void parallel_matmul(float *C, float *A, float *B, uint32_t m, uint32_t n, uint32_t p) {
-    // TODO: Implement parallel matrix multiplication using OpenMP
-    // A is m x n, B is n x p, C is m x p
+void blocked_matmul(double *A, double *B, double *C, int m, int n, int p, int bs) {
+    for (int i = 0; i < m*p; i++) C[i] = 0.0;
+    for (int ii = 0; ii < m; ii += bs) {
+        for (int jj = 0; jj < p; jj += bs) {
+            for (int kk = 0; kk < n; kk += bs) {
+                int ie = min(ii+bs, m);
+                int je = min(jj+bs, p);
+                int ke = min(kk+bs, n);
+                for (int i = ii; i < ie; i++) {
+                    for (int j = jj; j < je; j++) {
+                        for (int k = kk; k < ke; k++) {
+                            C[i*p + j] += A[i*n + k] * B[k*p + j];
+                        }
+                    }
+                }
+            }
+        }
+    }
 }
 
-bool validate_result(const std::string &result_file, const std::string &reference_file) {
-   //TODO : Implement result validation
+void parallel_matmul(double *A, double *B, double *C, int m, int n, int p) {
+    #pragma omp parallel for collapse(2)
+    for (int i = 0; i < m; i++) {
+        for (int j = 0; j < p; j++) {
+            C[i*p + j] = 0.0;
+            for (int k = 0; k < n; k++) {
+                C[i*p + j] += A[i*n + k] * B[k*p + j];
+            }
+        }
+    }
 }
 
-int main(int argc, char *argv[]) {
-    if (argc != 2) {
-        std::cerr << "Usage: " << argv[0] << " <case_number>" << std::endl;
-        return 1;
+double* read_matrix(const string& filename, int& rows, int& cols) {
+    ifstream file(filename);
+    if (!file) {
+        cerr << "Cannot open file" << endl;
+        exit(1);
+    }
+    file >> rows >> cols;
+    double* M = (double*)malloc(rows * cols * sizeof(double));
+    for (int i = 0; i < rows * cols; i++) {
+        file >> M[i];
     }
+    return M;
+}
 
-    int case_number = std::atoi(argv[1]);
-    if (case_number < 0 || case_number > 9) {
-        std::cerr << "Case number must be between 0 and 9" << std::endl;
+int main(int argc, char *argv[]) {
+    if (argc < 2) {
+        cerr << "Usage: program <case 0-9>" << endl;
         return 1;
     }
+    int cn = atoi(argv[1]);
 
-    // Construct file paths
-    std::string folder = "data/" + std::to_string(case_number) + "/";
-    std::string input0_file = folder + "input0.raw";
-    std::string input1_file = folder + "input1.raw";
-    std::string result_file = folder + "result.raw";
-    std::string reference_file = folder + "output.raw";
+    string pA = "data/" + to_string(cn) + "/input0.raw";
+    string pB = "data/" + to_string(cn) + "/input1.raw";
+    string pC = "data/" + to_string(cn) + "/output.raw";
 
-    // TODO Read input0.raw (matrix A)
+    int m, nA, nB, p, mo, po;
+    double *A = read_matrix(pA, m, nA);
+    double *B = read_matrix(pB, nB, p);
+    double *Cref = read_matrix(pC, mo, po);
+    int n = nA;
 
-
-    // TODO Read input1.raw (matrix B)
-
-
-    // Allocate memory for result matrices
-    float *C_naive = new float[m * p];
-    float *C_blocked = new float[m * p];
-    float *C_parallel = new float[m * p];
-
-    // Measure performance of naive_matmul
-    double start_time = omp_get_wtime();
-    naive_matmul(C_naive, A, B, m, n, p);
-    double naive_time = omp_get_wtime() - start_time;
-
-    // TODO Write naive result to file
-
-
-    // Validate naive result
-    bool naive_correct = validate_result(result_file, reference_file);
-    if (!naive_correct) {
-        std::cerr << "Naive result validation failed for case " << case_number << std::endl;
+    if (nA != nB) {
+        cerr << "Dimension mismatch" << endl;
+        return 1;
     }
 
-    // Measure performance of blocked_matmul (use block_size = 32 as default)
-    start_time = omp_get_wtime();
-    blocked_matmul(C_blocked, A, B, m, n, p, 32);
-    double blocked_time = omp_get_wtime() - start_time;
+    cout << "Case " << cn << ": A(" << m << "x" << n << ") * B(" << n << "x" << p << ")" << endl;
 
-    // TODO Write blocked result to file
+    double *Cn = (double*)malloc(m * p * sizeof(double));
+    double *Cb = (double*)malloc(m * p * sizeof(double));
+    double *Cp = (double*)malloc(m * p * sizeof(double));
 
+    omp_set_num_threads(4);
 
-    // Validate blocked result
-    bool blocked_correct = validate_result(result_file, reference_file);
-    if (!blocked_correct) {
-        std::cerr << "Blocked result validation failed for case " << case_number << std::endl;
-    }
+    double t0 = omp_get_wtime();
+    naive_matmul(A, B, Cn, m, n, p);
+    double tn = omp_get_wtime() - t0;
 
-    // Measure performance of parallel_matmul
-    start_time = omp_get_wtime();
-    parallel_matmul(C_parallel, A, B, m, n, p);
-    double parallel_time = omp_get_wtime() - start_time;
+    t0 = omp_get_wtime();
+    blocked_matmul(A, B, Cb, m, n, p, 64);
+    double tb = omp_get_wtime() - t0;
 
-    // TODO Write parallel result to file
+    t0 = omp_get_wtime();
+    parallel_matmul(A, B, Cp, m, n, p);
+    double tp = omp_get_wtime() - t0;
 
-
-    // Validate parallel result
-    bool parallel_correct = validate_result(result_file, reference_file);
-    if (!parallel_correct) {
-        std::cerr << "Parallel result validation failed for case " << case_number << std::endl;
+    bool nok = true, bok = true, pok = true;
+    double eps = 1e-2;
+    for (int i = 0; i < m * p; i++) {
+        if (fabs(Cn[i] - Cref[i]) > eps) nok = false;
+        if (fabs(Cb[i] - Cref[i]) > eps) bok = false;
+        if (fabs(Cp[i] - Cref[i]) > eps) pok = false;
     }
 
-    // Print performance results
-    std::cout << "Case " << case_number << " (" << m << "x" << n << "x" << p << "):\n";
-    std::cout << "Naive time: " << naive_time << " seconds\n";
-    std::cout << "Blocked time: " << blocked_time << " seconds\n";
-    std::cout << "Parallel time: " << parallel_time << " seconds\n";
-    std::cout << "Blocked speedup: " << (naive_time / blocked_time) << "x\n";
-    std::cout << "Parallel speedup: " << (naive_time / parallel_time) << "x\n";
-
-    // Clean up
-    delete[] A;
-    delete[] B;
-    delete[] C_naive;
-    delete[] C_blocked;
-    delete[] C_parallel;
+    cout << "  Naive:    " << tn << " s (" << (nok ? "OK" : "FAIL") << ")" << endl;
+    cout << "  Blocked:  " << tb << " s (" << (bok ? "OK" : "FAIL") << ") speedup: " << tn/tb << "x" << endl;
+    cout << "  Parallel: " << tp << " s (" << (pok ? "OK" : "FAIL") << ") speedup: " << tn/tp << "x" << endl;
+
+    string rp = "data/" + to_string(cn) + "/result.raw";
+    ofstream out(rp);
+    out << m << " " << p << endl;
+    for (int i = 0; i < m * p; i++) {
+        out << Cp[i];
+        if ((i + 1) % p == 0) out << endl;
+        else out << " ";
+    }
 
+    free(A);
+    free(B);
+    free(Cref);
+    free(Cn);
+    free(Cb);
+    free(Cp);
     return 0;
 }