|
| 1 | +package kernels |
| 2 | + |
| 3 | +import ( |
| 4 | + "math" |
| 5 | + "testing" |
| 6 | +) |
| 7 | + |
| 8 | +// gemvReductionAbsTol and gemvReductionRelTol are the standing |
| 9 | +// numpy-allclose-style tolerance gate for the GEMV kernel family's fp32 |
| 10 | +// accumulation (sgemv_m1.cu, gemv_q4k.cu / gemv_q4k_sm121.cu): an element |
| 11 | +// passes when |
| 12 | +// |
| 13 | +// |got - want| <= gemvReductionAbsTol + gemvReductionRelTol*|want| |
| 14 | +// |
| 15 | +// Ported from zerfoo's fork of this kernel family (zerfoo#847, zerfoo PR |
| 16 | +// #934, T135.3); see zerfoo's docs/kernel-tolerances.md for the full per-op |
| 17 | +// tolerance table and rationale. |
| 18 | +// |
| 19 | +// Both kernels reduce K/N elements with a FIXED, deterministic order (each |
| 20 | +// warp lane sequentially accumulates a strided subset, then a 5-level |
| 21 | +// warp-shuffle tree combines the 32 lane partials) -- this is not |
| 22 | +// nondeterministic accumulation. It is, however, a DIFFERENT valid |
| 23 | +// parenthesization of the sum than the naive left-to-right CPU/float64 |
| 24 | +// reference (cpuSgemv / buildQ4KTestData's reference) used by these tests, so |
| 25 | +// fp32 rounding legitimately differs between the two orders. |
| 26 | +// |
| 27 | +// A pure RELATIVE bound is the wrong shape for this failure mode: the |
| 28 | +// synthetic sin()-based test data produces occasional near-zero row sums |
| 29 | +// (catastrophic cancellation), and for those rows the ABSOLUTE error stays a |
| 30 | +// few micro-units while the RELATIVE error explodes because the denominator |
| 31 | +// (the reference value) is itself tiny. Measured on the GB10 (2026-07-03, |
| 32 | +// T135.3, ref 08531b5f), the true (full-array, not first-failure) worst |
| 33 | +// cases were: |
| 34 | +// - TestSgemvM1_MultipleSizes/large_4096x4096: y[3862] rel err 7.32e-3, |
| 35 | +// but |diff| = 5e-6 against want=6.30e-4. |
| 36 | +// - TestSgemvM1_MultipleSizes/gemma3_1b_6144x1536: y[1791] rel err 3.96e-3, |
| 37 | +// |diff| = 6e-6 against want=-1.521e-3. |
| 38 | +// - TestGemvQ4KF32_MultipleSizes/medium_64x512: rel err 7.55e-4, |
| 39 | +// |diff| ~ 6e-7 against want=7.94e-4. |
| 40 | +// |
| 41 | +// In every case |diff| stayed at or below ~6e-6. gemvReductionAbsTol=1e-5 |
| 42 | +// covers all of them with margin while gemvReductionRelTol=1e-4 keeps the |
| 43 | +// original tight relative bound for normal-magnitude elements (at |want|~1, |
| 44 | +// the combined bound is ~1.1e-4, essentially unchanged from the original flat |
| 45 | +// 1e-4 test). A real kernel bug (wrong index, dropped term, a |
| 46 | +// fast-math-class blowup like the tanh overflow in ztensor#125) produces |
| 47 | +// absolute errors many orders of magnitude above 1e-5 and stays caught. |
| 48 | +const ( |
| 49 | + gemvReductionAbsTol = 1e-5 |
| 50 | + gemvReductionRelTol = 1e-4 |
| 51 | +) |
| 52 | + |
| 53 | +// checkGemvRelError scans the FULL output array against the reference using |
| 54 | +// the combined absolute+relative bound (see gemvReductionAbsTol / |
| 55 | +// gemvReductionRelTol above), reports the true maximum relative error found, |
| 56 | +// and asserts once at the end. |
| 57 | +// |
| 58 | +// The original per-element loop called t.Errorf + `break` at the first |
| 59 | +// offending index, which meant the logged "max relative error" only ever |
| 60 | +// reflected the error at (or before) that first-broken element -- NOT the |
| 61 | +// true dataset-wide max. That under-reporting bug surfaced directly during |
| 62 | +// T135.3 tolerance tuning: changing the bound changed WHICH element the loop |
| 63 | +// broke on, so the logged max jumped between runs in a way that looked like |
| 64 | +// kernel nondeterminism but was actually just the test giving up early at a |
| 65 | +// different point each time (the kernel itself is bit-reproducible). Scan to |
| 66 | +// completion so the reported max is honest regardless of where the |
| 67 | +// tolerance line sits. |
| 68 | +func checkGemvRelError(t *testing.T, got, ref []float32, absTol, relTol float64) { |
| 69 | + t.Helper() |
| 70 | + |
| 71 | + maxRelErr := 0.0 |
| 72 | + maxIdx := -1 |
| 73 | + badCount := 0 |
| 74 | + const maxReported = 5 |
| 75 | + for i := range got { |
| 76 | + absRef := math.Abs(float64(ref[i])) |
| 77 | + diff := math.Abs(float64(got[i] - ref[i])) |
| 78 | + |
| 79 | + var relErr float64 |
| 80 | + if absRef > 1e-6 { |
| 81 | + relErr = diff / absRef |
| 82 | + } else { |
| 83 | + relErr = diff |
| 84 | + } |
| 85 | + if relErr > maxRelErr { |
| 86 | + maxRelErr = relErr |
| 87 | + maxIdx = i |
| 88 | + } |
| 89 | + |
| 90 | + if diff > absTol+relTol*absRef { |
| 91 | + badCount++ |
| 92 | + if badCount <= maxReported { |
| 93 | + t.Errorf("y[%d] = %f, want %f (diff %e > %e + %e*%e)", |
| 94 | + i, got[i], ref[i], diff, absTol, relTol, absRef) |
| 95 | + } |
| 96 | + } |
| 97 | + } |
| 98 | + if badCount > maxReported { |
| 99 | + t.Errorf("... and %d more elements exceeded tol", badCount-maxReported) |
| 100 | + } |
| 101 | + t.Logf("max relative error: %e (at index %d)", maxRelErr, maxIdx) |
| 102 | +} |
0 commit comments