1/*
2 * One run, described the same way src/bench/runSpec.ts describes it.
3 *
4 * This is the native half of the comparison: the same grid rule, the same
5 * presets and defaults, the same seeded perturbation, so that
6 *
7 * npm run bench -- --preset schnak-spots --lmax 63 --steps 2000
8 * ./shtbench --preset schnak-spots --lmax 63 --steps 2000
9 *
10 * describe the same computation, one through the WGSL transforms and one
11 * through upstream SHTNS.
12 *
13 * The tables below are a hand copy of src/mgpu/registry.ts and
14 * src/sht/layout.ts, which is unavoidable — C cannot import the TypeScript.
15 * It is also the one place the two sides could silently drift apart, so both
16 * emit their resolved spec in --json and scripts/compare-native.mjs refuses to
17 * compare two runs whose specs do not match.
18 *
19 * Compiled as C++ (by g++ for the CPU benchmark, by nvcc for the CUDA one), so
20 * the reaction can be written once as a template and used at float and double.
21 */
39/* ------------------------------------------------------------------ models */
41enum shtb_model_id { SHTB_SCHNAKENBERG = 0, SHTB_BRUSSELATOR, SHTB_ALLENCAHN };
43/* Every tunable scalar of every model. A model uses the subset its .m names as
44 * arguments; the tables below say which, in the .m's declared order. */
45struct shtb_params {
46 double a, b; /* schnakenberg */
47 double A, B; /* brusselator */
48 double eps2; /* allen-cahn */
49 double D1, D2; /* diffusivities (allen-cahn uses eps2 instead) */
50 double dt;
51};
53struct shtb_param_def {
54 const char *key;
55 size_t off; /* offset into shtb_params */
56 double dflt;
57};
62static const struct shtb_param_def SHTB_SCHNAK_PARAMS[] = {
63 SHTB_P(a, 0.1), SHTB_P(b, 0.9), SHTB_P(D1, 4e-4), SHTB_P(D2, 8e-3), SHTB_P(dt, 0.05),
64};
65static const struct shtb_param_def SHTB_BRUSSEL_PARAMS[] = {
66 SHTB_P(A, 3.0), SHTB_P(B, 9.0), SHTB_P(D1, 3.33e-3), SHTB_P(D2, 1.67e-2), SHTB_P(dt, 0.02),
67};
68static const struct shtb_param_def SHTB_ALLENCAHN_PARAMS[] = {
69 SHTB_P(eps2, 1e-3), SHTB_P(dt, 0.02),
70};
72struct shtb_model_def {
73 int id;
74 const char *key;
75 const char *label;
76 int nspecies;
77 int pdeg; /* polynomial degree of the reaction, for dealiasing */
78 double seed_amp; /* amplitude of the seeded perturbation */
79 const struct shtb_param_def *params;
80 int nparams;
81};
83static const struct shtb_model_def SHTB_MODELS[] = {
84 {SHTB_SCHNAKENBERG, "schnakenberg", "Schnakenberg", 2, 3, 1e-2, SHTB_SCHNAK_PARAMS, 5},
85 {SHTB_BRUSSELATOR, "brusselator", "Brusselator", 2, 3, 1e-2, SHTB_BRUSSEL_PARAMS, 5},
86 {SHTB_ALLENCAHN, "allencahn", "Allen-Cahn", 1, 3, 1e-2, SHTB_ALLENCAHN_PARAMS, 2},
87};
90struct shtb_override {
91 const char *key;
92 double value;
93};
95struct shtb_preset {
96 const char *key;
97 const char *label;
98 int model;
99 struct shtb_override over[2]; /* .key == NULL terminates */
100};
102static const struct shtb_preset SHTB_PRESETS[] = {
103 {"schnak-spots", "Schnakenberg - spots", SHTB_SCHNAKENBERG, {{NULL, 0}, {NULL, 0}}},
104 {"schnak-coarse", "Schnakenberg - coarse spots", SHTB_SCHNAKENBERG,
105 {{"D1", 1e-3}, {"D2", 2e-2}}},
106 {"schnak-fine", "Schnakenberg - fine spots", SHTB_SCHNAKENBERG,
107 {{"D1", 1.6e-4}, {"D2", 3.2e-3}}},
108 {"brussel", "Brusselator - stripes & spots", SHTB_BRUSSELATOR, {{NULL, 0}, {NULL, 0}}},
109 {"allencahn", "Allen-Cahn - coarsening", SHTB_ALLENCAHN, {{NULL, 0}, {NULL, 0}}},
110};
113static inline double *shtb_field(struct shtb_params *p, size_t off) {
114 return (double *)((char *)p + off);
115}
116static inline const double *shtb_field_c(const struct shtb_params *p, size_t off) {
117 return (const double *)((const char *)p + off);
118}
120/* --------------------------------------------------- the arithmetic per step */
122/*
123 * Everything a timestep needs, in the working precision. Products the .m forms
124 * from parameters — `(dt * D1)`, `(B + 1)` — are formed here in that same
125 * precision, so a float run multiplies floats exactly as the WGSL kernel does.
126 */
127template <typename real>
128struct shtb_step_const {
129 int model;
130 real p0, p1; /* (a, b) | (A, B) | unused */
131 real dt;
132 real dtD[2]; /* dt * D_k, one per species */
133};
135/*
136 * The reaction, transcribed from models/<key>.m. One line of MATLAB per line
137 * here, and `u.^3` is written out as `u*u*u` because that is what the WGSL
138 * backend emits for it (pow_i3; see emitPower in src/mgpu/wgsl.ts).
139 *
140 * schnakenberg: Un = (U + dt*analys(a - u + uuv)) ./ (1 + (dt*D1)*lam)
141 * Vn = (V + dt*analys(b - uuv)) ./ (1 + (dt*D2)*lam)
142 * brusselator: Un = (U + dt*analys(A - (B+1)*u + uuv)) ./ ...
143 * Vn = (V + dt*analys(B*u - uuv)) ./ ...
144 * allencahn: Un = (U + dt*analys(u - u.^3)) ./ (1 + (dt*eps2)*lam)
145 */
146template <typename real>
147SHTB_HD inline void shtb_react(const shtb_step_const<real> &c, real u, real v, real *r1,
148 real *r2) {
149 if (c.model == SHTB_ALLENCAHN) {
150 *r1 = u - u * u * u;
151 return;
152 }
153 const real uuv = u * u * v;
154 if (c.model == SHTB_SCHNAKENBERG) {
155 *r1 = c.p0 - u + uuv;
156 *r2 = c.p1 - uuv;
157 } else {
158 *r1 = c.p0 - (c.p1 + (real)1) * u + uuv;
159 *r2 = c.p1 * u - uuv;
160 }
161}
163/* One element of the IMEX update, for species k. `lam` is l(l+1) of that
164 * coefficient; the array is 2*nlm long with the value duplicated across the
165 * real and imaginary halves, matching the 2 x nlm spectral layout the .m sees. */
166template <typename real>
167SHTB_HD inline real shtb_imex(const shtb_step_const<real> &c, int k, real U, real R, real lam) {
168 return (U + c.dt * R) / ((real)1 + c.dtD[k] * lam);
169}
171/* ------------------------------------------------------------------ the spec */
173enum shtb_mode { SHTB_MODE_SOLVER = 0, SHTB_MODE_TRANSFORM };
174enum shtb_layout { SHTB_LAYOUT_THETA = 0, SHTB_LAYOUT_PHI };
176/* Defaults, from src/bench/runSpec.ts. */
183struct shtb_spec {
184 const struct shtb_preset *preset;
185 const struct shtb_model_def *model;
186 struct shtb_params params;
187 int lmax, nlat, nphi;
188 int seed, steps, warmup, batch;
189 int mode;
190 int layout;
191 int fp32; /* GPU only: use SHTNS' single-precision transforms */
192 int threads; /* CPU only: OpenMP threads, 0 = library default */
193 int json, digest;
194 const char *dump_state;
195};
197/* Grid sizes for a given lmax, dealiased for a reaction of polynomial degree
198 * pdeg. Identical to gridForLmax() in src/sht/layout.ts, including rounding
199 * nphi up to a power of two — which the WGSL side needs for its FFT path and
200 * which SHTNS does not, but the grids have to match to compare anything. */
201static inline void shtb_grid_for_lmax(int lmax, int pdeg, int *nlat, int *nphi) {
202 double min_lat = ((double)(pdeg + 1) * lmax + 1) / 2.0;
203 if (min_lat < lmax + 1) min_lat = lmax + 1;
204 *nlat = 2 * (int)ceil(min_lat / 2.0);
205 int n = 1;
206 while (n < (pdeg + 1) * lmax + 1) n *= 2;
207 *nphi = n;
208}
210static inline long shtb_nlm_calc(int lmax, int mmax) {
211 return (long)(mmax + 1) * (lmax + 1) - (long)mmax * (mmax + 1) / 2;
212}
214/* -------------------------------------------------------------- the seeding */
216/*
217 * mulberry32 + Box-Muller, transcribed from src/mgpu/noise.ts so an integer
218 * seed means the same perturbation on both sides. The integer state evolves
219 * bit-identically; the Box-Muller step then goes through log/sqrt/sin/cos, so a
220 * libm that rounds differently from V8's can differ in the last bit. Both sides
221 * report the perturbation's RMS for that reason.
222 */
223struct shtb_rng {
224 uint32_t s;
225 int have_spare;
226 double spare;
227};
229static inline void shtb_rng_init(struct shtb_rng *r, uint32_t seed) {
230 r->s = seed;
231 r->have_spare = 0;
232 r->spare = 0;
233}
235/* uniform in [0, 1) */
236static inline double shtb_rand(struct shtb_rng *r) {
237 r->s = r->s + 0x6d2b79f5u;
238 uint32_t t = r->s;
239 t = (t ^ (t >> 15)) * (t | 1u);
240 t ^= t + (t ^ (t >> 7)) * (t | 61u);
241 return (double)(t ^ (t >> 14)) / 4294967296.0;
242}
244static inline double shtb_randn(struct shtb_rng *r) {
245 if (r->have_spare) {
246 r->have_spare = 0;
247 return r->spare;
248 }
249 double u = 0;
250 while (u == 0) u = shtb_rand(r);
251 double rad = sqrt(-2 * log(u));
252 double th = 2 * M_PI * shtb_rand(r);
253 r->spare = rad * sin(th);
254 r->have_spare = 1;
255 return rad * cos(th);
256}
258/* amp-scaled normal deviates, one per grid point, in [ilat*nphi + iphi] order —
259 * the order src/mgpu/noise.ts produces them in. Rounded to float, because the
260 * browser stores them in a Float32Array; the fp64 run then differs from the
261 * fp32 one only in the arithmetic, not in the initial condition. */
262static inline void shtb_seeded_noise(long npts, double amp, uint32_t seed, float *out) {
263 struct shtb_rng r;
264 shtb_rng_init(&r, seed);
265 for (long i = 0; i < npts; i++) out[i] = (float)(amp * shtb_randn(&r));
266}
268/*
269 * A seeded spectrum for the transform-only benchmark: uniform in [-1, 1), and
270 * bit-identical to the TypeScript side because it never leaves integer
271 * arithmetic and exactly-representable doubles. The m = 0 imaginary parts are
272 * zeroed, since a real field has none and the two libraries need not agree on
273 * what to do with a coefficient that cannot occur.
274 *
275 * qlm is interleaved [re, im] per coefficient, SHTNS LM ordering.
276 */
277static inline void shtb_seeded_spectrum(int lmax, int mmax, uint32_t seed, float *qlm) {
278 struct shtb_rng r;
279 shtb_rng_init(&r, seed);
280 long lm = 0;
281 for (int m = 0; m <= mmax; m++) {
282 for (int l = m; l <= lmax; l++, lm++) {
283 qlm[2 * lm] = (float)(2 * shtb_rand(&r) - 1);
284 float im = (float)(2 * shtb_rand(&r) - 1);
285 qlm[2 * lm + 1] = (m == 0) ? 0.0f : im;
286 }
287 }
288}
290/* -------------------------------------------------------------- statistics */
292static inline double shtb_now_ms(void) {
293 struct timespec ts;
294 clock_gettime(CLOCK_MONOTONIC, &ts);
295 return (double)ts.tv_sec * 1e3 + (double)ts.tv_nsec * 1e-6;
296}
298struct shtb_timing {
299 double mean_ms, median_ms, p05_ms, p95_ms, min_ms;
300};
302static int shtb_cmp_double(const void *a, const void *b) {
303 double x = *(const double *)a, y = *(const double *)b;
304 return (x > y) - (x < y);
305}
307static inline struct shtb_timing shtb_stats(double *samples, int n) {
308 qsort(samples, (size_t)n, sizeof(double), shtb_cmp_double);
309 double total = 0;
310 for (int i = 0; i < n; i++) total += samples[i];
311 int i50 = (int)(0.50 * n), i05 = (int)(0.05 * n), i95 = (int)(0.95 * n);
312 if (i50 >= n) i50 = n - 1;
313 if (i05 >= n) i05 = n - 1;
314 if (i95 >= n) i95 = n - 1;
315 struct shtb_timing t;
316 t.mean_ms = total / n;
317 t.median_ms = samples[i50];
318 t.p05_ms = samples[i05];
319 t.p95_ms = samples[i95];
320 t.min_ms = samples[0];
321 return t;
322}
324/* The same five numbers digestOf() computes in src/mgpu/digest.ts. */
325struct shtb_digest {
326 long n;
327 double min, max, mean, rms;
328};
330static inline struct shtb_digest shtb_digest_of(const float *v, long n) {
331 struct shtb_digest d;
332 d.n = n;
333 d.min = INFINITY;
334 d.max = -INFINITY;
335 double sum = 0, sumsq = 0;
336 for (long i = 0; i < n; i++) {
337 double x = v[i];
338 if (x < d.min) d.min = x;
339 if (x > d.max) d.max = x;
340 sum += x;
341 sumsq += x * x;
342 }
343 d.mean = sum / (double)n;
344 d.rms = sqrt(sumsq / (double)n);
345 return d;
346}
348static inline int shtb_all_finite(const float *v, long n) {
349 for (long i = 0; i < n; i++)
350 if (!isfinite(v[i])) return 0;
351 return 1;
352}
354/* ----------------------------------------------------------- spec resolution */
356static inline const struct shtb_model_def *shtb_model_by_id(int id) {
357 for (int i = 0; i < SHTB_NMODELS; i++)
358 if (SHTB_MODELS[i].id == id) return &SHTB_MODELS[i];
359 return NULL;
360}
362static inline void shtb_default_params(const struct shtb_model_def *m, struct shtb_params *p) {
363 memset(p, 0, sizeof(*p));
364 for (int i = 0; i < m->nparams; i++) *shtb_field(p, m->params[i].off) = m->params[i].dflt;
365}
367/* Homogeneous background each species starts from, and the diffusivity each is
368 * advanced with. From the init/step functions of models/<key>.m; only species 0
369 * gets the seeded perturbation. */
370static inline void shtb_background(const struct shtb_model_def *m, const struct shtb_params *p,
371 double base[2]) {
372 base[0] = base[1] = 0;
373 if (m->id == SHTB_SCHNAKENBERG) {
374 double us = p->a + p->b;
375 base[0] = us;
376 base[1] = p->b / (us * us);
377 } else if (m->id == SHTB_BRUSSELATOR) {
378 base[0] = p->A;
379 base[1] = p->B / p->A;
380 }
381}
383static inline void shtb_diffusivity(const struct shtb_model_def *m, const struct shtb_params *p,
384 double d[2]) {
385 if (m->id == SHTB_ALLENCAHN) {
386 d[0] = p->eps2;
387 d[1] = 0;
388 } else {
389 d[0] = p->D1;
390 d[1] = p->D2;
391 }
392}
394template <typename real>
395static inline shtb_step_const<real> shtb_make_step_const(const struct shtb_model_def *m,
396 const struct shtb_params *p) {
397 double d[2];
398 shtb_diffusivity(m, p, d);
399 shtb_step_const<real> c;
400 c.model = m->id;
401 c.p0 = (real)(m->id == SHTB_BRUSSELATOR ? p->A : p->a);
402 c.p1 = (real)(m->id == SHTB_BRUSSELATOR ? p->B : p->b);
403 c.dt = (real)p->dt;
404 /* (dt * D_k) as one product in the working precision, as the .m writes it */
405 c.dtD[0] = (real)p->dt * (real)d[0];
406 c.dtD[1] = (real)p->dt * (real)d[1];
407 return c;
408}
410/* ------------------------------------------------------------ argument parsing
411 *
412 * `--key value` or `--key=value`, in any order — the same grammar parseArgs()
413 * accepts in src/bench/runSpec.ts, plus the flags only a native run has.
414 */
416static inline int shtb_parse_spec(int argc, char **argv, struct shtb_spec *s, const char *usage) {
417 const struct shtb_preset *preset = &SHTB_PRESETS[0];
419 /* --preset first: it decides which parameter names are legal. */
420 for (int i = 1; i < argc; i++) {
421 const char *a = argv[i];
422 const char *v = NULL;
423 if (strcmp(a, "--preset") == 0 && i + 1 < argc)
424 v = argv[i + 1];
425 else if (strncmp(a, "--preset=", 9) == 0)
426 v = a + 9;
427 if (!v) continue;
428 const struct shtb_preset *found = NULL;
429 for (int k = 0; k < SHTB_NPRESETS; k++)
430 if (strcmp(SHTB_PRESETS[k].key, v) == 0) found = &SHTB_PRESETS[k];
431 if (!found) {
432 fprintf(stderr, "shtbench: unknown preset '%s' (have:", v);
433 for (int k = 0; k < SHTB_NPRESETS; k++) fprintf(stderr, " %s", SHTB_PRESETS[k].key);
434 fprintf(stderr, ")\n");
435 return 2;
436 }
437 preset = found;
438 }
440 memset(s, 0, sizeof(*s));
441 s->preset = preset;
442 s->model = shtb_model_by_id(preset->model);
443 shtb_default_params(s->model, &s->params);
444 for (int k = 0; k < 2; k++)
445 if (preset->over[k].key) {
446 for (int i = 0; i < s->model->nparams; i++)
447 if (strcmp(s->model->params[i].key, preset->over[k].key) == 0)
448 *shtb_field(&s->params, s->model->params[i].off) = preset->over[k].value;
449 }
450 s->lmax = SHTB_DEFAULT_LMAX;
451 s->seed = SHTB_DEFAULT_SEED;
452 s->steps = SHTB_DEFAULT_STEPS;
453 s->warmup = SHTB_DEFAULT_WARMUP;
454 s->batch = SHTB_DEFAULT_BATCH;
455 s->mode = SHTB_MODE_SOLVER;
456 s->layout = SHTB_LAYOUT_THETA;
457 s->fp32 = 1;
458 s->threads = 0;
460 for (int i = 1; i < argc; i++) {
461 const char *a = argv[i];
462 if (strncmp(a, "--", 2) != 0) {
463 fprintf(stderr, "shtbench: unexpected argument '%s'\n\n%s\n", a, usage);
464 return 2;
465 }
466 if (strcmp(a, "--help") == 0 || strcmp(a, "-h") == 0) {
467 printf("%s\n", usage);
468 return 1;
469 }
470 if (strcmp(a, "--json") == 0) {
471 s->json = 1;
472 continue;
473 }
474 if (strcmp(a, "--digest") == 0) {
475 s->digest = 1;
476 continue;
477 }
478 if (strcmp(a, "--fp64") == 0) {
479 s->fp32 = 0;
480 continue;
481 }
483 /* split key / value */
484 char key[64];
485 const char *val = NULL;
486 const char *eq = strchr(a, '=');
487 if (eq) {
488 size_t n = (size_t)(eq - a - 2);
489 if (n >= sizeof(key)) n = sizeof(key) - 1;
490 memcpy(key, a + 2, n);
491 key[n] = 0;
492 val = eq + 1;
493 } else {
494 snprintf(key, sizeof(key), "%s", a + 2);
495 if (i + 1 >= argc) {
496 fprintf(stderr, "shtbench: --%s needs a value\n", key);
497 return 2;
498 }
499 val = argv[++i];
500 }
502 if (strcmp(key, "preset") == 0) continue; /* handled above */
503 if (strcmp(key, "lmax") == 0) {
504 s->lmax = atoi(val);
505 continue;
506 }
507 if (strcmp(key, "seed") == 0) {
508 s->seed = atoi(val);
509 continue;
510 }
511 if (strcmp(key, "steps") == 0) {
512 s->steps = atoi(val);
513 continue;
514 }
515 if (strcmp(key, "warmup") == 0) {
516 s->warmup = atoi(val);
517 continue;
518 }
519 if (strcmp(key, "batch") == 0) {
520 s->batch = atoi(val);
521 continue;
522 }
523 if (strcmp(key, "threads") == 0) {
524 s->threads = atoi(val);
525 continue;
526 }
527 if (strcmp(key, "dump-state") == 0) {
528 s->dump_state = val;
529 s->digest = 1;
530 continue;
531 }
532 if (strcmp(key, "mode") == 0) {
533 if (strcmp(val, "solver") == 0)
534 s->mode = SHTB_MODE_SOLVER;
535 else if (strcmp(val, "transform") == 0)
536 s->mode = SHTB_MODE_TRANSFORM;
537 else {
538 fprintf(stderr, "shtbench: --mode must be 'solver' or 'transform' (got '%s')\n", val);
539 return 2;
540 }
541 continue;
542 }
543 if (strcmp(key, "layout") == 0) {
544 if (strcmp(val, "theta") == 0)
545 s->layout = SHTB_LAYOUT_THETA;
546 else if (strcmp(val, "phi") == 0)
547 s->layout = SHTB_LAYOUT_PHI;
548 else {
549 fprintf(stderr, "shtbench: --layout must be 'theta' or 'phi' (got '%s')\n", val);
550 return 2;
551 }
552 continue;
553 }
555 int matched = 0;
556 for (int p = 0; p < s->model->nparams; p++)
557 if (strcmp(s->model->params[p].key, key) == 0) {
558 *shtb_field(&s->params, s->model->params[p].off) = atof(val);
559 matched = 1;
560 }
561 if (!matched) {
562 fprintf(stderr, "shtbench: unknown option --%s\nparameters of %s:", key, s->model->label);
563 for (int p = 0; p < s->model->nparams; p++)
564 fprintf(stderr, " --%s", s->model->params[p].key);
565 fprintf(stderr, "\n");
566 return 2;
567 }
568 }
570 if (s->lmax < 1) {
571 fprintf(stderr, "shtbench: --lmax must be >= 1\n");
572 return 2;
573 }
574 if (s->steps < 1 || s->warmup < 0 || s->batch < 1) {
575 fprintf(stderr, "shtbench: --steps and --batch must be >= 1, --warmup >= 0\n");
576 return 2;
577 }
578 shtb_grid_for_lmax(s->lmax, s->model->pdeg, &s->nlat, &s->nphi);
579 return 0;
580}
582/* ------------------------------------------------------------- JSON emission
583 *
584 * Deliberately shaped like the object scripts/bench.ts prints with --json, so
585 * scripts/compare-native.mjs can read a WebGPU run and a native run the same
586 * way. "step" means one solver timestep in solver mode and one
587 * spectral->grid->spectral round trip in transform mode.
588 */
590/* Strings from the library go into JSON, so strip anything that would break it.
591 * Returns `dst`, for use inline. */
592static inline char *shtb_json_safe(char *dst, size_t cap, const char *src) {
593 size_t j = 0;
594 for (size_t i = 0; src && src[i] && j + 1 < cap; i++) {
595 unsigned char c = (unsigned char)src[i];
596 dst[j++] = (c < 0x20 || c == '"' || c == '\\' || c == 0x7f) ? ' ' : (char)c;
597 }
598 dst[j] = 0;
599 return dst;
600}
602struct shtb_report {
603 const char *library; /* e.g. "SHTNS 3.7.5" */
604 const char *runtime; /* e.g. "cuda 12.4, vkfft" */
605 const char *adapter; /* GPU name, or the CPU's thread count */
606 const char *precision;/* "fp32" | "fp64" */
607 const char *fourier; /* which FFT the library used */
608 long nlm;
609 int ops_per_step;
610 double ms_per_step;
611 double encode_ms_per_step;
612 struct shtb_timing latency;
613 int have_latency;
614 struct shtb_digest digest;
615 int have_digest;
616 struct shtb_digest input_digest; /* transform mode: the seeded spectrum */
617 int have_input_digest;
618 double field_min, field_max;
619 int finite;
620 double model_t;
621 int steps_run;
622};
624static inline void shtb_print_json(const struct shtb_spec *s, const struct shtb_report *r) {
625 printf("{\n");
626 printf(" \"mode\": \"%s\",\n", s->mode == SHTB_MODE_SOLVER ? "solver" : "transform");
627 printf(" \"spec\": {\n");
628 printf(" \"preset\": \"%s\",\n", s->preset->key);
629 printf(" \"lmax\": %d,\n", s->lmax);
630 printf(" \"seed\": %d,\n", s->seed);
631 printf(" \"steps\": %d,\n", s->steps);
632 printf(" \"warmup\": %d,\n", s->warmup);
633 printf(" \"params\": {");
634 for (int i = 0; i < s->model->nparams; i++)
635 printf("%s\"%s\": %.17g", i ? ", " : "", s->model->params[i].key,
636 *shtb_field_c(&s->params, s->model->params[i].off));
637 printf("}\n },\n");
638 printf(" \"model\": \"%s\",\n", s->model->key);
639 printf(" \"backend\": {\"library\": \"%s\", \"runtime\": \"%s\", \"adapter\": \"%s\", "
640 "\"precision\": \"%s\", \"layout\": \"%s\"},\n",
641 r->library, r->runtime, r->adapter, r->precision,
642 s->layout == SHTB_LAYOUT_PHI ? "phi-contiguous" : "theta-contiguous");
643 printf(" \"grid\": {\"lmax\": %d, \"nlat\": %d, \"nphi\": %d, \"nlm\": %ld},\n", s->lmax,
644 s->nlat, s->nphi, r->nlm);
645 printf(" \"compiled\": {\"opsPerStep\": %d},\n", r->ops_per_step);
646 printf(" \"throughput\": {\"batch\": %d, \"msPerStep\": %.17g, \"stepsPerSec\": %.17g, "
647 "\"encodeMsPerStep\": %.17g},\n",
648 s->batch, r->ms_per_step, 1000.0 / r->ms_per_step, r->encode_ms_per_step);
649 if (r->have_latency)
650 printf(" \"latency\": {\"meanMs\": %.17g, \"medianMs\": %.17g, \"p05Ms\": %.17g, "
651 "\"p95Ms\": %.17g, \"minMs\": %.17g},\n",
652 r->latency.mean_ms, r->latency.median_ms, r->latency.p05_ms, r->latency.p95_ms,
653 r->latency.min_ms);
654 else
655 printf(" \"latency\": null,\n");
656 if (r->have_digest)
657 printf(" \"digest\": {\"n\": %ld, \"min\": %.17g, \"max\": %.17g, \"mean\": %.17g, "
658 "\"rms\": %.17g, \"fourier\": \"%s\", \"adapter\": \"%s\"},\n",
659 r->digest.n, r->digest.min, r->digest.max, r->digest.mean, r->digest.rms, r->fourier,
660 r->adapter);
661 else
662 printf(" \"digest\": null,\n");
663 if (r->have_input_digest)
664 printf(" \"input\": {\"n\": %ld, \"min\": %.17g, \"max\": %.17g, \"mean\": %.17g, "
665 "\"rms\": %.17g},\n",
666 r->input_digest.n, r->input_digest.min, r->input_digest.max, r->input_digest.mean,
667 r->input_digest.rms);
668 else
669 printf(" \"input\": null,\n");
670 printf(" \"state\": {\"t\": %.17g, \"steps\": %d, \"min\": %.17g, \"max\": %.17g, "
671 "\"contrast\": %.17g, \"finite\": %s}\n",
672 r->model_t, r->steps_run, r->field_min, r->field_max, r->field_max - r->field_min,
673 r->finite ? "true" : "false");
674 printf("}\n");
675}
677/* The state file scripts/compare-native.mjs diffs, in the same shape
678 * `npm run bench -- --dump-state` writes. */
679static inline int shtb_dump_state(const char *path, const struct shtb_spec *s,
680 const struct shtb_report *r, const float *state, long n) {
681 FILE *f = fopen(path, "w");
682 if (!f) return -1;
683 fprintf(f, "{\"spec\":{\"preset\":\"%s\",\"lmax\":%d,\"seed\":%d,\"steps\":%d,\"warmup\":%d,"
684 "\"params\":{",
685 s->preset->key, s->lmax, s->seed, s->steps, s->warmup);
686 for (int i = 0; i < s->model->nparams; i++)
687 fprintf(f, "%s\"%s\":%.17g", i ? "," : "", s->model->params[i].key,
688 *shtb_field_c(&s->params, s->model->params[i].off));
689 fprintf(f, "}},\"mode\":\"%s\",\"backend\":{\"library\":\"%s\",\"adapter\":\"%s\","
690 "\"precision\":\"%s\"},",
691 s->mode == SHTB_MODE_SOLVER ? "solver" : "transform", r->library, r->adapter,
692 r->precision);
693 fprintf(f, "\"digest\":{\"n\":%ld,\"min\":%.17g,\"max\":%.17g,\"mean\":%.17g,\"rms\":%.17g,"
694 "\"fourier\":\"%s\",\"adapter\":\"%s\"},",
695 r->digest.n, r->digest.min, r->digest.max, r->digest.mean, r->digest.rms, r->fourier,
696 r->adapter);
697 if (r->have_input_digest)
698 fprintf(f, "\"input\":{\"n\":%ld,\"min\":%.17g,\"max\":%.17g,\"mean\":%.17g,\"rms\":%.17g},",
699 r->input_digest.n, r->input_digest.min, r->input_digest.max, r->input_digest.mean,
700 r->input_digest.rms);
701 fprintf(f, "\"state\":[");
702 for (long i = 0; i < n; i++) fprintf(f, "%s%.9g", i ? "," : "", (double)state[i]);
703 fprintf(f, "]}\n");
704 fclose(f);
705 return 0;
706}