1/*
2 * Double-precision n×n matrix multiply, compiled to WASM with emscripten.
3 * Row-major C = A * B. Kernels:
4 * - matmul_naive: triple loop (sanity baseline)
5 * - matmul_blocked: blocked + SIMD (the "with optimizations" method)
6 * - matmul_blocked_mt: the blocked+SIMD kernel parallelized over rows with
7 * pthreads (only in the -DMATMUL_MT / -pthread build)
8 */
16EMSCRIPTEN_KEEPALIVE
17void matmul_naive(const double* a, const double* b, double* c, int n) {
18 for (int i = 0; i < n; i++) {
19 for (int j = 0; j < n; j++) {
20 double sum = 0.0;
21 for (int k = 0; k < n; k++) sum += a[i * n + k] * b[k * n + j];
22 c[i * n + j] = sum;
23 }
24 }
25}
27/* Block size tuned for L1/L2 cache residency of three double blocks. */
30/* Blocked+SIMD kernel over rows [i0, i1) of C. Each row range is independent
31 * (disjoint output rows), so threads can own disjoint ranges without locking. */
32static void blocked_range(const double* a, const double* b, double* c, int n,
33 int i0, int i1) {
34 memset(c + (size_t)i0 * n, 0, (size_t)(i1 - i0) * n * sizeof(double));
36 for (int jj = 0; jj < n; jj += BS) {
37 int jmax = jj + BS < n ? jj + BS : n;
38 for (int kk = 0; kk < n; kk += BS) {
39 int kmax = kk + BS < n ? kk + BS : n;
40 for (int i = i0; i < i1; i++) {
41 const double* arow = a + (size_t)i * n;
42 double* crow = c + (size_t)i * n;
43 for (int k = kk; k < kmax; k++) {
44 double aik = arow[k];
45 const double* brow = b + (size_t)k * n;
47 v128_t vaik = wasm_f64x2_splat(aik);
48 int j = jj;
49 for (; j + 1 < jmax; j += 2) {
50 v128_t vb = wasm_v128_load(brow + j);
51 v128_t vc = wasm_v128_load(crow + j);
52 vc = wasm_f64x2_add(vc, wasm_f64x2_mul(vaik, vb));
53 wasm_v128_store(crow + j, vc);
54 }
55 for (; j < jmax; j++) crow[j] += aik * brow[j];
57 for (int j = jj; j < jmax; j++) crow[j] += aik * brow[j];
59 }
60 }
61 }
62 }
63}
65EMSCRIPTEN_KEEPALIVE
66void matmul_blocked(const double* a, const double* b, double* c, int n) {
67 blocked_range(a, b, c, n, 0, n);
68}
73typedef struct {
74 const double* a;
75 const double* b;
76 double* c;
77 int n, i0, i1;
78} mm_task;
80static void* mm_worker(void* arg) {
81 mm_task* t = (mm_task*) arg;
82 blocked_range(t->a, t->b, t->c, t->n, t->i0, t->i1);
83 return NULL;
84}
86/* Same blocked+SIMD kernel, but the rows of C are split across nthreads
87 * pthreads (reused from emscripten's pthread pool). */
88EMSCRIPTEN_KEEPALIVE
89void matmul_blocked_mt(const double* a, const double* b, double* c, int n,
90 int nthreads) {
91 if (nthreads < 1) nthreads = 1;
92 if (nthreads > n) nthreads = n;
94 pthread_t th[nthreads];
95 mm_task tasks[nthreads];
96 int rows = (n + nthreads - 1) / nthreads;
97 int nt = 0;
99 for (int i = 0; i < nthreads; i++) {
100 int i0 = i * rows;
101 if (i0 >= n) break;
102 int i1 = i0 + rows;
103 if (i1 > n) i1 = n;
104 tasks[nt] = (mm_task){ a, b, c, n, i0, i1 };
105 if (pthread_create(&th[nt], NULL, mm_worker, &tasks[nt]) == 0) {
106 nt++;
107 } else {
108 /* Pool exhausted — run this chunk inline rather than dropping it. */
109 blocked_range(a, b, c, n, i0, i1);
110 }
111 }
113 for (int i = 0; i < nt; i++) pthread_join(th[i], NULL);
114}