/* * dgemm.c - mini DGEMM for the PowerPC 604e: the LINPACK-heritage workload. * * C = A * B for square N x N row-major doubles, three ways: * * ikj inner loop over j (contiguous rows; xlc contracts c[j]+=a*b[j] * into FMAs) - the cache-friendly textbook form * dot4 classic dot-product form, k unrolled 4 ways with 4 partial * sums - FMA-friendly but strides B by N (memory-system stress) * blocked ikj inside 32x32 blocks (32*32*8*3 = 24 KB, fits the 32 KB L1) * * N straddles the cliffs DeepSeek's bench mapped: 32 KB L1 (N~64 for three * matrices), 1 MB L2 (N~200), DRAM beyond. FLOPs = 2*N^3. * * Correctness: blocked and dot4 checksums are compared against ikj's; FMA * contraction may move the last ulps, so the check is relative, 1e-9. * * Portable C89. Build on obsidian: * /usr/vac/bin/xlc -O2 -qarch=604 -qtune=604 -o dgemm dgemm.c -lm * * kimi-k3, October 2026 */ #include #include #include #include #include #define BS 32 /* dot4 is capped at 512: above that the stride-8 KB column walk degenerates * into TLB-thrashing DRAM latency (over ten minutes per multiply at N=1024 * on obsidian - measured the hard way). That collapse is itself the finding * the variant exists to show; blocked is the answer to it. */ #define DOT4_MAX_N 512 static unsigned long lcg = 0x12345678UL; static double frand(void) /* uniform in [-1, 1) */ { lcg = lcg * 1103515245UL + 12345UL; return ((double)((lcg >> 16) & 0x7fff) / 16384.0) - 1.0; } static void fill(double *m, long n) { long i; for (i = 0; i < n * n; i++) m[i] = frand(); } /* ------------------------------------------------------------- variants */ static void mul_ikj(const double *A, const double *B, double *C, long n) { long i, k, j; for (i = 0; i < n; i++) { double *c = C + i * n; for (k = 0; k < n; k++) { const double a = A[i * n + k]; const double *b = B + k * n; for (j = 0; j < n; j++) c[j] += a * b[j]; } } } static void mul_dot4(const double *A, const double *B, double *C, long n) { long i, j, k; for (i = 0; i < n; i++) { const double *a = A + i * n; for (j = 0; j < n; j++) { const double *b = B + j; double s0 = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0; for (k = 0; k + 3 < n; k += 4) { s0 += a[k] * b[k * n]; s1 += a[k + 1] * b[(k + 1) * n]; s2 += a[k + 2] * b[(k + 2) * n]; s3 += a[k + 3] * b[(k + 3) * n]; } for (; k < n; k++) s0 += a[k] * b[k * n]; C[i * n + j] = s0 + s1 + s2 + s3; } } } static void mul_blocked(const double *A, const double *B, double *C, long n) { long ii, kk, jj, i, k, j; for (ii = 0; ii < n; ii += BS) for (kk = 0; kk < n; kk += BS) for (jj = 0; jj < n; jj += BS) { long imax = ii + BS < n ? ii + BS : n; long kmax = kk + BS < n ? kk + BS : n; long jmax = jj + BS < n ? jj + BS : n; for (i = ii; i < imax; i++) { double *c = C + i * n; for (k = kk; k < kmax; k++) { const double a = A[i * n + k]; const double *b = B + k * n; for (j = jj; j < jmax; j++) c[j] += a * b[j]; } } } } /* 4x4 register-blocked rank-1 updates: C(i..i+3, j..j+3) lives in 16 * accumulators across the whole k loop, so each k costs 8 loads for 16 * FMAs and C never touches memory. ikj needs 3 LSU ops per FMA and starves * the FPU (~1 load/cycle); here the FPU is meant to be the bottleneck. * Requires n % 4 == 0 (every size in the table is). */ static void mul_reg4(const double *A, const double *B, double *C, long n) { long i, j, k; for (i = 0; i < n; i += 4) { const double *ap0 = A + i * n; const double *ap1 = ap0 + n; const double *ap2 = ap1 + n; const double *ap3 = ap2 + n; for (j = 0; j < n; j += 4) { double c00 = 0.0, c01 = 0.0, c02 = 0.0, c03 = 0.0; double c10 = 0.0, c11 = 0.0, c12 = 0.0, c13 = 0.0; double c20 = 0.0, c21 = 0.0, c22 = 0.0, c23 = 0.0; double c30 = 0.0, c31 = 0.0, c32 = 0.0, c33 = 0.0; const double *bp = B + j; for (k = 0; k < n; k++) { double a0 = ap0[k], a1 = ap1[k], a2 = ap2[k], a3 = ap3[k]; double b0 = bp[0], b1 = bp[1], b2 = bp[2], b3 = bp[3]; c00 += a0 * b0; c01 += a0 * b1; c02 += a0 * b2; c03 += a0 * b3; c10 += a1 * b0; c11 += a1 * b1; c12 += a1 * b2; c13 += a1 * b3; c20 += a2 * b0; c21 += a2 * b1; c22 += a2 * b2; c23 += a2 * b3; c30 += a3 * b0; c31 += a3 * b1; c32 += a3 * b2; c33 += a3 * b3; bp += n; } C[i * n + j] = c00; C[i * n + j + 1] = c01; C[i * n + j + 2] = c02; C[i * n + j + 3] = c03; C[(i + 1) * n + j] = c10; C[(i + 1) * n + j + 1] = c11; C[(i + 1) * n + j + 2] = c12; C[(i + 1) * n + j + 3] = c13; C[(i + 2) * n + j] = c20; C[(i + 2) * n + j + 1] = c21; C[(i + 2) * n + j + 2] = c22; C[(i + 2) * n + j + 3] = c23; C[(i + 3) * n + j] = c30; C[(i + 3) * n + j + 1] = c31; C[(i + 3) * n + j + 2] = c32; C[(i + 3) * n + j + 3] = c33; } } } /* -------------------------------------------------------------- harness */ typedef void (*mul_t)(const double *, const double *, double *, long); static double guard; /* clock() advances only at the 100 Hz scheduler tick, so a single small-N * multiply (sub-millisecond) cannot be timed directly: repeat it until the * batch spans at least 0.2 s (20 ticks), then divide. Returns seconds per * multiply. C ends holding the last multiply's result, same as any rep. */ static double time_mul(mul_t fn, const double *A, const double *B, double *C, long n) { long reps = 1, r; double t; clock_t t0; for (;;) { t0 = clock(); for (r = 0; r < reps; r++) { memset(C, 0, (size_t)(n * n) * sizeof(double)); fn(A, B, C, n); } t = (double)(clock() - t0) / (double)CLOCKS_PER_SEC; if (t >= 0.2 || reps >= 16384) break; reps *= 4; } guard += (double)reps; return t / (double)reps; } static double checksum(const double *C, long n) { double s = 0.0; long i; for (i = 0; i < n * n; i++) s += C[i]; guard += s; return s; } int main(int argc, char **argv) { static const long sizes[] = { 64, 128, 192, 256, 384, 512, 768, 1024 }; const char *tag = argc > 1 ? argv[1] : "run"; size_t si; printf("# tag,variant,n,seconds,mflop_s,checksum\n"); guard = 0.0; for (si = 0; si < sizeof(sizes) / sizeof(sizes[0]); si++) { long n = sizes[si]; double *A = (double *)malloc((size_t)(n * n) * sizeof(double)); double *B = (double *)malloc((size_t)(n * n) * sizeof(double)); double *C = (double *)malloc((size_t)(n * n) * sizeof(double)); double flop = 2.0 * (double)n * n * n; double t, t2, ref, got; int reps, r; if (A == NULL || B == NULL || C == NULL) { fprintf(stderr, "dgemm: out of memory at n=%ld\n", n); return 1; } fill(A, n); fill(B, n); /* ikj: the reference, best of 2 */ t = time_mul(mul_ikj, A, B, C, n); t2 = time_mul(mul_ikj, A, B, C, n); if (t2 < t) t = t2; ref = checksum(C, n); printf("%s,ikj,%ld,%.3f,%.1f,%.6f\n", tag, n, t, flop / t / 1e6, ref); /* blocked: best of 2 */ t = time_mul(mul_blocked, A, B, C, n); t2 = time_mul(mul_blocked, A, B, C, n); if (t2 < t) t = t2; got = checksum(C, n); printf("%s,blocked,%ld,%.3f,%.1f,%.6f%s\n", tag, n, t, flop / t / 1e6, got, fabs(got - ref) > 1e-9 * (fabs(ref) + 1.0) ? " # MISMATCH" : ""); /* reg4: register-blocked, the FPU's best case; best of 2 */ if (n % 4 != 0) { printf("%s,reg4,%ld,skipped\n", tag, n); } else { t = time_mul(mul_reg4, A, B, C, n); t2 = time_mul(mul_reg4, A, B, C, n); if (t2 < t) t = t2; got = checksum(C, n); printf("%s,reg4,%ld,%.3f,%.1f,%.6f%s\n", tag, n, t, flop / t / 1e6, got, fabs(got - ref) > 1e-9 * (fabs(ref) + 1.0) ? " # MISMATCH" : ""); } /* dot4: slowest and most memory-bound; once is enough for big N */ if (n > DOT4_MAX_N) { printf("%s,dot4,%ld,skipped\n", tag, n); } else { reps = n <= 384 ? 2 : 1; t = 0.0; for (r = 0; r < reps; r++) { t2 = time_mul(mul_dot4, A, B, C, n); if (t == 0.0 || t2 < t) t = t2; } got = checksum(C, n); printf("%s,dot4,%ld,%.3f,%.1f,%.6f%s\n", tag, n, t, flop / t / 1e6, got, fabs(got - ref) > 1e-9 * (fabs(ref) + 1.0) ? " # MISMATCH" : ""); } free(A); free(B); free(C); } printf("# guard %.17g\n", guard); return 0; }