diff --git a/internal/cuda/kernels/gemv_q4k_sm121_test.go b/internal/cuda/kernels/gemv_q4k_sm121_test.go index 1bb64b1..30016c0 100644 --- a/internal/cuda/kernels/gemv_q4k_sm121_test.go +++ b/internal/cuda/kernels/gemv_q4k_sm121_test.go @@ -277,27 +277,8 @@ func TestQ4KGEMVOptimized(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e (sm_121=%v)", maxRelErr, IsQ4KSm121Supported()) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) + t.Logf("sm_121=%v", IsQ4KSm121Supported()) }) } } diff --git a/internal/cuda/kernels/gemv_q4k_test.go b/internal/cuda/kernels/gemv_q4k_test.go index 1fd4835..d64feca 100644 --- a/internal/cuda/kernels/gemv_q4k_test.go +++ b/internal/cuda/kernels/gemv_q4k_test.go @@ -253,27 +253,7 @@ func TestGemvQ4KF32_Parity(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e", maxRelErr) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) } func TestGemvQ4KF32_LargerMatrix(t *testing.T) { @@ -328,27 +308,7 @@ func TestGemvQ4KF32_LargerMatrix(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e", maxRelErr) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) } func TestGemvQ4KF32_MultipleSizes(t *testing.T) { @@ -416,27 +376,7 @@ func TestGemvQ4KF32_MultipleSizes(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e", maxRelErr) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) }) } } diff --git a/internal/cuda/kernels/sgemv_m1.cu b/internal/cuda/kernels/sgemv_m1.cu index 5fcb829..ae540b6 100644 --- a/internal/cuda/kernels/sgemv_m1.cu +++ b/internal/cuda/kernels/sgemv_m1.cu @@ -11,6 +11,7 @@ */ #include +#include #define BLOCK_SIZE 256 #define WARP_SIZE 32 @@ -40,21 +41,35 @@ __global__ void sgemv_m1_kernel( const float* row_ptr = A + (size_t)row * N; float acc = 0.0f; - /* Vectorized float4 path for the bulk of the row. */ - int n4 = N / 4; - const float4* row4 = (const float4*)row_ptr; - const float4* sx4 = (const float4*)sx; - - for (int i = lane_id; i < n4; i += WARP_SIZE) { - float4 a4 = __ldg(&row4[i]); - float4 x4 = sx4[i]; - acc = __fmaf_rn(a4.x, x4.x, acc); - acc = __fmaf_rn(a4.y, x4.y, acc); - acc = __fmaf_rn(a4.z, x4.z, acc); - acc = __fmaf_rn(a4.w, x4.w, acc); + /* Vectorized float4 path for the bulk of the row -- ONLY when row_ptr is + * 16-byte aligned. row_ptr = A + row*N; A itself is at least 16-byte + * aligned (cudaMalloc), but when N is not a multiple of 4 the per-row + * byte offset (row*N*4) is not a multiple of 16 for most rows, so + * reinterpreting row_ptr as float4* and __ldg-loading it is a misaligned + * vector load -- a hard fault ("misaligned address"), not just slow, and + * it poisons the whole CUDA context for the rest of the process (#847 + * tail, T135.3). Gate the vectorized path on actual pointer alignment + * instead of assuming N % 4 == 0; misaligned rows fall back to the + * scalar loop below, which is always correct regardless of N or row. */ + bool row_aligned = ((reinterpret_cast(row_ptr) & 0xF) == 0); + + int n4 = row_aligned ? (N / 4) : 0; + if (row_aligned) { + const float4* row4 = (const float4*)row_ptr; + const float4* sx4 = (const float4*)sx; + + for (int i = lane_id; i < n4; i += WARP_SIZE) { + float4 a4 = __ldg(&row4[i]); + float4 x4 = sx4[i]; + acc = __fmaf_rn(a4.x, x4.x, acc); + acc = __fmaf_rn(a4.y, x4.y, acc); + acc = __fmaf_rn(a4.z, x4.z, acc); + acc = __fmaf_rn(a4.w, x4.w, acc); + } } - /* Handle remainder elements (N not divisible by 4). */ + /* Scalar remainder: the tail when row_aligned (N not divisible by 4), or + * the entire row when !row_aligned. */ int rem_start = n4 * 4; for (int i = rem_start + lane_id; i < N; i += WARP_SIZE) { acc = __fmaf_rn(row_ptr[i], sx[i], acc); diff --git a/internal/cuda/kernels/sgemv_m1_test.go b/internal/cuda/kernels/sgemv_m1_test.go index 047d16e..22e4fdf 100644 --- a/internal/cuda/kernels/sgemv_m1_test.go +++ b/internal/cuda/kernels/sgemv_m1_test.go @@ -89,27 +89,7 @@ func TestSgemvM1_Parity(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e", maxRelErr) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) } func TestSgemvM1_MultipleSizes(t *testing.T) { @@ -178,27 +158,7 @@ func TestSgemvM1_MultipleSizes(t *testing.T) { t.Fatalf("Memcpy y: %v", err) } - maxRelErr := 0.0 - for i := range got { - absRef := math.Abs(float64(ref[i])) - diff := math.Abs(float64(got[i] - ref[i])) - var relErr float64 - if absRef > 1e-6 { - relErr = diff / absRef - } else { - relErr = diff - } - if relErr > maxRelErr { - maxRelErr = relErr - } - if relErr > 1e-4 { - t.Errorf("y[%d] = %f, want %f (rel err %e)", i, got[i], ref[i], relErr) - if t.Failed() { - break - } - } - } - t.Logf("max relative error: %e", maxRelErr) + checkGemvRelError(t, got, ref, gemvReductionAbsTol, gemvReductionRelTol) }) } } diff --git a/internal/cuda/kernels/tolerance_test.go b/internal/cuda/kernels/tolerance_test.go new file mode 100644 index 0000000..86a1f92 --- /dev/null +++ b/internal/cuda/kernels/tolerance_test.go @@ -0,0 +1,102 @@ +package kernels + +import ( + "math" + "testing" +) + +// gemvReductionAbsTol and gemvReductionRelTol are the standing +// numpy-allclose-style tolerance gate for the GEMV kernel family's fp32 +// accumulation (sgemv_m1.cu, gemv_q4k.cu / gemv_q4k_sm121.cu): an element +// passes when +// +// |got - want| <= gemvReductionAbsTol + gemvReductionRelTol*|want| +// +// Ported from zerfoo's fork of this kernel family (zerfoo#847, zerfoo PR +// #934, T135.3); see zerfoo's docs/kernel-tolerances.md for the full per-op +// tolerance table and rationale. +// +// Both kernels reduce K/N elements with a FIXED, deterministic order (each +// warp lane sequentially accumulates a strided subset, then a 5-level +// warp-shuffle tree combines the 32 lane partials) -- this is not +// nondeterministic accumulation. It is, however, a DIFFERENT valid +// parenthesization of the sum than the naive left-to-right CPU/float64 +// reference (cpuSgemv / buildQ4KTestData's reference) used by these tests, so +// fp32 rounding legitimately differs between the two orders. +// +// A pure RELATIVE bound is the wrong shape for this failure mode: the +// synthetic sin()-based test data produces occasional near-zero row sums +// (catastrophic cancellation), and for those rows the ABSOLUTE error stays a +// few micro-units while the RELATIVE error explodes because the denominator +// (the reference value) is itself tiny. Measured on the GB10 (2026-07-03, +// T135.3, ref 08531b5f), the true (full-array, not first-failure) worst +// cases were: +// - TestSgemvM1_MultipleSizes/large_4096x4096: y[3862] rel err 7.32e-3, +// but |diff| = 5e-6 against want=6.30e-4. +// - TestSgemvM1_MultipleSizes/gemma3_1b_6144x1536: y[1791] rel err 3.96e-3, +// |diff| = 6e-6 against want=-1.521e-3. +// - TestGemvQ4KF32_MultipleSizes/medium_64x512: rel err 7.55e-4, +// |diff| ~ 6e-7 against want=7.94e-4. +// +// In every case |diff| stayed at or below ~6e-6. gemvReductionAbsTol=1e-5 +// covers all of them with margin while gemvReductionRelTol=1e-4 keeps the +// original tight relative bound for normal-magnitude elements (at |want|~1, +// the combined bound is ~1.1e-4, essentially unchanged from the original flat +// 1e-4 test). A real kernel bug (wrong index, dropped term, a +// fast-math-class blowup like the tanh overflow in ztensor#125) produces +// absolute errors many orders of magnitude above 1e-5 and stays caught. +const ( + gemvReductionAbsTol = 1e-5 + gemvReductionRelTol = 1e-4 +) + +// checkGemvRelError scans the FULL output array against the reference using +// the combined absolute+relative bound (see gemvReductionAbsTol / +// gemvReductionRelTol above), reports the true maximum relative error found, +// and asserts once at the end. +// +// The original per-element loop called t.Errorf + `break` at the first +// offending index, which meant the logged "max relative error" only ever +// reflected the error at (or before) that first-broken element -- NOT the +// true dataset-wide max. That under-reporting bug surfaced directly during +// T135.3 tolerance tuning: changing the bound changed WHICH element the loop +// broke on, so the logged max jumped between runs in a way that looked like +// kernel nondeterminism but was actually just the test giving up early at a +// different point each time (the kernel itself is bit-reproducible). Scan to +// completion so the reported max is honest regardless of where the +// tolerance line sits. +func checkGemvRelError(t *testing.T, got, ref []float32, absTol, relTol float64) { + t.Helper() + + maxRelErr := 0.0 + maxIdx := -1 + badCount := 0 + const maxReported = 5 + for i := range got { + absRef := math.Abs(float64(ref[i])) + diff := math.Abs(float64(got[i] - ref[i])) + + var relErr float64 + if absRef > 1e-6 { + relErr = diff / absRef + } else { + relErr = diff + } + if relErr > maxRelErr { + maxRelErr = relErr + maxIdx = i + } + + if diff > absTol+relTol*absRef { + badCount++ + if badCount <= maxReported { + t.Errorf("y[%d] = %f, want %f (diff %e > %e + %e*%e)", + i, got[i], ref[i], diff, absTol, relTol, absRef) + } + } + } + if badCount > maxReported { + t.Errorf("... and %d more elements exceeded tol", badCount-maxReported) + } + t.Logf("max relative error: %e (at index %d)", maxRelErr, maxIdx) +}