1/* The gemm native program's inputs and reference (mesh-ref lib.rs `gemm`). */
2#include "mesh/ref/gemm.h"
3
4#include <math.h>
5#include <stdlib.h>
6
7float mesh_gemm_a(size_t i, size_t k) { return (float)((i * 7 + k * 13) % 17) / 8.0f - 1.0f; }
8float mesh_gemm_b(size_t k, size_t j) { return (float)((k * 5 + j * 11) % 19) / 9.0f - 1.0f; }
9
10void mesh_gemm_inputs(size_t n, float *a, float *b) {
11    for (size_t x = 0; x < n * n; x++) {
12        a[x] = mesh_gemm_a(x / n, x % n);
13        b[x] = mesh_gemm_b(x / n, x % n);
14    }
15}
16
17void mesh_gemm_reference(size_t n, double *c, double *s) {
18    for (size_t x = 0; x < n * n; x++) c[x] = s[x] = 0.0;
19    for (size_t i = 0; i < n; i++)
20        for (size_t j = 0; j < n; j++)
21            for (size_t k = 0; k < n; k++) {
22                double p = (double)mesh_gemm_a(i, k) * (double)mesh_gemm_b(k, j);
23                c[i * n + j] += p;
24                s[i * n + j] += fabs(p);
25            }
26}
27
28double mesh_gemm_max_rel_err(size_t n, const float *got) {
29    double *c = malloc(n * n * sizeof *c), *s = malloc(n * n * sizeof *s);
30    if (!c || !s) {
31        free(c), free(s);
32        return nan("");
33    }
34    mesh_gemm_reference(n, c, s);
35    double worst = 0.0;
36    for (size_t x = 0; x < n * n; x++) {
37        double e = fabs((double)got[x] - c[x]) / (s[x] + 1e-30);
38        worst = fmax(worst, e); /* f64::max: a NaN error is ignored, as in the Rust fold */
39    }
40    free(c), free(s);
41    return worst;
42}