1/* mesh-ref matmul.rs's test, then a golden digest from the Rust: matmul and matvec at four shapes and
2 * a thousand values. */
3#include "mesh/ref/matmul.h"
4
5#include "mesh/ftz.h"
6
7#include <stdio.h>
8#include <stdlib.h>
9
10static int failures;
11#define CHECK(c)                                                                                   \
12    do {                                                                                           \
13        if (!(c)) {                                                                                \
14            fprintf(stderr, "%s:%d: CHECK(%s) failed\n", __FILE__, __LINE__, #c);                 \
15            failures++;                                                                            \
16        }                                                                                          \
17    } while (0)
18
19static uint64_t fnv(uint64_t h, uint64_t w) { return (h ^ w) * 0x100000001b3ull; }
20
21static void rust_exact_and_ordered(void) {
22    size_t m = 5, k = 37, n = 3;
23    float a[5 * 37], b[37 * 3], c[15], c1[15], y[5];
24    mesh_matmul_matrix(1, m * k, a);
25    mesh_matmul_matrix(2, k * n, b);
26    for (size_t i = 0; i < m * k; i++) CHECK(a[i] >= -1.0f && a[i] < 1.0f);
27    mesh_matmul(a, b, m, k, n, 3, c);
28    mesh_matmul(a, b, m, k, n, 1, c1);
29    for (size_t i = 0; i < m * n; i++) CHECK(mesh_ftz_bits(c[i]) == mesh_ftz_bits(c1[i]));
30    /* The sequential order, as plain f32 (no subnormals here: FTZ changes nothing). */
31    float acc = 0.0f;
32    for (size_t kk = 0; kk < k; kk++) acc += a[2 * k + kk] * b[kk * n + 1];
33    CHECK(mesh_ftz_bits(c[2 * n + 1]) == mesh_ftz_bits(acc));
34    mesh_matvec(a, b, m, k, y);
35}
36
37static void rust_golden(void) {
38    const size_t shapes[4][3] = {{5, 37, 3}, {16, 64, 16}, {33, 100, 7}, {64, 64, 64}};
39    uint64_t h = 0xcbf29ce484222325ull;
40    for (int s = 0; s < 4; s++) {
41        size_t m = shapes[s][0], k = shapes[s][1], n = shapes[s][2];
42        float *a = malloc(m * k * sizeof *a), *b = malloc(k * n * sizeof *b), *c = malloc(m * n * sizeof *c),
43              *y = malloc(m * sizeof *y);
44        mesh_matmul_matrix(1, m * k, a);
45        mesh_matmul_matrix(2, k * n, b);
46        mesh_matmul(a, b, m, k, n, 3, c);
47        h = fnv(h, mesh_matmul_digest(c, m * n));
48        mesh_matvec(a, b, m, k, y);
49        for (size_t i = 0; i < m; i++) h = fnv(h, mesh_ftz_bits(y[i]));
50        free(a), free(b), free(c), free(y);
51    }
52    for (uint32_t i = 0; i < 1000; i++) h = fnv(h, mesh_ftz_bits(mesh_matmul_value(i * 7919, i)));
53    if (h != 0x03c20f795c51b1ecull) {
54        fprintf(stderr, "matmul digest %#018llx, want 0x03c20f795c51b1ec\n", (unsigned long long)h);
55        failures++;
56    }
57}
58
59int main(void) {
60    rust_exact_and_ordered();
61    rust_golden();
62    return failures != 0;
63}