David
/
Foundry
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Foundry/cuda/tests/bench_matvec.c

948 lines
38 KiB
C
Raw Permalink Blame History

/* SPDX-License-Identifier: AGPL-3.0-or-later */
/* Copyright (C) 2026 David (Corey) Herrera */
/*
* bench_matvec.c — micro-benchmark + bit-exact verifier for the Foundry quant
* matvec + grouped MoE kernels on one K80 die. No model load: weights are
* random bytes (fp16 block scales clamped to a sane exponent so sums stay
* finite), shapes mirror the Nex-N2-mini hot paths. Reports ms/launch and
* effective GB/s of weight-bytes so kernel variants can be compared in
* seconds instead of a 100 s model load per iteration.
*
* VERIFY: for each quant type a CPU replica recomputes the kernel output with
* the exact device op order (per-lane fmaf accumulation chains + the 5-round
* __shfl_down tree, mirroring the PTX: q3k/q6k partial = fma(acc, d, partial);
* q5k term = (d*sc)*nib - dmin*mn, partial += acc; q80 partial =
* fma(x, d*q, partial)). The kernel result must match BIT-EXACTLY — this is
* how the upload-time repack (fc_repack_weight) is proven numerically
* identical to the original GGUF-layout path.
*
* fc_matvec_q4k is the exception: it is MIN-LIFTED (see matvec_kquants.cu),
* which reassociates the sum, so it cannot be bit-identical to the previous
* term order. verify_mv_reassoc gates it BIT-EXACTLY against a replica of
* its new order (ref_row_q4k_minlift) and additionally bounds its fp32
* relative error against an fp64 reference, both absolutely and relative to
* the old order's error. The other quants keep the zero-tolerance memcmp.
*
* Build (from foundry/):
* gcc -O2 -I. -o cuda/tests/bench_matvec cuda/tests/bench_matvec.c \
* cuda/fc_cuda.c -lm -ldl
* -DUSE_REPACK=0 benches/verifies the raw-GGUF-layout kernels (the
* pre-repack ones); default 1 repacks via fc_repack_weight first.
* Run: CUDA_VISIBLE_DEVICES=0 ./cuda/tests/bench_matvec
*/
#include "cuda/fc_cuda.h"
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifndef USE_REPACK
#define USE_REPACK 1
#endif
static double now_s(void)
{
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9;
}
/* ------------------------------------------------------------------
* Random weight/activation generation
* ------------------------------------------------------------------ */
static size_t raw_block_bytes(int t)
{
switch (t) {
case 8: return 34;
case 11: return 110;
case 12: return 144;
case 13: return 176;
case 14: return 210;
default: return 0;
}
}
/* Clamp every fp16 block scale (d, and dmin for Q4_K/Q5_K) to exponent
* 2^-4..2^-3 with a random mantissa: random halves can be inf/nan which
* breaks the bit-exact compare (nan != nan). Positions per gguf_dequant.h. */
static void sanitize_quant(int t, unsigned char *p, size_t bytes)
{
size_t bs = raw_block_bytes(t);
size_t off;
for (off = 0; bs && off + bs <= bytes; off += bs) {
unsigned char *b = p + off;
unsigned short h = (unsigned short)(0x2C00 | (rand() & 0x3FF));
unsigned short h2 = (unsigned short)(0x2C00 | (rand() & 0x3FF));
switch (t) {
case 11: b[108] = h & 0xFF; b[109] = h >> 8; break;
case 14: b[208] = h & 0xFF; b[209] = h >> 8; break;
case 12: case 13:
b[0] = h & 0xFF; b[1] = h >> 8;
b[2] = h2 & 0xFF; b[3] = h2 >> 8;
break;
case 8: b[0] = h & 0xFF; b[1] = h >> 8; break;
}
}
}
static unsigned char *rand_bytes(size_t bytes)
{
unsigned char *p = malloc(bytes);
size_t i;
for (i = 0; i < bytes; i++)
p[i] = (unsigned char)(rand() & 0xFF);
return p;
}
static float *rand_floats(size_t n)
{
float *p = malloc(n * sizeof(float));
size_t i;
for (i = 0; i < n; i++)
p[i] = (float)rand() / (float)RAND_MAX - 0.5f;
return p;
}
/* Upload a raw quant weight for the kernels: repacked (default) or as-is.
* Returns the device buffer and writes the device row stride. */
static fc_buffer *upload_weight(int t, const unsigned char *raw, long rows,
int cols, int *dev_stride)
{
size_t raw_stride = raw_block_bytes(t) * (size_t)((t == 8) ? cols / 32
: cols / 256);
#if USE_REPACK
size_t rs = fc_repack_row_bytes(t, cols);
unsigned char *tmp = malloc(rs * (size_t)rows);
fc_buffer *b;
if (!tmp || fc_repack_weight(t, raw, tmp, rows, cols) != 0) {
free(tmp);
return NULL;
}
b = fc_buffer_from_ptr(tmp, rs * (size_t)rows);
free(tmp);
*dev_stride = (int)rs;
return b;
#else
*dev_stride = (int)raw_stride;
return fc_buffer_from_ptr(raw, raw_stride * (size_t)rows);
#endif
}
/* ------------------------------------------------------------------
* CPU replica — decodes the ORIGINAL GGUF layout (source of truth for
* both kernel generations) with the device's exact fp op order.
* ------------------------------------------------------------------ */
static float half2f(unsigned short h)
{
unsigned int sign = ((unsigned int)(h & 0x8000)) << 16;
int exp = (h >> 10) & 0x1F;
unsigned int mant = h & 0x3FF;
unsigned int bits;
float out;
if (exp == 0) {
if (mant == 0) {
bits = sign;
} else {
while (!(mant & 0x400)) { mant <<= 1; exp--; }
exp++;
mant &= ~0x400u;
bits = sign | ((unsigned int)(exp + 127 - 15) << 23) | (mant << 13);
}
} else if (exp == 31) {
bits = sign | 0x7F800000u | (mant << 13);
} else {
bits = sign | ((unsigned int)(exp + 127 - 15) << 23) | (mant << 13);
}
memcpy(&out, &bits, sizeof(out));
return out;
}
/* 5-round __shfl_down(16,8,4,2,1) + add tree over the 32 lane partials
* (out-of-range shfl returns the caller's own value). Result = lane 0. */
static float tree32(float v[32])
{
float t[32];
int off, l;
for (off = 16; off >= 1; off >>= 1) {
for (l = 0; l < 32; l++)
t[l] = v[l] + ((l + off < 32) ? v[l + off] : v[l]);
memcpy(v, t, sizeof(t));
}
return v[0];
}
static void kq_scale_min(const unsigned char *sr, int s, float *sc, float *mn)
{
if (s < 4) {
*sc = (float)(sr[s] & 63);
*mn = (float)(sr[s + 4] & 63);
} else {
*sc = (float)((sr[s + 4] & 0xF) | ((sr[s - 4] >> 6) << 4));
*mn = (float)((sr[s + 4] >> 4) | ((sr[s] >> 6) << 4));
}
}
/* One output row of each kernel, replicated per-lane. The _lanes variants
* accumulate the 32 per-lane fmaf chains into p[] (caller zeroes) without
* the final shfl tree — the grouped-MoE replicas below need the raw lane
* partials to reproduce the down kernels' per-slot weighted chains. */
static void ref_lanes_q3k(const unsigned char *row, const float *x, int cols,
float p[32])
{
int n_sb = cols >> 8, sb, lane, n, j;
for (sb = 0; sb < n_sb; sb++) {
const unsigned char *blk = row + (size_t)sb * 110;
const unsigned char *hm = blk, *qs = blk + 32, *sc = blk + 96;
float d = half2f((unsigned short)(blk[108] | (blk[109] << 8)));
const float *xb = x + (sb << 8);
for (lane = 0; lane < 32; lane++) {
int lhalf = lane >> 4;
float acc = 0.0f;
for (n = 0; n < 2; n++) {
unsigned char qb = qs[32 * n + lane];
for (j = 0; j < 4; j++) {
int s = 8 * n + 2 * j + lhalf;
int i = s & 3, a = s >> 2;
unsigned char lowb = (a & 1) ? sc[i + 4] : sc[i];
int low = (a & 2) ? (lowb >> 4) : (lowb & 0xF);
int high = (sc[8 + i] >> (2 * a)) & 3;
int scale = (low | (high << 4)) - 32;
int qv = (qb >> (2 * j)) & 3;
int hb = (hm[lane] & (1u << (4 * n + j))) ? 0 : 4;
acc = fmaf(xb[128 * n + 32 * j + lane],
(float)(scale * (qv - hb)), acc);
}
}
p[lane] = fmaf(acc, d, p[lane]);
}
}
}
static float ref_row_q3k(const unsigned char *row, const float *x, int cols)
{
float p[32] = { 0 };
ref_lanes_q3k(row, x, cols, p);
return tree32(p);
}
static void ref_lanes_q4k(const unsigned char *row, const float *x, int cols,
float p[32])
{
int n_sb = cols >> 8, sb, lane, j, h;
for (sb = 0; sb < n_sb; sb++) {
const unsigned char *blk = row + (size_t)sb * 144;
float d = half2f((unsigned short)(blk[0] | (blk[1] << 8)));
float dmin = half2f((unsigned short)(blk[2] | (blk[3] << 8)));
const float *xb = x + (sb << 8);
for (lane = 0; lane < 32; lane++) {
float acc = 0.0f;
for (j = 0; j < 4; j++) {
unsigned char qb = blk[16 + 32 * j + lane];
for (h = 0; h < 2; h++) {
float sc, mn;
int nib = h ? (qb >> 4) : (qb & 0xF);
kq_scale_min(blk + 4, 2 * j + h, &sc, &mn);
acc = fmaf(xb[64 * j + 32 * h + lane],
(d * sc) * (float)nib - dmin * mn, acc);
}
}
p[lane] += acc;
}
}
}
static float ref_row_q4k(const unsigned char *row, const float *x, int cols)
{
float p[32] = { 0 };
ref_lanes_q4k(row, x, cols, p);
return tree32(p);
}
/* Replica of the kernel's sub-block activation sum: four threads own each
* 32-element group, thread `quarter` summing the 8 elements
* 8*quarter + (k ^ (group & 7)) in k order (the XOR is the kernel's LDS
* bank-conflict swizzle), then __shfl_down(2) + __shfl_down(1) folds the
* quarters as (q0 + q2) + (q1 + q3). */
static float ref_xsum_group(const float *xg, int group)
{
float q[4];
int quarter, k;
for (quarter = 0; quarter < 4; quarter++) {
float s = 0.0f;
for (k = 0; k < 8; k++)
s += xg[8 * quarter + (k ^ (group & 7))];
q[quarter] = s;
}
return (q[0] + q[2]) + (q[1] + q[3]);
}
/* Replica of the MIN-LIFTED fc_matvec_q4k (matvec_kquants.cu), in the
* kernel's exact op order:
* - the SCALE half accumulates sc*nib*x per lane and applies d ONCE per
* super-block (fmaf(d, acc, p));
* - the MIN half accumulates dmin*mn*xsum at one (super-block, sub-block)
* pair per lane per pass, where xsum is the block's __shfl_down tree sum
* of the 32 staged activations of that sub-block;
* - the two are subtracted per lane, then ONE tree32.
* The 2048-element staging chunk is load-bearing: it fixes both the xsum
* grouping and the pair -> lane assignment, so it is replicated here. */
static float ref_row_q4k_minlift(const unsigned char *row, const float *x,
int cols)
{
float p[32] = { 0 }, m[32] = { 0 }, f[32];
int cbase, lane, sb, j, h, g, pi;
for (cbase = 0; cbase < cols; cbase += 2048) {
int cend = cols - cbase < 2048 ? cols - cbase : 2048;
int nsb = cend >> 8;
float xsum[64];
for (g = 0; g < (cend >> 5); g++)
xsum[g] = ref_xsum_group(x + cbase + (g << 5), g);
for (sb = 0; sb < nsb; sb++) {
const unsigned char *blk = row + (size_t)((cbase >> 8) + sb) * 144;
float d = half2f((unsigned short)(blk[0] | (blk[1] << 8)));
const float *xb = x + cbase + (sb << 8);
for (lane = 0; lane < 32; lane++) {
float a = 0.0f;
for (j = 0; j < 4; j++) {
unsigned char qb = blk[16 + 32 * j + lane];
for (h = 0; h < 2; h++) {
float sc, mn;
int nib = h ? (qb >> 4) : (qb & 0xF);
kq_scale_min(blk + 4, 2 * j + h, &sc, &mn);
a = fmaf(sc * (float)nib,
xb[64 * j + 32 * h + lane], a);
}
}
p[lane] = fmaf(d, a, p[lane]);
}
}
for (lane = 0; lane < 32; lane++) {
for (pi = lane; pi < (nsb << 3); pi += 32) {
const unsigned char *blk = row +
(size_t)((cbase >> 8) + (pi >> 3)) * 144;
float dmin = half2f((unsigned short)(blk[2] | (blk[3] << 8)));
float sc, mn;
kq_scale_min(blk + 4, pi & 7, &sc, &mn);
m[lane] = fmaf(dmin * mn, xsum[pi], m[lane]);
}
}
}
for (lane = 0; lane < 32; lane++)
f[lane] = p[lane] - m[lane];
return tree32(f);
}
/* fp64 reference for the same dot — the yardstick both fp32 associations are
* measured against, so "the reassociation is numerically fine" is a number
* and not an assertion. Association-independent to within fp64 rounding. */
static double ref_row_q4k_f64(const unsigned char *row, const float *x,
int cols)
{
double acc = 0.0;
int n_sb = cols >> 8, sb, lane, j, h;
for (sb = 0; sb < n_sb; sb++) {
const unsigned char *blk = row + (size_t)sb * 144;
double d = half2f((unsigned short)(blk[0] | (blk[1] << 8)));
double dmin = half2f((unsigned short)(blk[2] | (blk[3] << 8)));
const float *xb = x + (sb << 8);
for (lane = 0; lane < 32; lane++)
for (j = 0; j < 4; j++) {
unsigned char qb = blk[16 + 32 * j + lane];
for (h = 0; h < 2; h++) {
float sc, mn;
int nib = h ? (qb >> 4) : (qb & 0xF);
kq_scale_min(blk + 4, 2 * j + h, &sc, &mn);
acc += ((d * sc) * nib - dmin * mn) *
(double)xb[64 * j + 32 * h + lane];
}
}
}
return acc;
}
static float ref_row_q5k(const unsigned char *row, const float *x, int cols)
{
float p[32] = { 0 };
int n_sb = cols >> 8, sb, lane, j, h;
for (sb = 0; sb < n_sb; sb++) {
const unsigned char *blk = row + (size_t)sb * 176;
float d = half2f((unsigned short)(blk[0] | (blk[1] << 8)));
float dmin = half2f((unsigned short)(blk[2] | (blk[3] << 8)));
const float *xb = x + (sb << 8);
for (lane = 0; lane < 32; lane++) {
unsigned char qh = blk[16 + lane];
float acc = 0.0f;
for (j = 0; j < 4; j++) {
unsigned char qb = blk[48 + 32 * j + lane];
for (h = 0; h < 2; h++) {
float sc, mn;
int nib = h ? (qb >> 4) : (qb & 0xF);
int hi = (qh & (1u << (2 * j + h))) ? 16 : 0;
kq_scale_min(blk + 4, 2 * j + h, &sc, &mn);
acc = fmaf(xb[64 * j + 32 * h + lane],
(d * sc) * (float)(nib + hi) - dmin * mn, acc);
}
}
p[lane] += acc;
}
}
return tree32(p);
}
static float ref_row_q6k(const unsigned char *row, const float *x, int cols)
{
float p[32] = { 0 };
int n_sb = cols >> 8, sb, lane, half, quad;
for (sb = 0; sb < n_sb; sb++) {
const unsigned char *blk = row + (size_t)sb * 210;
float d = half2f((unsigned short)(blk[208] | (blk[209] << 8)));
const float *xb = x + (sb << 8);
for (lane = 0; lane < 32; lane++) {
int lhalf = lane >> 4;
float acc = 0.0f;
for (half = 0; half < 2; half++) {
unsigned char qh = blk[128 + 32 * half + lane];
unsigned char ql0 = blk[64 * half + lane];
unsigned char ql1 = blk[64 * half + 32 + lane];
for (quad = 0; quad < 4; quad++) {
unsigned char qlb = (quad & 1) ? ql1 : ql0;
int nib = (quad & 2) ? (qlb >> 4) : (qlb & 0xF);
int hb = (qh >> (2 * quad)) & 3;
int q = (nib | (hb << 4)) - 32;
signed char sc = (signed char)blk[192 + half * 8 +
quad * 2 + lhalf];
acc = fmaf(xb[128 * half + 32 * quad + lane],
(float)(sc * q), acc);
}
}
p[lane] = fmaf(acc, d, p[lane]);
}
}
return tree32(p);
}
static float ref_row_q80(const unsigned char *row, const float *x, int cols)
{
float p[32] = { 0 };
int n_b = cols >> 5, b, lane;
for (b = 0; b < n_b; b++) {
const unsigned char *blk = row + (size_t)b * 34;
float d = half2f((unsigned short)(blk[0] | (blk[1] << 8)));
for (lane = 0; lane < 32; lane++)
p[lane] = fmaf(x[(b << 5) + lane],
d * (float)((signed char)blk[2 + lane]), p[lane]);
}
return tree32(p);
}
typedef float (*ref_row_fn)(const unsigned char *, const float *, int);
typedef void (*mv_fn)(fc_buffer *, fc_buffer *, fc_buffer *, int, int, int);
/* Only fc_matvec_q3k has a default-ctx wrapper; shim the rest. */
#define MV_SHIM(name) \
static void name##_dflt(fc_buffer *w, fc_buffer *x, fc_buffer *y, \
int r, int c, int rs) \
{ name##_ctx(fc_default_ctx(), w, x, y, r, c, rs); }
MV_SHIM(fc_matvec_q4k)
MV_SHIM(fc_matvec_q5k)
MV_SHIM(fc_matvec_q6k)
MV_SHIM(fc_matvec_q80)
MV_SHIM(fc_matvec_f32)
static int ulp_diff(float a, float b)
{
int32_t ia, ib;
int64_t d;
memcpy(&ia, &a, 4);
memcpy(&ib, &b, 4);
d = (int64_t)ia - ib;
if (d < 0) d = -d;
return (d > 1000000) ? 1000000 : (int)d;
}
/* Verify one quant matvec kernel vs the CPU replica, bit-exact. */
static int verify_mv(const char *name, int t, mv_fn fn, ref_row_fn ref,
int rows, int cols)
{
size_t raw_stride = raw_block_bytes(t) * (size_t)((t == 8) ? cols / 32
: cols / 256);
unsigned char *raw = rand_bytes(raw_stride * (size_t)rows);
float *xh = rand_floats((size_t)cols);
float *yh = malloc((size_t)rows * sizeof(float));
int dev_stride = 0, r, nmis = 0, maxulp = 0, ret;
fc_buffer *w, *x, *y;
sanitize_quant(t, raw, raw_stride * (size_t)rows);
w = upload_weight(t, raw, rows, cols, &dev_stride);
x = fc_buffer_from_ptr(xh, (size_t)cols * sizeof(float));
y = fc_buffer_alloc((size_t)rows * sizeof(float));
fn(w, x, y, rows, cols, dev_stride);
fc_sync();
fc_buffer_to_cpu(y, yh);
for (r = 0; r < rows; r++) {
float want = ref(raw + (size_t)r * raw_stride, xh, cols);
if (memcmp(&want, &yh[r], 4) != 0) {
int u = ulp_diff(want, yh[r]);
if (u > maxulp) maxulp = u;
if (nmis < 3)
printf(" row %d: ref %.9g (%08x) gpu %.9g (%08x)\n",
r, (double)want, *(unsigned *)&want,
(double)yh[r], *(unsigned *)&yh[r]);
nmis++;
}
}
printf("verify %-6s %6dx%-5d %s (mismatches %d, max ulp %d)\n",
name, rows, cols, nmis ? "FAIL" : "bit-exact", nmis, maxulp);
ret = nmis ? -1 : 0;
fc_buffer_free(w); fc_buffer_free(x); fc_buffer_free(y);
free(raw); free(xh); free(yh);
return ret;
}
/* Verify a REASSOCIATED kernel — one whose op order was deliberately changed
* for speed, so it cannot be bit-identical to the previous order.
*
* Two gates, both hard:
* 1. BIT-EXACT against a replica of the kernel's NEW op order. This is the
* real proof: it says the restructured algebra is exactly what the
* kernel computes, not merely "close". A dropped or misindexed
* correction term fails here immediately.
* 2. The fp32 relative error against an fp64 reference must stay under
* tol_rel AND within REASSOC_ERR_FACTOR of the OLD order's error. The
* second half is what catches a restructure that is algebraically right
* but numerically worse; an absolute cap alone would not.
* The old-vs-new ULP spread is reported for context (it is the amount real
* model outputs will move), not gated — it is expected to be nonzero. */
#define REASSOC_ERR_FACTOR 4.0
static int verify_mv_reassoc(const char *name, int t, mv_fn fn,
ref_row_fn ref_new, ref_row_fn ref_old,
double (*ref_f64)(const unsigned char *,
const float *, int),
int rows, int cols, double tol_rel)
{
size_t raw_stride = raw_block_bytes(t) * (size_t)(cols / 256);
unsigned char *raw = rand_bytes(raw_stride * (size_t)rows);
float *xh = rand_floats((size_t)cols);
float *yh = malloc((size_t)rows * sizeof(float));
int dev_stride = 0, r, nmis = 0, maxulp = 0, drift_ulp = 0, ret;
double rel_new = 0.0, rel_old = 0.0;
fc_buffer *w, *x, *y;
sanitize_quant(t, raw, raw_stride * (size_t)rows);
w = upload_weight(t, raw, rows, cols, &dev_stride);
x = fc_buffer_from_ptr(xh, (size_t)cols * sizeof(float));
y = fc_buffer_alloc((size_t)rows * sizeof(float));
fn(w, x, y, rows, cols, dev_stride);
fc_sync();
fc_buffer_to_cpu(y, yh);
for (r = 0; r < rows; r++) {
const unsigned char *rp = raw + (size_t)r * raw_stride;
float want = ref_new(rp, xh, cols);
float prev = ref_old(rp, xh, cols);
double exact = ref_f64(rp, xh, cols);
double scale = fabs(exact) > 1e-6 ? fabs(exact) : 1e-6;
double en = fabs((double)want - exact) / scale;
double eo = fabs((double)prev - exact) / scale;
int d = ulp_diff(prev, want);
if (en > rel_new) rel_new = en;
if (eo > rel_old) rel_old = eo;
if (d > drift_ulp) drift_ulp = d;
if (memcmp(&want, &yh[r], 4) != 0) {
int u = ulp_diff(want, yh[r]);
if (u > maxulp) maxulp = u;
if (nmis < 3)
printf(" row %d: ref %.9g (%08x) gpu %.9g (%08x)\n",
r, (double)want, *(unsigned *)&want,
(double)yh[r], *(unsigned *)&yh[r]);
nmis++;
}
}
ret = (nmis || rel_new > tol_rel ||
rel_new > REASSOC_ERR_FACTOR * rel_old) ? -1 : 0;
printf("verify %-6s %6dx%-5d %s (mismatches %d, max ulp %d)\n",
name, rows, cols, nmis ? "FAIL" : "bit-exact vs new order",
nmis, maxulp);
printf(" reassoc: max |new-old| %d ulp; rel err vs fp64 "
"new %.3g old %.3g (tol %.3g, factor cap %.1fx) %s\n",
drift_ulp, rel_new, rel_old, tol_rel, REASSOC_ERR_FACTOR,
ret ? "FAIL" : "ok");
fc_buffer_free(w); fc_buffer_free(x); fc_buffer_free(y);
free(raw); free(xh); free(yh);
return ret;
}
/* ------------------------------------------------------------------
* Grouped-MoE replicas (moe_expert_group.cu)
* ------------------------------------------------------------------ */
/* The down kernels' wpp launch arg — must mirror fc_moe_down_accum_ctx's
* env/default resolution so the replica reproduces the launched shape. */
static int moe_down_wpp(void)
{
const char *e = getenv("FC_MOE_DOWN_WPP");
if (e) {
int v = atoi(e);
if (v == 1 || v == 2 || v == 4 || v == 8)
return v;
}
return FC_MOE_DOWN_WPP_DEFAULT;
}
/* Replica of one fc_moe_gateup_q3k output: silu(gate.row . x) * (up.row . x).
* glibc expf and the device expf differ by a couple of ulp, so the caller
* compares with a small ulp tolerance instead of bit-exact. */
static float ref_moe_gateup(const unsigned char *gate_raw,
const unsigned char *up_raw, const float *x,
int cols)
{
float g = ref_row_q3k(gate_raw, x, cols);
float u = ref_row_q3k(up_raw, x, cols);
return (g / (1.0f + expf(-g))) * u;
}
/* Replica of one fc_moe_down_* output row in the kernel's exact op order:
* warp wk chains t = fmaf(top_wt[s], dot_s, t) per lane over slots wk,
* wk+wpp, ..., tree-reduces its lanes, and the wpp per-warp partials are
* added in ascending wk order. Bit-exact vs the kernel. hid should be the
* DEVICE-computed hid (downloaded) so this isolates the down kernel. */
static float ref_moe_down(int t, const unsigned char *down_raw,
size_t raw_row_stride, size_t raw_expert_stride,
const float *hid, const int32_t *idx,
const float *wt, int k, int row, int ff, int wpp)
{
float total = 0.0f;
int wk, s, l;
for (wk = 0; wk < wpp; wk++) {
float tl[32] = { 0 };
for (s = wk; s < k; s += wpp) {
const unsigned char *rp = down_raw +
(size_t)idx[s] * raw_expert_stride +
(size_t)row * raw_row_stride;
float p[32] = { 0 };
if (t == 11)
ref_lanes_q3k(rp, hid + (size_t)s * ff, ff, p);
else
ref_lanes_q4k(rp, hid + (size_t)s * ff, ff, p);
for (l = 0; l < 32; l++)
tl[l] = fmaf(wt[s], p[l], tl[l]);
}
total += tree32(tl);
}
return total;
}
/* ------------------------------------------------------------------
* Timing
* ------------------------------------------------------------------ */
static double bench_mv(const char *name, int t, mv_fn fn, int rows, int cols,
int iters)
{
size_t raw_stride = raw_block_bytes(t) * (size_t)((t == 8) ? cols / 32
: cols / 256);
unsigned char *raw = rand_bytes(raw_stride * (size_t)rows);
int dev_stride = 0, i;
fc_buffer *w, *x, *y;
double t0, dt, gb;
sanitize_quant(t, raw, raw_stride * (size_t)rows);
w = upload_weight(t, raw, rows, cols, &dev_stride);
free(raw);
{
float *xh = rand_floats((size_t)cols);
x = fc_buffer_from_ptr(xh, (size_t)cols * sizeof(float));
free(xh);
}
y = fc_buffer_alloc((size_t)rows * sizeof(float));
fn(w, x, y, rows, cols, dev_stride); /* warm up */
fc_sync();
t0 = now_s();
for (i = 0; i < iters; i++)
fn(w, x, y, rows, cols, dev_stride);
fc_sync();
dt = (now_s() - t0) / iters;
gb = (double)rows * dev_stride / dt / 1e9;
printf("%-24s %6dx%-5d %8.3f ms %7.2f GB/s\n", name, rows, cols,
dt * 1e3, gb);
fc_buffer_free(w); fc_buffer_free(x); fc_buffer_free(y);
return dt * 1e3;
}
/* Matched-shape per-quant sweep — the apples-to-apples question the per-type
* benches above cannot answer, because each of them uses a different shape.
* Same rows x cols for every type, so the ms column IS the cost of the quant
* CHOICE for a tensor of that shape. Device bytes per 256 weights after
* fc_repack_weight: q3k 152, q4k 152 (identical!), q5k 184, q6k 216, q80 288.
* Normalised to q4k so "what does staying on q6k cost me" reads directly. */
static void sweep_quants(int rows, int cols, int iters)
{
static const char *nm[5] = { "q3k", "q4k", "q5k", "q6k", "q80" };
static const int ty[5] = { 11, 12, 13, 14, 8 };
static const double bpb[5] = { 152, 152, 184, 216, 288 }; /* B / 256 wts */
mv_fn fn[5];
double ms[5];
int i;
fn[0] = fc_matvec_q3k; fn[1] = fc_matvec_q4k_dflt;
fn[2] = fc_matvec_q5k_dflt; fn[3] = fc_matvec_q6k_dflt;
fn[4] = fc_matvec_q80_dflt;
printf("-- matched shape %d x %d, %d iters --\n", rows, cols, iters);
for (i = 0; i < 5; i++)
ms[i] = bench_mv(nm[i], ty[i], fn[i], rows, cols, iters);
printf(" %-6s %9s %9s %9s %9s\n",
"type", "ms", "vs q4k", "dev B/wt", "ms/B ratio");
for (i = 0; i < 5; i++)
printf(" %-6s %9.3f %8.3fx %9.4f %9.3f\n", nm[i], ms[i],
ms[i] / ms[1], bpb[i] / 256.0,
(ms[i] / ms[1]) / (bpb[i] / bpb[1]));
printf(" (ms/B ratio > 1 = that quant's DECODE is more expensive per\n"
" byte than q4k's; = 1 = pure bandwidth/issue scaling)\n");
}
/* F32 matvec has no repack; bench it raw. */
static void bench_f32(int rows, int cols, int iters)
{
float *wh = rand_floats((size_t)rows * cols);
float *xh = rand_floats((size_t)cols);
fc_buffer *w = fc_buffer_from_ptr(wh, (size_t)rows * cols * sizeof(float));
fc_buffer *x = fc_buffer_from_ptr(xh, (size_t)cols * sizeof(float));
fc_buffer *y = fc_buffer_alloc((size_t)rows * sizeof(float));
int stride = cols * 4, i;
double t0, dt;
free(wh); free(xh);
fc_matvec_f32_dflt(w, x, y, rows, cols, stride);
fc_sync();
t0 = now_s();
for (i = 0; i < iters; i++)
fc_matvec_f32_dflt(w, x, y, rows, cols, stride);
fc_sync();
dt = (now_s() - t0) / iters;
printf("%-24s %6dx%-5d %8.3f ms %7.2f GB/s\n", "f32 router", rows, cols,
dt * 1e3, (double)rows * stride / dt / 1e9);
fc_buffer_free(w); fc_buffer_free(x); fc_buffer_free(y);
}
int main(void)
{
fc_cuda_ctx *ctx;
int it = 50, fails = 0;
if (fc_init() != 0) {
fprintf(stderr, "no CUDA\n");
return 1;
}
ctx = fc_default_ctx();
srand(42);
printf("== verify vs CPU replica (repack=%d) ==\n", USE_REPACK);
fails += verify_mv("q3k", 11, fc_matvec_q3k, ref_row_q3k, 512, 2048);
/* q4k is MIN-LIFTED and therefore reassociated — bit-exact against its
* own order's replica, plus an fp64 error gate. tol 1e-4: a 2048-term
* fp32 dot of random signed values cancels hard, so the ratio of the
* fp32 rounding error to the (small) result is orders of magnitude
* larger than eps; both orders sit far below this, and the factor cap
* is what actually holds the new order to the old one's accuracy. */
fails += verify_mv_reassoc("q4k", 12, fc_matvec_q4k_dflt,
ref_row_q4k_minlift, ref_row_q4k,
ref_row_q4k_f64, 512, 2048, 1e-4);
/* 512 cols: 2 super-blocks per chunk, so only lanes 0-15 run the min
* correction pass. 4096 cols: TWO staging chunks, so the xsum table is
* rebuilt and the pair -> lane assignment restarts mid-row. Both shapes
* exercise indexing the 2048 case cannot. */
fails += verify_mv_reassoc("q4k512", 12, fc_matvec_q4k_dflt,
ref_row_q4k_minlift, ref_row_q4k,
ref_row_q4k_f64, 2048, 512, 1e-4);
fails += verify_mv_reassoc("q4k4k", 12, fc_matvec_q4k_dflt,
ref_row_q4k_minlift, ref_row_q4k,
ref_row_q4k_f64, 512, 4096, 1e-4);
fails += verify_mv("q5k", 13, fc_matvec_q5k_dflt, ref_row_q5k, 512, 2048);
fails += verify_mv("q6k", 14, fc_matvec_q6k_dflt, ref_row_q6k, 512, 2048);
fails += verify_mv("q80", 8, fc_matvec_q80_dflt, ref_row_q80, 512, 2048);
fails += verify_mv("q3k512", 11, fc_matvec_q3k, ref_row_q3k, 2048, 512);
if (fails) {
printf("VERIFY FAILED — skipping timing\n");
fc_free();
return 1;
}
printf("== plain matvecs (rows x cols, device-weight GB/s) ==\n");
bench_mv("q3k gate/up (expert)", 11, fc_matvec_q3k, 512, 2048, it);
bench_mv("q3k down (expert)", 11, fc_matvec_q3k, 2048, 512, it);
bench_mv("q3k qkv-ish", 11, fc_matvec_q3k, 8192, 2048, it);
bench_mv("q4k 2048x512", 12, fc_matvec_q4k_dflt, 2048, 512, it);
bench_mv("q5k qkv 8192x2048", 13, fc_matvec_q5k_dflt, 8192, 2048, it);
bench_mv("q6k head 248320x2048", 14, fc_matvec_q6k_dflt, 248320, 2048, 10);
bench_mv("q80 4096x2048", 8, fc_matvec_q80_dflt, 4096, 2048, it);
bench_f32(256, 2048, it);
printf("== per-quant MATCHED-SHAPE sweep (what the quant choice costs) ==\n");
sweep_quants(2048, 2048, it); /* generic dense / ssm projection */
sweep_quants(8192, 2048, it); /* attn qkv */
sweep_quants(248320, 2048, 10); /* LM head (saber/ornith vocab) */
printf("== grouped MoE (per-layer shapes, k=8) ==\n");
{
int n_expert = 256, ff = 512, cols = 2048, k = 8;
int gu_stride = 0, dn3_stride = 0, dn4_stride = 0;
size_t graw = raw_block_bytes(11) * (cols / 256);
size_t d3raw = raw_block_bytes(11) * (ff / 256);
size_t d4raw = raw_block_bytes(12) * (ff / 256);
/* Only k experts are touched per launch; allocating all 256 experts
* matches the real memory layout / cache behavior. */
unsigned char *gr = rand_bytes((size_t)n_expert * ff * graw);
unsigned char *ur = rand_bytes((size_t)n_expert * ff * graw);
unsigned char *d3r = rand_bytes((size_t)n_expert * cols * d3raw);
unsigned char *d4r = rand_bytes((size_t)n_expert * cols * d4raw);
int32_t idx[8] = { 3, 17, 42, 77, 101, 155, 201, 250 };
float wt[8] = { .2f, .2f, .15f, .15f, .1f, .1f, .05f, .05f };
fc_buffer *gate, *up, *dn3, *dn4, *tidx, *twt, *x, *hid, *out;
float *xh = rand_floats((size_t)cols);
double t0, dt, gb;
int i;
sanitize_quant(11, gr, (size_t)n_expert * ff * graw);
sanitize_quant(11, ur, (size_t)n_expert * ff * graw);
sanitize_quant(11, d3r, (size_t)n_expert * cols * d3raw);
sanitize_quant(12, d4r, (size_t)n_expert * cols * d4raw);
gate = upload_weight(11, gr, (long)n_expert * ff, cols, &gu_stride);
up = upload_weight(11, ur, (long)n_expert * ff, cols, &gu_stride);
dn3 = upload_weight(11, d3r, (long)n_expert * cols, ff, &dn3_stride);
dn4 = upload_weight(12, d4r, (long)n_expert * cols, ff, &dn4_stride);
tidx = fc_buffer_from_ptr(idx, sizeof(idx));
twt = fc_buffer_from_ptr(wt, sizeof(wt));
x = fc_buffer_from_ptr(xh, (size_t)cols * sizeof(float));
hid = fc_buffer_alloc((size_t)k * ff * sizeof(float));
out = fc_buffer_alloc((size_t)cols * sizeof(float));
/* -- verify vs the CPU replicas before timing: gateup with a small
* ulp tolerance (expf), the cooperative down kernels bit-exact
* against the downloaded device hid -- */
{
float *hid_h = malloc((size_t)k * ff * sizeof(float));
float *out_h = malloc((size_t)cols * sizeof(float));
int wpp = moe_down_wpp();
int s, r, q, moe_fails = 0;
fc_moe_gateup_q3k_ctx(ctx, gate, up, x, tidx, hid, k, ff, cols,
gu_stride);
fc_sync();
fc_buffer_to_cpu(hid, hid_h);
{
int maxulp = 0, nmis = 0;
for (s = 0; s < k; s++)
for (r = 0; r < ff; r++) {
const unsigned char *grow =
gr + ((size_t)idx[s] * ff + r) * graw;
const unsigned char *urow =
ur + ((size_t)idx[s] * ff + r) * graw;
float want = ref_moe_gateup(grow, urow, xh, cols);
int u = ulp_diff(want, hid_h[(size_t)s * ff + r]);
if (u > maxulp) maxulp = u;
if (u > 8) nmis++;
}
printf("verify moe gateup q3k %s (mismatches %d, "
"max ulp %d, tol 8)\n",
nmis ? "FAIL" : "ok", nmis, maxulp);
if (nmis) moe_fails++;
}
for (q = 0; q < 2; q++) {
int t = q ? 12 : 11;
size_t rstride = q ? d4raw : d3raw;
const unsigned char *draw = q ? d4r : d3r;
int maxulp = 0, nmis = 0;
fc_moe_down_accum_ctx(ctx, t, q ? dn4 : dn3, hid, tidx, twt,
out, k, cols, ff,
q ? dn4_stride : dn3_stride);
fc_sync();
fc_buffer_to_cpu(out, out_h);
for (r = 0; r < cols; r++) {
float want = ref_moe_down(t, draw, rstride,
(size_t)cols * rstride, hid_h,
idx, wt, k, r, ff, wpp);
if (memcmp(&want, &out_h[r], 4) != 0) {
int u = ulp_diff(want, out_h[r]);
if (u > maxulp) maxulp = u;
if (nmis < 3)
printf(" row %d: ref %.9g (%08x) gpu %.9g "
"(%08x)\n", r, (double)want,
*(unsigned *)&want, (double)out_h[r],
*(unsigned *)&out_h[r]);
nmis++;
}
}
printf("verify moe down %s (wpp=%d) %s (mismatches %d, "
"max ulp %d)\n", q ? "q4k" : "q3k", wpp,
nmis ? "FAIL" : "bit-exact", nmis, maxulp);
if (nmis) moe_fails++;
}
free(hid_h); free(out_h);
if (moe_fails) {
printf("MOE VERIFY FAILED — skipping timing\n");
fc_free();
return 1;
}
}
free(gr); free(ur); free(d3r); free(d4r);
free(xh);
fc_moe_gateup_q3k_ctx(ctx, gate, up, x, tidx, hid, k, ff, cols,
gu_stride);
fc_sync();
t0 = now_s();
for (i = 0; i < it; i++)
fc_moe_gateup_q3k_ctx(ctx, gate, up, x, tidx, hid, k, ff, cols,
gu_stride);
fc_sync();
dt = (now_s() - t0) / it;
gb = 2.0 * k * ff * gu_stride / dt / 1e9;
printf("%-24s %6s %8.3f ms %7.2f GB/s\n", "moe gateup q3k",
"8x512x2048", dt * 1e3, gb);
fc_moe_down_accum_ctx(ctx, 11, dn3, hid, tidx, twt, out, k, cols, ff,
dn3_stride);
fc_sync();
t0 = now_s();
for (i = 0; i < it; i++)
fc_moe_down_accum_ctx(ctx, 11, dn3, hid, tidx, twt, out, k, cols,
ff, dn3_stride);
fc_sync();
dt = (now_s() - t0) / it;
gb = (double)k * cols * dn3_stride / dt / 1e9;
printf("%-24s %6s %8.3f ms %7.2f GB/s\n", "moe down q3k",
"8x2048x512", dt * 1e3, gb);
fc_moe_down_accum_ctx(ctx, 12, dn4, hid, tidx, twt, out, k, cols, ff,
dn4_stride);
fc_sync();
t0 = now_s();
for (i = 0; i < it; i++)
fc_moe_down_accum_ctx(ctx, 12, dn4, hid, tidx, twt, out, k, cols,
ff, dn4_stride);
fc_sync();
dt = (now_s() - t0) / it;
gb = (double)k * cols * dn4_stride / dt / 1e9;
printf("%-24s %6s %8.3f ms %7.2f GB/s\n", "moe down q4k",
"8x2048x512", dt * 1e3, gb);
fc_buffer_free(gate); fc_buffer_free(up);
fc_buffer_free(dn3); fc_buffer_free(dn4);
fc_buffer_free(tidx); fc_buffer_free(twt);
fc_buffer_free(x); fc_buffer_free(hid); fc_buffer_free(out);
}
fc_free();
return 0;
}
Reference in New Issue View Git Blame Copy Permalink