From f82646005d01061cce72e558ce8989b840b23dd4 Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Wed, 4 Nov 2026 14:27:18 -0500 Subject: [PATCH 1/8] HIP: two-array PFPL stream, AMD kernels, and an FMA-safe ABS bound check Work in progress, on a branch for testing. - fz::pfpl (modules/fused/chunk_fusion/pfpl_stream.{h,cu}): compress a float field from one device buffer into another and back, in one kernel each, ordered on a stream with no host synchronization. The archive is byte-identical to the PFPL preset's in ABS mode. - AMD wave64 kernels for it (pfpl_stream_amd.inl), used on HIP. CUDA uses the generic chunk harness. - include/backend/fp.h: fmulRounded(). On HIP, clang fused the bound check's multiply into its subtraction, so the check tested the exact product and passed values that decode up to half an ulp past the bound. - backend/warp.h whole-wavefront primitives, chunk_fusion.cuh on the backend wrappers, and no NVRTC fusion on HIP. Tested on MI250X with the HIP backend on development only. main configures for CUDA alone, so this branch has not been built yet. Still missing for a PR: CHANGELOG entries, an in-tree test, a CN-* note, and removing two device-side printf calls. Co-Authored-By: Claude Opus 5.5 --- include/backend/api.h | 2 + include/backend/fp.h | 33 ++ include/backend/warp.h | 50 ++ modules/fused/chunk_fusion/CMakeLists.txt | 1 + modules/fused/chunk_fusion/chunk_fusion.cuh | 134 ++++-- modules/fused/chunk_fusion/pfpl_stream.cu | 362 ++++++++++++++ modules/fused/chunk_fusion/pfpl_stream.h | 101 ++++ .../fused/chunk_fusion/pfpl_stream_amd.inl | 450 ++++++++++++++++++ modules/quantizers/quantizer/quantizer.cu | 10 +- .../shufflers/bitshuffle/bitshuffle_stage.cu | 4 +- src/pipeline/compressor.cpp | 4 + src/pipeline/fusion_registry.cpp | 8 + 12 files changed, 1107 insertions(+), 52 deletions(-) create mode 100644 include/backend/fp.h create mode 100644 modules/fused/chunk_fusion/pfpl_stream.cu create mode 100644 modules/fused/chunk_fusion/pfpl_stream.h create mode 100644 modules/fused/chunk_fusion/pfpl_stream_amd.inl diff --git a/include/backend/api.h b/include/backend/api.h index 613cb7bd..c867bb10 100644 --- a/include/backend/api.h +++ b/include/backend/api.h @@ -46,6 +46,7 @@ #define cudaDevAttrMaxSharedMemoryPerBlock hipDeviceAttributeMaxSharedMemoryPerBlock #define cudaDevAttrMaxSharedMemoryPerBlockOptin hipDeviceAttributeSharedMemPerBlockOptin #define cudaDevAttrMultiProcessorCount hipDeviceAttributeMultiprocessorCount +#define cudaDevAttrCooperativeLaunch hipDeviceAttributeCooperativeLaunch // ── Handle types #define cudaError_t hipError_t @@ -75,6 +76,7 @@ #define cudaDeviceGetAttribute hipDeviceGetAttribute #define cudaDeviceSynchronize hipDeviceSynchronize #define cudaFuncSetAttribute hipFuncSetAttribute +#define cudaLaunchCooperativeKernel hipLaunchCooperativeKernel #define cudaFuncAttributeMaxDynamicSharedMemorySize hipFuncAttributeMaxDynamicSharedMemorySize #define cudaOccupancyMaxActiveBlocksPerMultiprocessor hipOccupancyMaxActiveBlocksPerMultiprocessor #define cudaProfilerStart hipProfilerStart diff --git a/include/backend/fp.h b/include/backend/fp.h new file mode 100644 index 00000000..96f96c8e --- /dev/null +++ b/include/backend/fp.h @@ -0,0 +1,33 @@ +#pragma once + +/** + * @file include/backend/fp.h + * @brief Floating-point operations whose rounding must not depend on the + * compiler's contraction policy. + * + * Dependency-free (no runtime headers): also compiled by NVRTC. + */ + +namespace fz { +namespace backend { + +/** + * a * b rounded to float on its own, never fused into the operation that + * consumes it. An error-bound check `|x - q*step| <= eb` must test the product + * the decoder will store, which is the rounded one. nvcc honors __fmul_rn; + * HIP's __fmul_rn is a plain multiply that clang's default + * -ffp-contract=fast-honor-pragmas fuses with a following subtraction into one + * FMA (gfx90a: v_fma_f32), which tests the exact product and passes values + * whose stored reconstruction misses the bound by up to half an ulp. + */ +__device__ __forceinline__ float fmulRounded(float a, float b) { +#if defined(FZGMOD_BACKEND_HIP) + #pragma clang fp contract(off) + return a * b; +#else + return __fmul_rn(a, b); +#endif +} + +} // namespace backend +} // namespace fz diff --git a/include/backend/warp.h b/include/backend/warp.h index 1f8b80e6..b875b33c 100644 --- a/include/backend/warp.h +++ b/include/backend/warp.h @@ -57,6 +57,37 @@ using warp_mask_t = unsigned; inline constexpr warp_mask_t kFullMask = 0xffffffffu; #endif +// ── Whole-wavefront primitives ─────────────────────────────────────────── +// For code that works per hardware wavefront rather than per 32-lane group: +// kWaveSize lanes (64 on CDNA, 32 on NVIDIA), lane i at bit i of a ballot. +// The host pass of a HIP compile emits no device code, so it only needs some +// value; the device pass takes the target's. +#if defined(FZGMOD_BACKEND_HIP) +# if defined(__GFX9__) || !defined(__HIP_DEVICE_COMPILE__) +inline constexpr int kWaveSize = 64; // GCN/CDNA wavefronts are always 64 lanes +# else +inline constexpr int kWaveSize = __AMDGCN_WAVEFRONT_SIZE__; // RDNA: wave32 unless -mwavefrontsize64 +# endif +using wave_mask_t = unsigned long long; +#else +inline constexpr int kWaveSize = 32; +using wave_mask_t = unsigned; +#endif + +/** Ballot over the whole wavefront (all lanes must take part). */ +__device__ inline wave_mask_t ballotWave(int pred) { + return __ballot_sync(kFullMask, pred); +} + +/** 1-based index of the lowest set bit of a wavefront mask, 0 if none. */ +__device__ inline int ffsWave(wave_mask_t m) { +#if defined(FZGMOD_BACKEND_HIP) + return __ffsll(static_cast(m)); +#else + return __ffs(static_cast(m)); +#endif +} + template __device__ inline T shflUp(T val, unsigned delta, int width) { return __shfl_up_sync(kFullMask, val, delta, width); @@ -77,6 +108,14 @@ __device__ inline T shfl(T val, int srcLane, int width) { return __shfl_sync(kFullMask, val, srcLane, width); } +/** Sum over the whole wavefront, returned to every lane (all lanes must take part). */ +template +__device__ inline T waveSum(T v) { + #pragma unroll + for (int off = kWaveSize / 2; off > 0; off >>= 1) v += shflXor(v, off, kWaveSize); + return v; +} + /** 32-lane-sub-group ballot, normalized to bits [0:31] regardless of which * hardware half of a 64-lane wavefront the calling lane sits in. */ __device__ inline uint32_t ballotSync32(int pred) { @@ -90,6 +129,17 @@ __device__ inline uint32_t ballotSync32(int pred) { #endif } +/** 32-lane-sub-group __activemask, normalized to bits [0:31] like ballotSync32: + * the active lanes of the caller's own 32-lane half. */ +__device__ inline uint32_t activeMask32() { +#if defined(FZGMOD_BACKEND_HIP) + const unsigned long long m = __activemask(); + return static_cast((__lane_id() & 32u) ? (m >> 32) : m); +#else + return __activemask(); +#endif +} + /** 32-lane-sub-group __any_sync — true if `pred` is true for any lane in the * caller's own 32-lane half of the wavefront. */ __device__ inline bool anySync32(int pred) { diff --git a/modules/fused/chunk_fusion/CMakeLists.txt b/modules/fused/chunk_fusion/CMakeLists.txt index f1120759..62ffe768 100644 --- a/modules/fused/chunk_fusion/CMakeLists.txt +++ b/modules/fused/chunk_fusion/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(fzgmod_modules PRIVATE ${CMAKE_CURRENT_LIST_DIR}/chunk_fusion.cu + ${CMAKE_CURRENT_LIST_DIR}/pfpl_stream.cu ${CMAKE_CURRENT_LIST_DIR}/nvrtc_chunk_fusion.cpp ) diff --git a/modules/fused/chunk_fusion/chunk_fusion.cuh b/modules/fused/chunk_fusion/chunk_fusion.cuh index 92c371b7..bdf3ee63 100644 --- a/modules/fused/chunk_fusion/chunk_fusion.cuh +++ b/modules/fused/chunk_fusion/chunk_fusion.cuh @@ -40,6 +40,8 @@ #include "coders/lc_common/lc_chunk_components.cuh" #include "coders/lc_common/lc_clog_components.cuh" // d_CLOG / d_HCLOG #include "transforms/negabinary/negabinary.h" +#include "backend/fp.h" +#include "backend/warp.h" #include "transforms/zigzag/zigzag.h" #include @@ -51,17 +53,17 @@ using byte = uint8_t; // 5-stage register butterfly (copy of bitshuffle_stage.cu butterfly32). __device__ __forceinline__ unsigned butterfly32(unsigned a, int sublane) { - unsigned q = __shfl_xor_sync(0xffffffffu, a, 16, 32); + unsigned q = fz::backend::shflXor(a, 16, 32); a = ((sublane&16)==0) ? __byte_perm(a,q,(3u<<12)|(2u<<8)|(7u<<4)|6u) : __byte_perm(a,q,(5u<<12)|(4u<<8)|(1u<<4)|0u); - q = __shfl_xor_sync(0xffffffffu, a, 8, 32); + q = fz::backend::shflXor(a, 8, 32); a = ((sublane&8)==0) ? __byte_perm(a,q,(3u<<12)|(7u<<8)|(1u<<4)|5u) : __byte_perm(a,q,(6u<<12)|(2u<<8)|(4u<<4)|0u); - q = __shfl_xor_sync(0xffffffffu, a, 4, 32); unsigned m=0x0F0F0F0Fu; + q = fz::backend::shflXor(a, 4, 32); unsigned m=0x0F0F0F0Fu; a = ((sublane&4)==0) ? ((a&~m)|((q>>4)&m)) : (((q<<4)&~m)|(a&m)); - q = __shfl_xor_sync(0xffffffffu, a, 2, 32); m=0x33333333u; + q = fz::backend::shflXor(a, 2, 32); m=0x33333333u; a = ((sublane&2)==0) ? ((a&~m)|((q>>2)&m)) : (((q<<2)&~m)|(a&m)); - q = __shfl_xor_sync(0xffffffffu, a, 1, 32); m=0x55555555u; + q = fz::backend::shflXor(a, 1, 32); m=0x55555555u; a = ((sublane&1)==0) ? ((a&~m)|((q>>1)&m)) : (((q<<1)&~m)|(a&m)); return a; } @@ -98,7 +100,7 @@ struct QuantInplaceZigzag { bool representable = q > -(int)p.radius && q < (int)p.radius && fabsf(x) < p.threshold; if (representable && p.verify_reconstruction) - representable = fabsf(x - __fmul_rn(static_cast(q), p.ebx2)) <= p.abs_eb; + representable = fabsf(x - fz::backend::fmulRounded(static_cast(q), p.ebx2)) <= p.abs_eb; if (representable) c = (uint32_t)((q << 1) ^ (q >> 31)); else @@ -136,13 +138,17 @@ struct QuantSplitOutlier { // atomicAdd per outlier. Cuts global-counter contention up to 32x. Only // outlier lanes reach here, so __activemask() is exactly the outlier set // (tail-safe: partial final chunks just have fewer active lanes). - const unsigned outmask = __activemask(); + const unsigned outmask = fz::backend::activeMask32(); // this 32-lane group's outlier lanes const int leader = __ffs(outmask) - 1; const unsigned rank = __popc(outmask & ((1u << lane) - 1)); uint32_t warp_base = 0; if (lane == leader) warp_base = atomicAdd(side.out_count, __popc(outmask)); +#if defined(FZGMOD_BACKEND_HIP) + warp_base = fz::backend::shfl(warp_base, leader, 32); // broadcast base slot (mask inert on HIP) +#else warp_base = __shfl_sync(outmask, warp_base, leader); // broadcast base slot +#endif const uint32_t slot = warp_base + rank; if (slot < side.max) { side.out_idxs[slot] = (uint32_t)(base + i); // global element index @@ -385,7 +391,7 @@ struct GolombRiceCoder { for (uint32_t k = 0; k <= kMaxK; ++k) { uint32_t c = local_cost[k]; #pragma unroll - for (int off = 16; off > 0; off >>= 1) c += __shfl_down_sync(0xffffffffu, c, off, 32); + for (int off = 16; off > 0; off >>= 1) c += fz::backend::shflDown(c, off, 32); if (lane == 0) s_warp_cost[warp * (kMaxK + 1) + k] = c; } __syncthreads(); @@ -394,7 +400,7 @@ struct GolombRiceCoder { for (uint32_t k = 0; k <= kMaxK; ++k) { uint32_t c = (lane < nwarps) ? s_warp_cost[lane * (kMaxK + 1) + k] : 0u; #pragma unroll - for (int off = 16; off > 0; off >>= 1) c += __shfl_down_sync(0xffffffffu, c, off, 32); + for (int off = 16; off > 0; off >>= 1) c += fz::backend::shflDown(c, off, 32); if (lane == 0) s_total_cost[k] = c; } if (lane == 0) { @@ -420,10 +426,10 @@ struct GolombRiceCoder { uint32_t warp_excl = thread_total; #pragma unroll for (int d = 1; d < 32; d <<= 1) { - const uint32_t up = __shfl_up_sync(0xffffffffu, warp_excl, d, 32); + const uint32_t up = fz::backend::shflUp(warp_excl, d, 32); if (lane >= d) warp_excl += up; } - const uint32_t warp_bit_total = __shfl_sync(0xffffffffu, warp_excl, 31, 32); + const uint32_t warp_bit_total = fz::backend::shfl(warp_excl, 31, 32); warp_excl -= thread_total; const uint32_t warp_byte_len = (warp_bit_total + 7u) >> 3; if (lane == 0) s_warp_bytelen[warp] = warp_byte_len; @@ -433,10 +439,10 @@ struct GolombRiceCoder { uint32_t excl = v; #pragma unroll for (int d = 1; d < 32; d <<= 1) { - const uint32_t up = __shfl_up_sync(0xffffffffu, excl, d, 32); + const uint32_t up = fz::backend::shflUp(excl, d, 32); if (lane >= d) excl += up; } - const uint32_t grand_total = __shfl_sync(0xffffffffu, excl, nwarps - 1, 32); + const uint32_t grand_total = fz::backend::shfl(excl, nwarps - 1, 32); excl -= v; if (lane < nwarps) s_warp_bytebase[lane] = excl; if (lane == 0) *s_total_bytes = grand_total; @@ -604,7 +610,11 @@ constexpr unsigned long long kLookbackInclusive = 2ull << 62; constexpr unsigned long long kLookbackValueMask = (1ull << 62) - 1ull; __device__ __forceinline__ unsigned long long lookbackLoad(const unsigned long long* p) { +#if defined(FZGMOD_BACKEND_HIP) + return __hip_atomic_load(p, __ATOMIC_RELAXED, __HIP_MEMORY_SCOPE_AGENT); +#else return *reinterpret_cast(p); +#endif } // Block-wide shared->global copy to an arbitrarily aligned destination (after @@ -654,12 +664,12 @@ chunk_emit_direct(uint32_t cid, const ChunkEncoded& e, byte* __restrict__ archiv if (j >= 0) { while (((v = lookbackLoad(&state[j])) >> 62) == 0ull) {} } - const unsigned incl = __ballot_sync(0xffffffffu, (v >> 62) == 2ull); + const unsigned incl = fz::backend::ballotSync32((v >> 62) == 2ull); const unsigned upto = incl ? (__ffs(incl) - 1) : 31u; // nearest inclusive lane unsigned long long part = (lane <= upto) ? (v & kLookbackValueMask) : 0ull; for (int off = 16; off > 0; off >>= 1) - part += __shfl_down_sync(0xffffffffu, part, off); - prefix += __shfl_sync(0xffffffffu, part, 0); + part += fz::backend::shflDown(part, off, 32); + prefix += fz::backend::shfl(part, 0, 32); if (incl) break; } if (lane == 0) { @@ -703,41 +713,48 @@ chunk_fused_kernel(const float* __restrict__ in, size_t n, } // ── Inverse chunk harness (initial RZE evidence-gated path). ───────────────── -// One CTA consumes one packed LC chunk, keeps every intermediate in the same -// shared-memory ping-pong buffers as compression, and writes reconstructed -// floats once. `offsets` are payload-relative exclusive offsets computed from -// the archive's flagged size table by the launcher. +// Block-wide global->shared copy of `len` bytes from an arbitrarily aligned +// `src` into the 16-byte aligned shared `dst` (the mirror of storeChunkBytes): +// aligned word loads realigned with __funnelshift_r. `src` must lie in a buffer +// whose start is 4-byte aligned; nothing past the word holding src[len-1] is read. +__device__ __forceinline__ void loadChunkBytes(byte* __restrict__ dst, + const byte* __restrict__ src, uint32_t len) { + const uint32_t mis = (uint32_t)((size_t)src & 3u); + const uint32_t words = (len + 3u) / 4u; + const uint32_t* __restrict__ sw = reinterpret_cast(src - mis); + uint32_t* __restrict__ dw = reinterpret_cast(dst); + if (mis == 0) { + for (uint32_t w = threadIdx.x; w < words; w += TPB) dw[w] = sw[w]; + } else { + const uint32_t need = mis + len; // bytes from the aligned base + const uint32_t shift = mis * 8u; + for (uint32_t w = threadIdx.x; w < words; w += TPB) { + const uint32_t hi = (4u * (w + 1u) < need) ? sw[w + 1u] : 0u; + dw[w] = __funnelshift_r(sw[w], hi, shift); + } + } +} + +// Decode one packed chunk into floats. `payload` is the chunk's bytes in the +// archive, `entry` its size-table word (bit 31 = stored raw). `coded_full`: the +// chunk was encoded whole (ChunkBytes, bitshuffled) -- a chunk padded with zeros +// past the field's end still was; `cnt` floats are reconstructed and written to +// out_chunk[0, cnt). Every intermediate stays in shared memory. template __device__ __forceinline__ void -chunk_inverse_pfpl_body(const byte* __restrict__ archive, - const uint32_t* __restrict__ entries, - const uint32_t* __restrict__ offsets, - size_t output_bytes, float ebx2, - bool inplace_outliers, uint32_t quant_radius, - float* __restrict__ out) { +chunk_inverse_pfpl_at(const byte* __restrict__ payload, uint32_t entry, bool coded_full, int cnt, + float ebx2, bool inplace_outliers, uint32_t quant_radius, + float* __restrict__ out_chunk) { constexpr int NELEM = Geom::NELEM; constexpr int NPP = Geom::NPP; __shared__ __align__(16) uint32_t sA[NELEM]; __shared__ __align__(16) uint32_t sB[NELEM]; __shared__ __align__(16) byte sTemp[TEMP_BYTES]; - const uint32_t cid = blockIdx.x; - const size_t base_bytes = static_cast(cid) * ChunkBytes; - if (base_bytes >= output_bytes) return; - const int out_bytes = static_cast(min(static_cast(ChunkBytes), - output_bytes - base_bytes)); - const int cnt = out_bytes / static_cast(sizeof(uint32_t)); - const bool full = out_bytes == ChunkBytes; - const uint32_t entry = entries[cid]; const uint32_t stored = entry & 0x7fffffffu; const bool raw = (entry & 0x80000000u) != 0; - const uint32_t nchunks = static_cast( - (output_bytes + ChunkBytes - 1) / ChunkBytes); - const size_t header = 8u + static_cast(nchunks) * sizeof(uint32_t); - const byte* payload = archive + header + offsets[cid]; - for (uint32_t i = threadIdx.x; i < stored; i += blockDim.x) - reinterpret_cast(sA)[i] = payload[i]; + loadChunkBytes(reinterpret_cast(sA), payload, stored); __syncthreads(); uint32_t* cur = sA; @@ -752,7 +769,7 @@ chunk_inverse_pfpl_body(const byte* __restrict__ archive, } // Bitshuffle is self-inverse but its read/write indexing is reversed. - if (full) { + if (coded_full) { const int lane = threadIdx.x & 31; for (int i = threadIdx.x; i < NELEM; i += blockDim.x) { const unsigned a = cur[i / 32 + lane * NPP]; @@ -776,7 +793,7 @@ chunk_inverse_pfpl_body(const byte* __restrict__ archive, const int lane = threadIdx.x & 31; const int warp = threadIdx.x >> 5; for (int delta = 1; delta < 32; delta <<= 1) { - const int32_t up = __shfl_up_sync(0xffffffffu, v, delta); + const int32_t up = fz::backend::shflUp(v, delta, 32); if (lane >= delta) v += up; } if (lane == 31) scan[warp] = v; @@ -784,7 +801,7 @@ chunk_inverse_pfpl_body(const byte* __restrict__ archive, if (warp == 0) { int32_t w = lane < 16 ? scan[lane] : 0; for (int delta = 1; delta < 32; delta <<= 1) { - const int32_t up = __shfl_up_sync(0xffffffffu, w, delta); + const int32_t up = fz::backend::shflUp(w, delta, 32); if (lane >= delta) w += up; } if (lane < 16) scan[lane] = w; @@ -797,21 +814,46 @@ chunk_inverse_pfpl_body(const byte* __restrict__ archive, __syncthreads(); } - const size_t base_elem = base_bytes / sizeof(uint32_t); for (int i = threadIdx.x; i < cnt; i += blockDim.x) { // Difference reconstructs the quantizer's *zigzag code* stream. Undo // that final map before scaling back to the reconstructed float. const uint32_t code = static_cast(qout[i]); if (inplace_outliers && (code >> 1) >= quant_radius) { - out[base_elem + static_cast(i)] = __uint_as_float(code); + out_chunk[i] = __uint_as_float(code); } else { const int32_t q = static_cast( (code >> 1) ^ (0u - (code & 1u))); - out[base_elem + static_cast(i)] = static_cast(q) * ebx2; + out_chunk[i] = static_cast(q) * ebx2; } } } +// One CTA consumes one packed LC chunk, keeps every intermediate in the same +// shared-memory ping-pong buffers as compression, and writes reconstructed +// floats once. `offsets` are payload-relative exclusive offsets computed from +// the archive's flagged size table by the launcher. +template +__device__ __forceinline__ void +chunk_inverse_pfpl_body(const byte* __restrict__ archive, + const uint32_t* __restrict__ entries, + const uint32_t* __restrict__ offsets, + size_t output_bytes, float ebx2, + bool inplace_outliers, uint32_t quant_radius, + float* __restrict__ out) { + const uint32_t cid = blockIdx.x; + const size_t base_bytes = static_cast(cid) * ChunkBytes; + if (base_bytes >= output_bytes) return; + const int out_bytes = static_cast(min(static_cast(ChunkBytes), + output_bytes - base_bytes)); + const uint32_t nchunks = static_cast( + (output_bytes + ChunkBytes - 1) / ChunkBytes); + const size_t header = 8u + static_cast(nchunks) * sizeof(uint32_t); + chunk_inverse_pfpl_at( + archive + header + offsets[cid], entries[cid], out_bytes == ChunkBytes, + out_bytes / static_cast(sizeof(uint32_t)), ebx2, inplace_outliers, quant_radius, + out + base_bytes / sizeof(uint32_t)); +} + template __global__ void __launch_bounds__(TPB) chunk_inverse_pfpl_kernel(const byte* __restrict__ archive, diff --git a/modules/fused/chunk_fusion/pfpl_stream.cu b/modules/fused/chunk_fusion/pfpl_stream.cu new file mode 100644 index 00000000..2e424447 --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream.cu @@ -0,0 +1,362 @@ +/** + * @file modules/fused/chunk_fusion/pfpl_stream.cu + * @brief Two-array PFPL stream (pfpl_stream.h): the chunk-cooperative harness + * (chunk_fusion.cuh) driven by one encode and one decode kernel that each + * place their chunk in the archive by decoupled look-back. + */ + +#include "fused/chunk_fusion/pfpl_stream.h" +#include "fused/chunk_fusion/chunk_fusion.cuh" +#include "backend/api.h" +#include "backend/warp.h" +#include "cuda_check.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace fz { +namespace pfpl { +namespace { + +using fused::chunk::byte; +using fused::chunk::TPB; +using fused::chunk::Geom; +using fused::chunk::ChunkEncoded; +using fused::chunk::ChunkSideCtx; +using fused::chunk::QuantInplaceZigzag; +using fused::chunk::QuantInplaceZigzagParams; +using fused::chunk::DiffNegabinary; +using fused::chunk::Bitshuffle32; +using fused::chunk::RZECoder; + +constexpr uint32_t kElems = Geom::NELEM; // floats per chunk +constexpr uint32_t kLcHeader = 8; // LC archive: [u32 coded bytes][u32 chunks] +constexpr uint32_t kStreamHeader = sizeof(StreamHeader); +constexpr int kWave = fz::backend::kWaveSize; + +// Look-back word: [63:62] kind (0 empty, 1 the chunk's own size, 2 its inclusive +// prefix) | [61:32] epoch of the call that wrote it | [31:0] bytes. A word left by +// another call reads as empty, so the state is never cleared between calls. +constexpr unsigned long long kKindMask = 3ull << 62; +constexpr unsigned long long kAggregate = 1ull << 62; +constexpr unsigned long long kInclusive = 2ull << 62; +constexpr uint32_t kEpochMask = 0x3fffffffu; + +__device__ __forceinline__ unsigned long long lbWord(unsigned long long kind, uint32_t epoch, uint32_t bytes) { + return kind | (static_cast(epoch) << 32) | bytes; +} + +__device__ __forceinline__ bool lbCurrent(unsigned long long w, uint32_t epoch) { + return (w & kKindMask) != 0ull && static_cast((w >> 32) & kEpochMask) == epoch; +} + +__device__ __forceinline__ void lbPause() { +#if defined(FZGMOD_BACKEND_HIP) + __builtin_amdgcn_s_sleep(1); +#elif defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 700 + __nanosleep(32); +#endif +} + +__device__ __forceinline__ void streamTrap() { +#if defined(FZGMOD_BACKEND_HIP) + __builtin_trap(); +#else + __trap(); +#endif +} + +// Exclusive payload offset of chunk `cid`, found by one whole wavefront (every +// lane gets it): kWave predecessors at a time (lane k looks at chunk cid-1-k), +// stopping at the nearest inclusive prefix. Encode (`table` null) waits for a +// predecessor that has not published -- chunks are ticketed in start order, so +// it is running. Decode reads an unpublished predecessor's size from the +// archive's size table instead, so it never waits and needs no ticket. +__device__ uint32_t lookbackWave(uint32_t cid, uint32_t own, unsigned long long* __restrict__ state, + uint32_t epoch, const uint32_t* __restrict__ table) { + const int lane = static_cast(threadIdx.x) % kWave; + if (lane == 0 && table == nullptr) + atomicExch(&state[cid], lbWord(cid == 0 ? kInclusive : kAggregate, epoch, own)); + unsigned long long prefix = 0; + for (long long window = static_cast(cid) - 1; window >= 0; window -= kWave) { + const long long j = window - lane; + bool incl = true; // before chunk 0: an inclusive prefix of 0 + unsigned long long bytes = 0; + if (j >= 0) { + unsigned long long w = fused::chunk::lookbackLoad(&state[j]); + if (table == nullptr) { + while (!lbCurrent(w, epoch)) { + lbPause(); + w = fused::chunk::lookbackLoad(&state[j]); + } + } else if (!lbCurrent(w, epoch) || (w & kKindMask) != kInclusive) { + w = lbWord(kAggregate, epoch, table[j] & 0x7fffffffu); + } + incl = (w & kKindMask) == kInclusive; + bytes = w & 0xffffffffull; + } + const fz::backend::wave_mask_t inclMask = fz::backend::ballotWave(incl); + const int upto = inclMask ? fz::backend::ffsWave(inclMask) - 1 : kWave - 1; + prefix += fz::backend::waveSum(lane <= upto ? bytes : 0ull); + if (inclMask) break; + } + if (lane == 0 && (cid != 0 || table != nullptr)) + atomicExch(&state[cid], lbWord(kInclusive, epoch, static_cast(prefix) + own)); + return static_cast(prefix); +} + +// lookbackWave() run by the block's first wavefront, its result shared with the block. +__device__ uint32_t lookbackOffset(uint32_t cid, uint32_t own, unsigned long long* __restrict__ state, + uint32_t epoch, const uint32_t* __restrict__ table) { + __shared__ uint32_t s_offset; + if (threadIdx.x < kWave) { + const uint32_t offset = lookbackWave(cid, own, state, epoch, table); + if (threadIdx.x == 0) s_offset = offset; + } + __syncthreads(); + return s_offset; +} + +// One block per chunk: code it in shared memory, find its offset, write its +// size-table entry and payload into the stream; the last chunk adds the headers +// and the length. Layout: [StreamHeader][LC header][u32 sizes x chunks][payloads]. +__global__ void __launch_bounds__(TPB) +encodeKernel(const float* __restrict__ in, uint32_t n, uint32_t coded_n, QuantInplaceZigzagParams qp, + byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, + unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, + unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { + __shared__ uint32_t s_cid; + __shared__ QuantInplaceZigzagParams s_qp; + if (threadIdx.x == 0) { + s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); + s_qp = qp; + } + __syncthreads(); + const uint32_t cid = s_cid; + + const ChunkEncoded e = + fused::chunk::chunk_encode_at, + DiffNegabinary, Bitshuffle32>( + cid, in, coded_n, reinterpret_cast(&s_qp), ChunkSideCtx{}, n); + const uint32_t offset = lookbackOffset(cid, e.nbytes, state, epoch, nullptr); + + const unsigned long long table = kStreamHeader + kLcHeader; + const unsigned long long data = table + 4ull * chunks; + if (threadIdx.x == 0 && table + 4ull * (cid + 1u) <= capacity) + reinterpret_cast(out + table)[cid] = e.size_word; + if (data + offset + e.nbytes <= capacity) + fused::chunk::storeChunkBytes(out + data + offset, e.src, e.nbytes); + + if (cid == chunks - 1u && threadIdx.x == 0) { + const unsigned long long archive = kLcHeader + 4ull * chunks + offset + e.nbytes; + const unsigned long long padded = (archive + 3ull) & ~3ull; + const unsigned long long total = kStreamHeader + padded; + uint32_t* w = reinterpret_cast(out); + if (total <= capacity) { + for (unsigned long long b = kStreamHeader + archive; b < total; ++b) out[b] = 0; + w[4] = coded_n * 4u; + w[5] = chunks; + } + // The stream header goes in whenever it fits, so a stream too long for the + // buffer still records the length it needed. + if (kStreamHeader <= capacity) { + w[0] = kStreamMagic; + w[1] = static_cast(padded); + w[2] = __float_as_uint(qp.abs_eb); + w[3] = n; + } + if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); + } +} + +// A field of no elements: the headers alone. +__global__ void emptyKernel(byte* __restrict__ out, unsigned long long capacity, float abs_eb, + size_t* __restrict__ d_stream_bytes) { + const unsigned long long total = kStreamHeader + kLcHeader; + uint32_t* w = reinterpret_cast(out); + if (kStreamHeader <= capacity) { + w[0] = kStreamMagic; + w[1] = kLcHeader; + w[2] = __float_as_uint(abs_eb); + w[3] = 0u; + } + if (total <= capacity) { + w[4] = 0u; + w[5] = 0u; + } + if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); +} + +// One block per chunk (gridDim.x == chunks): find the chunk's payload, decode it. +__global__ void __launch_bounds__(TPB) +decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t n, + unsigned long long* __restrict__ state, uint32_t epoch) { + const uint32_t* w = reinterpret_cast(stream); + const uint32_t chunks = gridDim.x; + if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || + w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { + if (blockIdx.x == 0 && threadIdx.x == 0) + printf("[fzgmod] pfpl::decompress: not a PFPL stream of %u floats " + "(magic %08x, n %u, chunks %u, archive %u B)\n", n, w[0], w[3], w[5], w[1]); + streamTrap(); + } + const uint32_t cid = blockIdx.x; + const uint32_t* table = w + (kStreamHeader + kLcHeader) / 4u; + const uint32_t entry = table[cid]; + const uint32_t offset = lookbackOffset(cid, entry & 0x7fffffffu, state, epoch, table); + const byte* payload = stream + kStreamHeader + kLcHeader + 4ull * chunks + offset; + const uint32_t base = cid * kElems; + fused::chunk::chunk_inverse_pfpl_at>( + payload, entry, /*coded_full=*/true, static_cast(min(kElems, n - base)), + 2.0f * __uint_as_float(w[2]), /*inplace_outliers=*/true, kQuantRadius, out + base); +} + +#if defined(FZGMOD_BACKEND_HIP) +#include "fused/chunk_fusion/pfpl_stream_amd.inl" +#endif + +// FZ_PFPL_HARNESS=1 runs the generic harness kernels instead of the AMD ones (same stream; for comparison). +bool useHarness() { +#if defined(FZGMOD_BACKEND_HIP) + static const bool harness = [] { + const char* e = std::getenv("FZ_PFPL_HARNESS"); + return e != nullptr && e[0] == '1'; + }(); + return harness; +#else + return true; +#endif +} + +uint32_t chunksOf(size_t n) { return static_cast((n + kElems - 1) / kElems); } + +// The archive's size fields are 32-bit: even an all-raw stream must stay below 4 GiB. +bool fitsArchive(size_t n) { + return kStreamHeader + kLcHeader + static_cast(chunksOf(n)) * (4u + kChunkBytes) + < (static_cast(1) << 32); +} + +void requireAligned(const void* p, const char* what) { + if (reinterpret_cast(p) % 4u != 0) + throw std::invalid_argument(std::string("fz::pfpl: ") + what + " must be 4-byte aligned"); +} + +} // namespace + +size_t maxStreamBytes(size_t n) { + const size_t chunks = chunksOf(n); + return kStreamHeader + kLcHeader + chunks * (4u + static_cast(kChunkBytes)); +} + +size_t streamBytes(const void* h_stream, size_t capacity) { + if (h_stream == nullptr || capacity < kStreamHeader) return 0; + StreamHeader h; + std::memcpy(&h, h_stream, sizeof(h)); + if (h.magic != kStreamMagic || (h.archive_bytes & 3u) != 0 || + h.archive_bytes < kLcHeader + 4ull * chunksOf(h.n)) + return 0; + const size_t total = kStreamHeader + static_cast(h.archive_bytes); + return total <= capacity ? total : 0; +} + +Plan::~Plan() { + // Teardown: nothing useful to do with an error here. + if (state_ != nullptr) (void)cudaFree(state_); + if (ticket_ != nullptr) (void)cudaFree(ticket_); +} + +void Plan::reserve(size_t n) { + const size_t chunks = std::max(chunksOf(n), 1); + if (ticket_ == nullptr) { + FZ_CUDA_CHECK(cudaMalloc(&ticket_, sizeof(unsigned long long))); + FZ_CUDA_CHECK(cudaMemset(ticket_, 0, sizeof(unsigned long long))); + tickets_ = 0; + } + if (chunks <= chunks_) return; + if (state_ != nullptr) { // calls still in flight may read the old state + FZ_CUDA_CHECK(cudaDeviceSynchronize()); + FZ_CUDA_CHECK(cudaFree(state_)); + state_ = nullptr; + } + FZ_CUDA_CHECK(cudaMalloc(&state_, chunks * sizeof(unsigned long long))); + FZ_CUDA_CHECK(cudaMemset(state_, 0, chunks * sizeof(unsigned long long))); + chunks_ = chunks; + epoch_ = 0; // a cleared state: every epoch is new again +} + +uint32_t Plan::nextEpoch(fz::stream_t stream) { + if (++epoch_ > kEpochMask) { // after 2^30 calls, start over on a cleared state + FZ_CUDA_CHECK(cudaMemsetAsync(state_, 0, chunks_ * sizeof(unsigned long long), stream)); + epoch_ = 1; + } + return epoch_; +} + +void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { + if (!(abs_eb > 0.0f) || !std::isfinite(abs_eb) || !std::isfinite(0.5f / abs_eb)) + throw std::invalid_argument("fz::pfpl::compress: the error bound must be positive and finite, " + "with a finite quantizer scale 0.5/abs_eb"); + if (!fitsArchive(n)) + throw std::invalid_argument("fz::pfpl::compress: a field of 4 GiB or more (in whole chunks) exceeds " + "the LC archive's 32-bit sizes"); + requireAligned(d_out, "the stream buffer"); + plan.reserve(n); + const uint32_t epoch = plan.nextEpoch(stream); + + // examples/presets/pfpl_specialization.toml in ABS mode, as QuantizerStage::getFusedOp() packs it. + QuantInplaceZigzagParams qp{}; + qp.ebx2_r = 0.5f / abs_eb; // verify_reconstruction's operation order + qp.ebx2 = 2.0f * abs_eb; + qp.abs_eb = abs_eb; + qp.radius = kQuantRadius; + qp.threshold = std::numeric_limits::infinity(); + qp.verify_reconstruction = 1; + + const uint32_t chunks = chunksOf(n); + if (chunks == 0) { + emptyKernel<<<1, 1, 0, stream>>>(static_cast(d_out), capacity, abs_eb, d_stream_bytes); + FZ_CUDA_CHECK(cudaGetLastError()); + return; + } +#if defined(FZGMOD_BACKEND_HIP) + if (!useHarness()) + amd::encodeKernel<<>>(d_in, static_cast(n), qp, + static_cast(d_out), capacity, chunks, + plan.state_, plan.ticket_, plan.tickets_, epoch, + d_stream_bytes); + else +#endif + encodeKernel<<>>(d_in, static_cast(n), chunks * kElems, qp, + static_cast(d_out), capacity, chunks, plan.state_, + plan.ticket_, plan.tickets_, epoch, d_stream_bytes); + FZ_CUDA_CHECK(cudaGetLastError()); + plan.tickets_ += chunks; +} + +void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream) { + if (n == 0) return; + if (!fitsArchive(n)) + throw std::invalid_argument("fz::pfpl::decompress: no stream holds a field of 4 GiB or more"); + requireAligned(d_stream, "the stream"); + requireAligned(d_out, "the output"); + plan.reserve(n); + const uint32_t epoch = plan.nextEpoch(stream); +#if defined(FZGMOD_BACKEND_HIP) + if (!useHarness()) + amd::decodeKernel<<>>(static_cast(d_stream), d_out, + static_cast(n), plan.state_, epoch); + else +#endif + decodeKernel<<>>(static_cast(d_stream), d_out, + static_cast(n), plan.state_, epoch); + FZ_CUDA_CHECK(cudaGetLastError()); +} + +} // namespace pfpl +} // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream.h b/modules/fused/chunk_fusion/pfpl_stream.h new file mode 100644 index 00000000..a8afe931 --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream.h @@ -0,0 +1,101 @@ +#pragma once + +/** + * @file modules/fused/chunk_fusion/pfpl_stream.h + * @brief Two-array PFPL: a float field compresses straight into one caller + * buffer and decompresses straight back, stream-ordered, with no host + * synchronization and no device temporaries. + * + * The stream is self-describing: a 16-byte StreamHeader, then the LC/RZE + * archive that the PFPL preset (examples/presets/pfpl_specialization.toml with + * error_bound_mode = "ABS": Quantizer in-place zigzag -> Difference negabinary + * -> Bitshuffle -> RZE over 16 KiB chunks) produces for the same input and + * bound, byte for byte. Encoding and decoding are one kernel each over the + * chunk-cooperative harness (chunk_fusion.cuh): one block per chunk keeps every + * intermediate in shared memory and finds the chunk's place in the archive by + * decoupled look-back, so neither direction touches device memory beyond the + * two arrays and a Plan's 8 bytes of look-back state per chunk. + * + * Archive bytes are a multiple of 4; a field whose length is not a whole number + * of chunks is coded padded with zeros (as the staged pipeline pads it) and + * decodes to exactly its own length. Both arrays must be 4-byte aligned. + */ + +#include "backend/types.h" + +#include +#include + +namespace fz { +namespace pfpl { + +/// Leading 16 bytes of every stream. +struct StreamHeader { + uint32_t magic; ///< kStreamMagic + uint32_t archive_bytes; ///< LC archive bytes after the header (a multiple of 4) + float abs_eb; ///< absolute error bound every value was coded to + uint32_t n; ///< float elements of the field +}; +static_assert(sizeof(StreamHeader) == 16, "the PFPL stream header is 16 bytes"); + +constexpr uint32_t kStreamMagic = 0x4C504650u; ///< "PFPL" in little-endian byte order +constexpr int kChunkBytes = 16384; ///< the preset's Difference/Bitshuffle/RZE chunk +constexpr uint32_t kQuantRadius = 4194304; ///< the preset's quant_radius + +/** Longest stream n floats can produce (every chunk stored raw). */ +size_t maxStreamBytes(size_t n); + +/** Length of the stream whose first bytes are at `h_stream` (host memory), or 0 + * if they are not a PFPL stream header or claim more than `capacity` bytes. */ +size_t streamBytes(const void* h_stream, size_t capacity = SIZE_MAX); + +/** + * Per-stream look-back state. Calls through one Plan must be ordered (one CUDA/HIP + * stream, or synchronized between streams): use one Plan per concurrently used + * stream. It grows to the largest field it has seen. + */ +class Plan { +public: + Plan() = default; + ~Plan(); + Plan(const Plan&) = delete; + Plan& operator=(const Plan&) = delete; + + /** Make room for fields of up to n floats. Allocating synchronizes the device. */ + void reserve(size_t n); + +private: + friend void compress(Plan&, const float*, size_t, float, void*, size_t, size_t*, fz::stream_t); + friend void decompress(Plan&, const void*, float*, size_t, fz::stream_t); + uint32_t nextEpoch(fz::stream_t stream); + + unsigned long long* state_ = nullptr; ///< one look-back word per chunk, tagged with its call's epoch + unsigned long long* ticket_ = nullptr; ///< encode chunk tickets handed out on the device + size_t chunks_ = 0; ///< chunks state_ holds + unsigned long long tickets_ = 0; ///< the ticket counter's value between calls + uint32_t epoch_ = 0; ///< last epoch used (0 = never: fresh state reads as empty) +}; + +/** + * Compress n floats at d_in to within abs_eb into d_out, which holds `capacity` + * bytes. Stream-ordered, no host synchronization. *d_stream_bytes (device + * memory; may be null) receives the stream's length -- greater than `capacity` + * when it does not fit. Nothing at or past d_out + capacity is ever written; a + * stream that does not fit leaves its 16-byte header (when that fits), which + * records the length it needed, so streamBytes() reports the overflow too. + * abs_eb must be positive and finite; the field must code to under 4 GiB (the + * LC archive's size fields are 32-bit). + */ +void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream); + +/** + * Decompress the compress() stream at d_stream, which must describe n floats, + * into d_out[0, n). Stream-ordered, no host synchronization. A stream that is + * not a PFPL stream of n floats stops the kernel (device trap) rather than + * decoding garbage. + */ +void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream); + +} // namespace pfpl +} // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream_amd.inl b/modules/fused/chunk_fusion/pfpl_stream_amd.inl new file mode 100644 index 00000000..8d6b1327 --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream_amd.inl @@ -0,0 +1,450 @@ +// Included by pfpl_stream.cu (inside fz::pfpl::): the AMD (wave64) encode/decode kernels. +// +// The generic harness kernels keep a chunk in two 16 KiB shared buffers plus 4 KiB of coder scratch (36 KiB: one +// 512-thread block per CU on MI200's 64 KiB LDS) and run d_RZE's bitmap recursion with block-wide scans (~20 +// barriers per chunk); profiled on MI250X they wait ~65% of the time at 8 waves per CU. These produce the same +// stream, byte for byte, from one padded 16.1 KiB buffer and 2 KiB of bitmap (3 blocks per CU), with 6-7 barriers: +// - quantize, difference and bitshuffle in registers (the difference's cross-wavefront neighbours through a +// 64-word exchange), scattering the bitshuffled planes into a buffer padded to 129 words per plane so the +// plane scatter hits 32 different banks (the unpadded layout puts all 32 lanes of a group in one bank); +// - zero elimination with one block scan, then the bitmap levels L1 -> L4 (d_RZE's repetition eliminations of +// 2048, 256 and 32 bytes) by wavefront 0 alone with wave scans and ballots, and the look-back; +// - the packed chunk compacted in place in shared memory and stored once. +// Decoding mirrors it: an element of a repetition-eliminated level is data[k - 1], k the set bitmap bits up to and +// including it, so each level decodes in parallel from popcounts; the chunk-local prefix sum runs as eight +// wave scans at once with one barrier, in the element order the inverse bitshuffle leaves in registers. +// Tuned for, and only used on, 64-lane wavefronts. + +namespace amd { + +constexpr int kTPB = 512; // 8 wavefronts; must equal the harness's TPB (storeChunkBytes) +constexpr int kPerT = kElems / kTPB; // 8 elements per thread, and 8 bitshuffled words +constexpr int kNPP = kElems / 32; // words per bitshuffle plane +constexpr int kPS = kNPP + 1; // padded plane stride (words) +constexpr int kBufWords = 32 * kPS; +constexpr int kL1 = kChunkBytes / 8; // zero-elimination bitmap bytes +constexpr uint32_t kAvail = kChunkBytes - 2; // d_RZE: room before its 2-byte trailer +static_assert(kTPB == TPB && kTPB == 8 * 64 && kPerT == 8, "one 512-thread block per 16 KiB chunk"); +static_assert(kBufWords * 4 >= kChunkBytes + 4, "the packed chunk is compacted in place in the plane buffer"); + +__device__ __forceinline__ uint32_t inclMask(int j) { return j == 31 ? 0xffffffffu : (2u << j) - 1u; } + +template +__device__ __forceinline__ T waveInclusiveScan(T v) { + const int lane = threadIdx.x & 63; + #pragma unroll + for (int d = 1; d < 64; d <<= 1) { + const T up = fz::backend::shflUp(v, d, 64); + if (lane >= d) v += up; + } + return v; +} + +// Exclusive prefix of v over the block in thread order, and the total. One barrier; `waves` (8 words) must not be +// rewritten before a later barrier. +__device__ __forceinline__ uint32_t blockExclusiveScan(uint32_t v, uint32_t* waves, uint32_t& total) { + const int lane = threadIdx.x & 63, wave = threadIdx.x >> 6; + const uint32_t incl = waveInclusiveScan(v); + if (lane == 63) waves[wave] = incl; + __syncthreads(); + uint32_t before = 0, all = 0; + #pragma unroll + for (int w = 0; w < 8; ++w) { + const uint32_t x = waves[w]; + before += w < wave ? x : 0u; + all += x; + } + total = all; + return before + incl - v; +} + +// 32 bytes at an arbitrarily aligned global address (in a buffer whose start is 4-byte aligned and which extends to +// the next 4-byte boundary past them). +__device__ __forceinline__ void load32(const byte* __restrict__ src, uint32_t v[8]) { + const uint32_t mis = static_cast(reinterpret_cast(src) & 3u); + const uint32_t* a = reinterpret_cast(src - mis); + if (mis == 0) { + #pragma unroll + for (int r = 0; r < 8; ++r) v[r] = a[r]; + } else { + uint32_t lo = a[0]; + #pragma unroll + for (int r = 0; r < 8; ++r) { + const uint32_t hi = a[r + 1]; + v[r] = __funnelshift_r(lo, hi, 8u * mis); + lo = hi; + } + } +} + +__device__ __forceinline__ void store32(byte* __restrict__ dst, const uint32_t v[8]) { + if ((reinterpret_cast(dst) & 3u) == 0) { + uint32_t* d = reinterpret_cast(dst); + #pragma unroll + for (int r = 0; r < 8; ++r) d[r] = v[r]; + } else { + #pragma unroll + for (int r = 0; r < 8; ++r) + #pragma unroll + for (int b = 0; b < 4; ++b) dst[4 * r + b] = static_cast(v[r] >> (8 * b)); + } +} + +__global__ void __launch_bounds__(kTPB) +encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams qp, + byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, + unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, + unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { + __shared__ __align__(16) uint32_t buf[kBufWords]; // bitshuffled planes (padded), then the packed chunk + __shared__ __align__(16) uint32_t l1w[kL1 / 4]; // L1 bitmap, then its repetition-eliminated bytes + __shared__ uint8_t re2[kL1 / 8]; // L2's repetition-eliminated bytes + __shared__ uint8_t re3[32]; // L3's + __shared__ uint32_t scratch[64]; + __shared__ uint32_t s_cid; + __shared__ uint32_t s_info[8]; // raw, csize, offset, n1, n2, n3, L4, wpos + + const int t = threadIdx.x, lane = t & 63, wave = t >> 6; + if (t == 0) s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); + __syncthreads(); + const uint32_t cid = s_cid; + const size_t base = static_cast(cid) * kElems; + + // Quantize elements i = t + 512k (QuantInplaceZigzag, verify_reconstruction): zigzag codes, raw bits for outliers. + uint32_t c[kPerT]; + #pragma unroll + for (int k = 0; k < kPerT; ++k) { + const size_t i = base + t + kTPB * k; + const float x = i < n ? in[i] : 0.0f; + const int q = static_cast(roundf(x * qp.ebx2_r)); + bool repr = q > -static_cast(qp.radius) && q < static_cast(qp.radius) && fabsf(x) < qp.threshold; + if (repr) repr = fabsf(x - fz::backend::fmulRounded(static_cast(q), qp.ebx2)) <= qp.abs_eb; + c[k] = repr ? static_cast((q << 1) ^ (q >> 31)) : __float_as_uint(x); + } + + // Difference with element i-1 (DiffNegabinary): lane-1 of this wavefront, or for lane 0 the previous + // wavefront's lane 63 (wavefront 0: wavefront 7's of slab k-1; element 0 keeps its code). + if (lane == 63) { + #pragma unroll + for (int k = 0; k < kPerT; ++k) scratch[wave * kPerT + k] = c[k]; + } + __syncthreads(); + #pragma unroll + for (int k = 0; k < kPerT; ++k) { + uint32_t prev = fz::backend::shflUp(c[k], 1, 64); + if (lane == 0) prev = wave > 0 ? scratch[(wave - 1) * kPerT + k] : (k > 0 ? scratch[7 * kPerT + k - 1] : 0u); + const int32_t d = (t == 0 && k == 0) ? static_cast(c[k]) : static_cast(c[k] - prev); + // Bitshuffle32: the butterfly leaves plane (t & 31) of group (t >> 5) + 16k in this lane. + buf[(t & 31) * kPS + (t >> 5) + 16 * k] = + fused::chunk::butterfly32(Negabinary::encode(d), t & 31); + } + __syncthreads(); + + // Zero elimination (d_ZEencodebyteshort): thread t owns bytes 32t..32t+31 of the bitshuffled chunk, i.e. plane + // t >> 4, words 8(t & 15)..+7. + uint32_t v[kPerT]; + uint32_t bmp = 0; + { + const uint32_t* row = buf + (t >> 4) * kPS + (t & 15) * 8; + #pragma unroll + for (int r = 0; r < kPerT; ++r) { + v[r] = row[r]; + #pragma unroll + for (int b = 0; b < 4; ++b) bmp |= static_cast(((v[r] >> (8 * b)) & 0xffu) != 0u) << (4 * r + b); + } + } + l1w[t] = bmp; + uint32_t wpos; + const uint32_t pos = blockExclusiveScan(__popc(bmp), scratch, wpos); + + // Wavefront 0: d_RZE's bitmap levels, its fall-back checks, then the look-back. + if (wave == 0) { + uint32_t raw = 0, csize = 2, n1 = 0, n2 = 0, n3 = 0, L4 = 0; + if (wpos >= kAvail) { + raw = 1; + } else if (wpos != 0) { + // L1 -> RE data + L2: lane owns L1 bytes 32*lane..+31; L2 bit e = (byte e != byte e-1), byte -1 = 0. + uint32_t b1[kPerT]; + #pragma unroll + for (int r = 0; r < kPerT; ++r) b1[r] = l1w[kPerT * lane + r]; + uint32_t prevByte = fz::backend::shflUp(b1[7], 1, 64) >> 24; + if (lane == 0) prevByte = 0; + uint32_t m1 = 0; + #pragma unroll + for (int r = 0; r < kPerT; ++r) + #pragma unroll + for (int b = 0; b < 4; ++b) { + const uint32_t cur = (b1[r] >> (8 * b)) & 0xffu; + m1 |= static_cast(cur != prevByte) << (4 * r + b); + prevByte = cur; + } + const uint32_t c1 = __popc(m1); + const uint32_t i1 = waveInclusiveScan(c1); + n1 = fz::backend::shfl(i1, 63, 64); + if (n1 > kAvail - wpos) { + raw = 1; + } else { + byte* l1b = reinterpret_cast(l1w); // compacted in place: positions only move down + uint32_t q = i1 - c1; + #pragma unroll + for (int r = 0; r < kPerT; ++r) + #pragma unroll + for (int b = 0; b < 4; ++b) + if ((m1 >> (4 * r + b)) & 1u) l1b[q++] = static_cast(b1[r] >> (8 * b)); + // L2 (lane owns bytes 4*lane..+3, i.e. m1) -> RE data + L3. + uint32_t prev2 = fz::backend::shflUp(m1, 1, 64) >> 24; + if (lane == 0) prev2 = 0; + uint32_t m2 = 0; + #pragma unroll + for (int b = 0; b < 4; ++b) { + const uint32_t cur = (m1 >> (8 * b)) & 0xffu; + m2 |= static_cast(cur != prev2) << b; + prev2 = cur; + } + const uint32_t c2 = __popc(m2); + const uint32_t i2 = waveInclusiveScan(c2); + n2 = fz::backend::shfl(i2, 63, 64); + if (n2 > kAvail - wpos - n1) { + raw = 1; + } else { + q = i2 - c2; + #pragma unroll + for (int b = 0; b < 4; ++b) + if ((m2 >> b) & 1u) re2[q++] = static_cast(m1 >> (8 * b)); + // L3 (32 bytes: byte j = the nibbles of lanes 2j, 2j+1) -> RE data + L4. + const uint32_t lo = fz::backend::shfl(m2, (2 * lane) & 63, 64); + const uint32_t hi = fz::backend::shfl(m2, (2 * lane + 1) & 63, 64); + const uint32_t l3 = lane < 32 ? (lo | (hi << 4)) : 0u; + uint32_t prev3 = fz::backend::shflUp(l3, 1, 64); + if (lane == 0) prev3 = 0; + const bool nr3 = lane < 32 && l3 != prev3; + L4 = static_cast(fz::backend::ballotWave(nr3)); + n3 = __popc(L4); + if (n3 > kAvail - wpos - n1 - n2) { + raw = 1; + } else { + if (nr3) re3[__popc(L4 & ((1u << lane) - 1u))] = static_cast(l3); + const uint32_t w3 = wpos + n1 + n2 + n3; + if (w3 >= kAvail - 4) raw = 1; + else csize = w3 + 4 + 2; + } + } + } + } + const uint32_t nbytes = raw ? static_cast(kChunkBytes) : csize; + const uint32_t offset = lookbackWave(cid, nbytes, state, epoch, nullptr); + if (lane == 0) { + s_info[0] = raw; s_info[1] = csize; s_info[2] = offset; s_info[3] = n1; + s_info[4] = n2; s_info[5] = n3; s_info[6] = L4; s_info[7] = wpos; + } + } + __syncthreads(); + const uint32_t raw = s_info[0], csize = s_info[1], offset = s_info[2]; + const uint32_t nbytes = raw ? static_cast(kChunkBytes) : csize; + + const unsigned long long table = kStreamHeader + kLcHeader; + const unsigned long long dst = table + 4ull * chunks + offset; + if (t == 0 && table + 4ull * (cid + 1u) <= capacity) + reinterpret_cast(out + table)[cid] = raw ? (0x80000000u | kChunkBytes) : csize; + const bool fits = dst + nbytes <= capacity; + if (raw) { + if (fits) store32(out + dst + 32u * t, v); // the bitshuffled chunk itself, as LC stores it + } else { + // Compact in place: the zero-eliminated bytes, then [RE1][RE2][RE3][L4][trailer] from wavefront 0. + byte* bb = reinterpret_cast(buf); + uint32_t q = pos; + #pragma unroll + for (int r = 0; r < kPerT; ++r) + #pragma unroll + for (int b = 0; b < 4; ++b) + if ((bmp >> (4 * r + b)) & 1u) bb[q++] = static_cast(v[r] >> (8 * b)); + if (wave == 0) { + const uint32_t n1 = s_info[3], n2 = s_info[4], n3 = s_info[5], L4 = s_info[6]; + const byte* l1b = reinterpret_cast(l1w); + for (uint32_t j = lane; j < n1; j += 64) bb[wpos + j] = l1b[j]; + for (uint32_t j = lane; j < n2; j += 64) bb[wpos + n1 + j] = re2[j]; + if (lane < static_cast(n3)) bb[wpos + n1 + n2 + lane] = re3[lane]; + if (wpos != 0 && lane < 4) bb[wpos + n1 + n2 + n3 + lane] = static_cast(L4 >> (8 * lane)); + if (lane == 0) { + bb[csize - 2] = static_cast(kChunkBytes & 0xff); + bb[csize - 1] = static_cast(kChunkBytes >> 8); + } + } + __syncthreads(); + if (fits) fused::chunk::storeChunkBytes(out + dst, bb, csize); + } + + if (cid == chunks - 1u && t == 0) { + const unsigned long long archive = kLcHeader + 4ull * chunks + offset + nbytes; + const unsigned long long padded = (archive + 3ull) & ~3ull; + const unsigned long long total = kStreamHeader + padded; + uint32_t* w = reinterpret_cast(out); + if (total <= capacity) { + for (unsigned long long b = kStreamHeader + archive; b < total; ++b) out[b] = 0; + w[4] = static_cast(chunks) * kElems * 4u; + w[5] = chunks; + } + if (kStreamHeader <= capacity) { + w[0] = kStreamMagic; + w[1] = static_cast(padded); + w[2] = __float_as_uint(qp.abs_eb); + w[3] = n; + } + if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); + } +} + +// Wavefront 0: rebuild the L1 bitmap of a packed (not raw) chunk from its tail -- [L4][trailer] at the end, RE3, +// RE2 and RE1 data before them -- into l1w, and return where the RE1 data starts (the zero-eliminated bytes' length). +// `payload` is the chunk's copy in shared memory: each level's reads depend on the one before, so they are latency. +__device__ __forceinline__ uint32_t decodeBitmaps(const byte* payload, uint32_t stored, uint32_t* __restrict__ l1w) { + const int lane = threadIdx.x & 63; + uint32_t rpos = stored - 2; + if (rpos == 0) { // an all-zero chunk: the trailer alone + for (int i = lane; i < kL1 / 4; i += 64) l1w[i] = 0; + return 0; + } + rpos -= 4; + const uint32_t L4 = static_cast(payload[rpos]) | static_cast(payload[rpos + 1]) << 8 | + static_cast(payload[rpos + 2]) << 16 | static_cast(payload[rpos + 3]) << 24; + rpos -= __popc(L4); + // L3: byte j (lanes < 32) = RE3[k - 1], k the set L4 bits in [0, j]; 0 before the first. + uint32_t l3 = 0; + if (lane < 32) { + const uint32_t k = __popc(L4 & inclMask(lane)); + l3 = k ? payload[rpos + k - 1] : 0u; + } + rpos -= fz::backend::waveSum(static_cast(__popc(l3))); + // L2: lane owns bytes 4*lane..+3, their L3 bits the nibble (lane & 1) of L3 byte lane >> 1. + const uint32_t nib = (fz::backend::shfl(l3, lane >> 1, 64) >> (4 * (lane & 1))) & 0xfu; + const uint32_t c2 = __popc(nib); + const uint32_t e2 = waveInclusiveScan(c2) - c2; + uint32_t l2 = 0; + #pragma unroll + for (int j = 0; j < 4; ++j) { + const uint32_t k = e2 + __popc(nib & inclMask(j)); + l2 |= (k ? static_cast(payload[rpos + k - 1]) : 0u) << (8 * j); + } + rpos -= fz::backend::waveSum(static_cast(__popc(l2))); + // L1: lane owns bytes 32*lane..+31, their L2 bits are l2. + const uint32_t c1 = __popc(l2); + const uint32_t e1 = waveInclusiveScan(c1) - c1; + #pragma unroll + for (int r = 0; r < kPerT; ++r) { + uint32_t word = 0; + #pragma unroll + for (int b = 0; b < 4; ++b) { + const uint32_t k = e1 + __popc(l2 & inclMask(4 * r + b)); + word |= (k ? static_cast(payload[rpos + k - 1]) : 0u) << (8 * b); + } + l1w[kPerT * lane + r] = word; + } + return rpos; +} + +__global__ void __launch_bounds__(kTPB) +decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t n, + unsigned long long* __restrict__ state, uint32_t epoch) { + __shared__ __align__(16) uint32_t buf[kBufWords]; // the packed chunk, then the bitshuffled planes (padded) + __shared__ __align__(16) uint32_t l1w[kL1 / 4]; + __shared__ uint32_t scratch[64]; + __shared__ uint32_t s_offset; + + const int t = threadIdx.x, lane = t & 63, wave = t >> 6; + const uint32_t* w = reinterpret_cast(stream); + const uint32_t chunks = gridDim.x; + if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || + w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { + if (blockIdx.x == 0 && t == 0) + printf("[fzgmod] pfpl::decompress: not a PFPL stream of %u floats " + "(magic %08x, n %u, chunks %u, archive %u B)\n", n, w[0], w[3], w[5], w[1]); + streamTrap(); + } + const uint32_t cid = blockIdx.x; + const uint32_t* table = w + (kStreamHeader + kLcHeader) / 4u; + const uint32_t entry = table[cid]; + const uint32_t stored = entry & 0x7fffffffu; + const bool raw = (entry & 0x80000000u) != 0u; + if (wave == 0) { + const uint32_t offset = lookbackWave(cid, stored, state, epoch, table); + if (lane == 0) s_offset = offset; + } + __syncthreads(); + const byte* payload = stream + kStreamHeader + kLcHeader + 4ull * chunks + s_offset; + + uint32_t v[kPerT]; + if (raw) { + load32(payload + 32u * t, v); + } else { + fused::chunk::loadChunkBytes(reinterpret_cast(buf), payload, stored); + __syncthreads(); + if (wave == 0) (void)decodeBitmaps(reinterpret_cast(buf), stored, l1w); + __syncthreads(); + const uint32_t bmp = l1w[t]; + uint32_t total; + uint32_t q = blockExclusiveScan(__popc(bmp), scratch, total); + const byte* bb = reinterpret_cast(buf); + #pragma unroll + for (int r = 0; r < kPerT; ++r) { + uint32_t word = 0; + #pragma unroll + for (int b = 0; b < 4; ++b) + if ((bmp >> (4 * r + b)) & 1u) word |= static_cast(bb[q++]) << (8 * b); + v[r] = word; + } + __syncthreads(); // every gather done before the planes overwrite the packed bytes + } + { + uint32_t* row = buf + (t >> 4) * kPS + (t & 15) * 8; + #pragma unroll + for (int r = 0; r < kPerT; ++r) row[r] = v[r]; + } + __syncthreads(); + + // Inverse bitshuffle and negabinary: element i = t + 512k; then its chunk-local inclusive prefix sum. + uint32_t d[kPerT]; + #pragma unroll + for (int k = 0; k < kPerT; ++k) + d[k] = static_cast(Negabinary::decode( + fused::chunk::butterfly32(buf[(t & 31) * kPS + (t >> 5) + 16 * k], t & 31))); + #pragma unroll + for (int dd = 1; dd < 64; dd <<= 1) + #pragma unroll + for (int k = 0; k < kPerT; ++k) { + const uint32_t up = fz::backend::shflUp(d[k], dd, 64); + if (lane >= dd) d[k] += up; + } + if (lane == 63) { + #pragma unroll + for (int k = 0; k < kPerT; ++k) scratch[k * 8 + wave] = d[k]; + } + __syncthreads(); + const size_t base = static_cast(cid) * kElems; + const uint32_t cnt = min(static_cast(kElems), n - cid * kElems); + const float ebx2 = 2.0f * __uint_as_float(w[2]); + uint32_t carry = 0; + #pragma unroll + for (int k = 0; k < kPerT; ++k) { + uint32_t before = 0, slab = 0; + #pragma unroll + for (int ww = 0; ww < 8; ++ww) { + const uint32_t x = scratch[k * 8 + ww]; + before += ww < wave ? x : 0u; + slab += x; + } + const uint32_t code = d[k] + before + carry; + carry += slab; + const uint32_t e = t + kTPB * k; + if (e < cnt) { + float y; + if ((code >> 1) >= kQuantRadius) { + y = __uint_as_float(code); + } else { + const int32_t qq = static_cast((code >> 1) ^ (0u - (code & 1u))); + y = static_cast(qq) * ebx2; + } + out[base + e] = y; + } + } +} + +} // namespace amd diff --git a/modules/quantizers/quantizer/quantizer.cu b/modules/quantizers/quantizer/quantizer.cu index 16b0a508..c49df4ee 100644 --- a/modules/quantizers/quantizer/quantizer.cu +++ b/modules/quantizers/quantizer/quantizer.cu @@ -20,6 +20,7 @@ #include "cuda_check.h" #include "log.h" #include +#include "backend/fp.h" namespace fz { @@ -168,11 +169,12 @@ __global__ void quantizer_abs_fwd_inplace_kernel( && quantAbs(x) < static_cast(threshold); if (representable && verify_reconstruction) { // Keep the multiply rounded as a standalone float operation, matching - // PFPL's `recon = bin * eb2` before its subsequent subtraction. Without - // the explicit intrinsic nvcc may contract the multiply into the error - // subtraction and classify a different set of boundary values. + // PFPL's `recon = bin * eb2` before its subsequent subtraction. Fused + // into the subtraction (nvcc may, and HIP's clang does by default), the + // check tests the exact product rather than the stored one and passes + // values whose reconstruction misses the bound (backend/fp.h). const TInput recon = std::is_same::value - ? static_cast(__fmul_rn(static_cast(q), static_cast(ebx2))) + ? static_cast(fz::backend::fmulRounded(static_cast(q), static_cast(ebx2))) : static_cast(q) * ebx2; representable = quantAbs(x - recon) <= abs_eb; } diff --git a/modules/shufflers/bitshuffle/bitshuffle_stage.cu b/modules/shufflers/bitshuffle/bitshuffle_stage.cu index a29b8327..459324a7 100644 --- a/modules/shufflers/bitshuffle/bitshuffle_stage.cu +++ b/modules/shufflers/bitshuffle/bitshuffle_stage.cu @@ -465,7 +465,7 @@ void BitshuffleStage::execute( const size_t smem = bitshuf32SmemBytes(N_chunk); if (bitshufUseSmem() && smem <= (size_t)bitshufMaxSmem()) { if (smem > 48u * 1024u) - cudaFuncSetAttribute(bitshuffleEncodeKernel32Smem, + cudaFuncSetAttribute(reinterpret_cast(bitshuffleEncodeKernel32Smem), cudaFuncAttributeMaxDynamicSharedMemorySize, (int)smem); bitshuffleEncodeKernel32Smem <<>>( @@ -515,7 +515,7 @@ void BitshuffleStage::execute( const size_t smem = bitshuf32SmemBytes(N_chunk); if (bitshufUseSmem() && smem <= (size_t)bitshufMaxSmem()) { if (smem > 48u * 1024u) - cudaFuncSetAttribute(bitshuffleDecodeKernel32Smem, + cudaFuncSetAttribute(reinterpret_cast(bitshuffleDecodeKernel32Smem), cudaFuncAttributeMaxDynamicSharedMemorySize, (int)smem); bitshuffleDecodeKernel32Smem <<>>( diff --git a/src/pipeline/compressor.cpp b/src/pipeline/compressor.cpp index 1ef7780f..7d2e5d71 100644 --- a/src/pipeline/compressor.cpp +++ b/src/pipeline/compressor.cpp @@ -206,6 +206,10 @@ void Pipeline::planAndInstallFusion() { fusion_info_.fallback_reason = "policy_off"; return; } +#if defined(FZGMOD_BACKEND_HIP) + fusion_info_.fallback_reason = "unsupported_backend"; + return; +#endif std::vector installed; const bool include_experimental = mode == FusionPolicy::Force; diff --git a/src/pipeline/fusion_registry.cpp b/src/pipeline/fusion_registry.cpp index 2e9d747e..f4675d3b 100644 --- a/src/pipeline/fusion_registry.cpp +++ b/src/pipeline/fusion_registry.cpp @@ -21,6 +21,7 @@ namespace fz { +#if !defined(FZGMOD_BACKEND_HIP) // the fused implementations are NVRTC/CUDA-only; see findFusedImpl() namespace { // ── Warp-register (cuSZp): Quantizer(linear,ABS) -> block-local predictor -> @@ -536,13 +537,20 @@ const FusedImpl kBuiltins[] = { }; } // namespace +#endif // !FZGMOD_BACKEND_HIP const FusedImpl* findFusedImpl( const std::vector& group, bool include_experimental) { +#if defined(FZGMOD_BACKEND_HIP) + // No fused implementation exists on HIP (they JIT with NVRTC), so specialization falls back to staged. + (void)group; (void)include_experimental; + return nullptr; +#else for (const auto& impl : kBuiltins) if ((impl.auto_enabled || include_experimental) && impl.matches(group)) return &impl; return nullptr; +#endif } } // namespace fz From 6d223860d80b21d84f92de6a4c76a022291d6bdc Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Wed, 4 Nov 2026 15:52:50 -0500 Subject: [PATCH 2/8] Restore the HIP backend on main -DFZGMOD_BACKEND=HIP configures and builds again. The root CMakeLists.txt, presets and package config take their HIP hunks from development, without its profiling target, which main does not have. Module and test CMake files already branch on the backend. NVRTC specialization is not built on HIP, so every stage runs staged. Four tests need it or include cuda_runtime.h directly, so they are CUDA-only. On one MI250X GCD, 50 of 53 tests pass. test_huffman, test_log_transform and test_stage_registry fail on HIP, and they are not changed here. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 3 + CMakeLists.txt | 130 ++++++++++++++++++++++++---- CMakePresets.json | 23 +++++ cmake/FZGPUModulesConfig.cmake.in | 15 +++- modules/fused/common/CMakeLists.txt | 15 ++++ tests/pipeline/CMakeLists.txt | 45 ++++++---- tests/stages/CMakeLists.txt | 15 ++-- 7 files changed, 204 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 964936da..96286d98 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ Version numbers follow [Semantic Versioning](https://semver.org/). ## [Unreleased] +### Added +- `-DFZGMOD_BACKEND=HIP` configures again, with `hip` CMake presets. NVRTC specialization is not built on HIP, so every stage runs staged, and `test_fusion_f64`, `test_fusion_fixed`, `test_fusion_planner` and `test_lorenzo_quant_shared_init` are CUDA-only. + ### Changed - CMake now locates the build-time Python interpreter with `find_package(Python3 3.7 REQUIRED COMPONENTS Interpreter)` and runs the header generators through `${Python3_EXECUTABLE}`, so a missing or too-old Python fails at configure time instead of mid-build. Override with `-DPython3_EXECUTABLE=...`. diff --git a/CMakeLists.txt b/CMakeLists.txt index 1f863647..93112a60 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,15 +1,26 @@ cmake_minimum_required(VERSION 3.24 FATAL_ERROR) -set(FZGMOD_BACKEND "CUDA" CACHE STRING "GPU backend (stable main supports CUDA)") -set_property(CACHE FZGMOD_BACKEND PROPERTY STRINGS CUDA) -if(NOT FZGMOD_BACKEND STREQUAL "CUDA") - message(FATAL_ERROR "Stable main supports CUDA only; backend development is on the development branch.") -endif() -if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) - set(CMAKE_CUDA_ARCHITECTURES "86" CACHE STRING "CUDA architectures to compile for") +set(FZGMOD_BACKEND "CUDA" CACHE STRING "GPU backend to build against: CUDA, HIP, or SYCL") +set_property(CACHE FZGMOD_BACKEND PROPERTY STRINGS CUDA HIP SYCL) + +if(FZGMOD_BACKEND STREQUAL "CUDA") + set(FZGMOD_PROJECT_LANGUAGES CUDA CXX) + if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) + set(CMAKE_CUDA_ARCHITECTURES "86" CACHE STRING "CUDA architectures to compile for") + endif() +elseif(FZGMOD_BACKEND STREQUAL "HIP") + set(FZGMOD_PROJECT_LANGUAGES CXX) + message(WARNING + "FZGMOD_BACKEND=HIP is experimental: backend coverage and install/consumer " + "compatibility are still being completed, and some stages are unavailable.") +elseif(FZGMOD_BACKEND STREQUAL "SYCL") + message(FATAL_ERROR "FZGMOD_BACKEND=SYCL is not implemented yet.") +else() + message(FATAL_ERROR "Unknown FZGMOD_BACKEND '${FZGMOD_BACKEND}' — must be CUDA, HIP, or SYCL.") endif() -project(FZGPUModules LANGUAGES CUDA CXX VERSION 2.0.1) + +project(FZGPUModules LANGUAGES ${FZGMOD_PROJECT_LANGUAGES} VERSION 2.0.1) # General settings @@ -39,18 +50,76 @@ endif() set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED True) -# CUDA toolchain -find_package(CUDAToolkit REQUIRED) -if(CUDAToolkit_VERSION VERSION_LESS "11.2") - message(FATAL_ERROR "FZGPUModules requires CUDA Toolkit 11.2 or newer for stream-ordered allocation.") +# GPU backend toolchain +# Link the selected runtime through FZGMOD_GPU_RUNTIME. + +if(FZGMOD_BACKEND STREQUAL "CUDA") + find_package(CUDAToolkit REQUIRED) + if(CUDAToolkit_VERSION VERSION_LESS "11.2") + message(FATAL_ERROR + "FZGPUModules requires CUDA Toolkit 11.2 or newer.\n" + " Reason: stream-ordered allocators (cudaMallocAsync / cudaMemPool_t) require 11.2+.\n" + " Detected: ${CUDAToolkit_VERSION}") + endif() + + # Set CUDA architectures before project() enables the language. + enable_language(CUDA) + + set(FZGMOD_GPU_RUNTIME CUDA::cudart) + +else() # HIP — SYCL and unknown values already rejected above + + # Set HIP architectures before enabling the language; default to MI100. + if(NOT DEFINED CMAKE_HIP_ARCHITECTURES) + set(CMAKE_HIP_ARCHITECTURES "gfx908" CACHE STRING + "AMD GPU architectures to compile for (gfx908=MI100, gfx90a=MI200, gfx942=MI300)") + endif() + + enable_language(HIP) + + # Pass explicit targets to hip-config, including on GPU-less build nodes. + set(GPU_TARGETS "${CMAKE_HIP_ARCHITECTURES}" CACHE STRING "AMD GPU targets for hip-config" FORCE) + + # Put the ROCm root on CMAKE_PREFIX_PATH for CONFIG-mode discovery. + find_package(hip CONFIG REQUIRED) + find_package(hipcub CONFIG REQUIRED) # cub facade — mirrors cub's API + find_package(rocthrust CONFIG REQUIRED) # thrust facade, used by chunk_fusion/rze/golomb_rice + + # Use hip::host: hip::device would apply -x hip to plain C++ sources. + # Details: docs/codebase_notes.md CN-BUILD-1. + set(FZGMOD_GPU_RUNTIME hip::host) + + # Use header-only include paths to avoid transitively linking hip::device. + set(FZGMOD_HIP_HEADER_ONLY_INCLUDES "") + foreach(_hdr_tgt hip::hipcub roc::rocprim roc::rocthrust) + if(TARGET ${_hdr_tgt}) + get_target_property(_inc ${_hdr_tgt} INTERFACE_INCLUDE_DIRECTORIES) + if(_inc) + list(APPEND FZGMOD_HIP_HEADER_ONLY_INCLUDES ${_inc}) + endif() + endif() + endforeach() + list(REMOVE_DUPLICATES FZGMOD_HIP_HEADER_ONLY_INCLUDES) + endif() -set(FZGMOD_GPU_RUNTIME CUDA::cudart) # Build-time header generators (stage fingerprints, embedded NVRTC headers). find_package(Python3 3.7 REQUIRED COMPONENTS Interpreter) -# Retain the module-local source hook; CUDA needs no retagging. +# Retag .cu sources as HIP; define this before subdirectories use it. +# Backend rationale: docs/codebase_notes.md CN-BUILD-1. function(fzgmod_set_gpu_source_language target) + if(NOT FZGMOD_BACKEND STREQUAL "HIP") + return() + endif() + get_target_property(_srcs ${target} SOURCES) + foreach(_src IN LISTS _srcs) + if(_src MATCHES "\\.cu$") + set_source_files_properties(${_src} + TARGET_DIRECTORY ${target} + PROPERTIES LANGUAGE HIP) + endif() + endforeach() endfunction() # WSL detection @@ -135,7 +204,7 @@ elseif (USE_SANITIZER STREQUAL "Compute") if(NOT FZGMOD_BACKEND STREQUAL "CUDA") message(FATAL_ERROR "USE_SANITIZER=Compute requires the CUDA backend (compute-sanitizer is an " - "NVIDIA tool).") + "NVIDIA tool). Sanitizer parity for HIP is not implemented.") endif() message(STATUS "Configuring for Compute Sanitizer") set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -lineinfo -Xcompiler -rdynamic") @@ -163,7 +232,20 @@ target_compile_options( ) target_compile_features(fzgmod_compile_settings INTERFACE cxx_std_17) -target_compile_features(fzgmod_compile_settings INTERFACE cuda_std_17) +if(FZGMOD_BACKEND STREQUAL "CUDA") + target_compile_features(fzgmod_compile_settings INTERFACE cuda_std_17) +else() + set(CMAKE_HIP_STANDARD 17) + set(CMAKE_HIP_STANDARD_REQUIRED True) + + # HIP teardown paths intentionally ignore [[nodiscard]] runtime status. + target_compile_options(fzgmod_compile_settings INTERFACE + $<$:-Wno-unused-result>) + + # ROCm requires this flag for the sync warp intrinsics. + target_compile_definitions(fzgmod_compile_settings INTERFACE + HIP_ENABLE_WARP_SYNC_BUILTINS) +endif() target_link_libraries(fzgmod_compile_settings INTERFACE ${FZGMOD_GPU_RUNTIME}) @@ -188,8 +270,18 @@ target_include_directories(fzgmod_compile_settings INTERFACE $ ) +# Treat HIP dependency headers as system headers. +if(FZGMOD_HIP_HEADER_ONLY_INCLUDES) + target_include_directories(fzgmod_compile_settings SYSTEM INTERFACE + ${FZGMOD_HIP_HEADER_ONLY_INCLUDES}) +endif() + +# Production annotations can be enabled without building research programs. option(ENABLE_NVTX "Enable NVTX annotations in production library code" OFF) if(ENABLE_NVTX) + if(NOT FZGMOD_BACKEND STREQUAL "CUDA") + message(FATAL_ERROR "ENABLE_NVTX requires the CUDA backend.") + endif() find_path(NVTX3_INCLUDE_DIR NAMES nvtx3/nvtx3.hpp HINTS ${CUDAToolkit_INCLUDE_DIRS} REQUIRED) target_compile_definitions(fzgmod_compile_settings INTERFACE @@ -456,6 +548,10 @@ endif() message(STATUS "FZGPUModules configuration complete") message(STATUS " Backend: ${FZGMOD_BACKEND}") message(STATUS " Build type: ${CMAKE_BUILD_TYPE}") -message(STATUS " CUDA architectures: ${CMAKE_CUDA_ARCHITECTURES}") +if(FZGMOD_BACKEND STREQUAL "CUDA") + message(STATUS " CUDA architectures: ${CMAKE_CUDA_ARCHITECTURES}") +else() + message(STATUS " HIP architectures: ${CMAKE_HIP_ARCHITECTURES}") +endif() message(STATUS " Install prefix: ${CMAKE_INSTALL_PREFIX}") message(STATUS " Build examples: ${BUILD_EXAMPLES}") diff --git a/CMakePresets.json b/CMakePresets.json index e3461b3e..37e0d4a5 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -44,6 +44,16 @@ } }, + { + "name": "hip", + "inherits": "base", + "displayName": "HIP (AMD ROCm, experimental)", + "description": "Experimental release build against the HIP backend. Defaults to gfx908 (MI100); override with -DCMAKE_HIP_ARCHITECTURES=. Requires the ROCm root on CMAKE_PREFIX_PATH, e.g. export CMAKE_PREFIX_PATH=/opt/rocm:$CMAKE_PREFIX_PATH", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "FZGMOD_BACKEND": "HIP" + } + }, { "name": "cuda-h200", @@ -87,6 +97,11 @@ "configurePreset": "asan", "displayName": "Build ASan + UBSan" }, + { + "name": "hip", + "configurePreset": "hip", + "displayName": "Build HIP (AMD ROCm, experimental)" + }, { "name": "cuda-h200", "configurePreset": "cuda-h200", @@ -129,6 +144,14 @@ "filter": { "include": { "label": "pipeline" } } }, + { + "name": "hip", + "configurePreset": "hip", + "displayName": "All tests (HIP, experimental)", + "description": "Run the full test suite against the experimental HIP build. ANSStage is excluded on HIP (see Stage::isSupportedOnBackend).", + "execution": { "noTestsAction": "error" }, + "output": { "outputOnFailure": true } + }, { "name": "cuda-h200", diff --git a/cmake/FZGPUModulesConfig.cmake.in b/cmake/FZGPUModulesConfig.cmake.in index b1743cef..d9d0119a 100644 --- a/cmake/FZGPUModulesConfig.cmake.in +++ b/cmake/FZGPUModulesConfig.cmake.in @@ -4,7 +4,20 @@ include(CMakeFindDependencyMacro) set(FZGPUModules_BACKEND "@FZGMOD_BACKEND@") -find_dependency(CUDAToolkit REQUIRED) +if(FZGPUModules_BACKEND STREQUAL "CUDA") + find_dependency(CUDAToolkit REQUIRED) +elseif(FZGPUModules_BACKEND STREQUAL "HIP") + # HIP support is experimental. These packages provide the runtime target and + # public headers referenced by the exported FZGPUModules targets. + find_dependency(hip CONFIG REQUIRED) + find_dependency(hipcub CONFIG REQUIRED) + find_dependency(rocthrust CONFIG REQUIRED) +else() + set(FZGPUModules_FOUND FALSE) + set(FZGPUModules_NOT_FOUND_MESSAGE + "This FZGPUModules package was built with unknown backend '${FZGPUModules_BACKEND}'.") + return() +endif() include("${CMAKE_CURRENT_LIST_DIR}/FZGPUModulesTargets.cmake") diff --git a/modules/fused/common/CMakeLists.txt b/modules/fused/common/CMakeLists.txt index cb6cf4fa..082b7c33 100644 --- a/modules/fused/common/CMakeLists.txt +++ b/modules/fused/common/CMakeLists.txt @@ -18,3 +18,18 @@ if(NOT FZGMOD_BACKEND STREQUAL "HIP") add_dependencies(fzgmod_modules fzgmod_nvrtc_headers) target_include_directories(fzgmod_modules PRIVATE "${PROJECT_BINARY_DIR}/generated") endif() + +# Pipeline specialization JIT-compiles its fused kernels with NVRTC and launches them through the CUDA driver API, and its +# harness assumes 32-lane warps: none of it is built on HIP, where the registry offers no fused implementation +# (src/pipeline/fusion_registry.cpp) and every stage runs staged. +if(FZGMOD_BACKEND STREQUAL "HIP") + set_source_files_properties( + ${CMAKE_CURRENT_LIST_DIR}/../fused_block/fused_block.cu + ${CMAKE_CURRENT_LIST_DIR}/../fused_block/nvrtc_warp_fusion.cu + ${CMAKE_CURRENT_LIST_DIR}/../chunk_fusion/chunk_fusion.cu + ${CMAKE_CURRENT_LIST_DIR}/../chunk_fusion/nvrtc_chunk_fusion.cpp + ${CMAKE_CURRENT_LIST_DIR}/nvrtc_jit.cpp + TARGET_DIRECTORY fzgmod_modules + PROPERTIES HEADER_FILE_ONLY TRUE + ) +endif() diff --git a/tests/pipeline/CMakeLists.txt b/tests/pipeline/CMakeLists.txt index a195f749..b0aaa3f7 100644 --- a/tests/pipeline/CMakeLists.txt +++ b/tests/pipeline/CMakeLists.txt @@ -58,28 +58,37 @@ fz_add_test( ) # ─── Fusion planner (identifies fusable stage chains in the DAG) ─────────────── -fz_add_test( - NAME test_fusion_planner - SOURCES test_fusion_planner.cpp - LIBS fzgmod - LABELS pipeline gpu -) +# Calls the NVRTC source generators, which only the CUDA backend builds. +if(FZGMOD_BACKEND STREQUAL "CUDA") + fz_add_test( + NAME test_fusion_planner + SOURCES test_fusion_planner.cpp + LIBS fzgmod + LABELS pipeline gpu + ) +endif() # ─── f64 Auto specialization correctness and fallback coverage ──────────────── -fz_add_test( - NAME test_fusion_f64 - SOURCES test_fusion_f64.cpp - LIBS fzgmod - LABELS pipeline gpu -) +# Calls the NVRTC source generators, which only the CUDA backend builds. +if(FZGMOD_BACKEND STREQUAL "CUDA") + fz_add_test( + NAME test_fusion_f64 + SOURCES test_fusion_f64.cpp + LIBS fzgmod + LABELS pipeline gpu + ) +endif() # ─── Predictor-free fixed 1D warp specialization ────────────────────────────── -fz_add_test( - NAME test_fusion_fixed - SOURCES test_fusion_fixed.cpp - LIBS fzgmod - LABELS pipeline gpu -) +# Requires an installed specialization, which only the CUDA backend has. +if(FZGMOD_BACKEND STREQUAL "CUDA") + fz_add_test( + NAME test_fusion_fixed + SOURCES test_fusion_fixed.cpp + LIBS fzgmod + LABELS pipeline gpu + ) +endif() # ─── Config file load/save (TOML) ───────────────────────────────────────────── fz_add_test( diff --git a/tests/stages/CMakeLists.txt b/tests/stages/CMakeLists.txt index 46453039..7357a5bb 100644 --- a/tests/stages/CMakeLists.txt +++ b/tests/stages/CMakeLists.txt @@ -189,12 +189,15 @@ fz_add_test( ) # 3-D LorenzoQuant forward kernel must not depend on stale shared memory -fz_add_test( - NAME test_lorenzo_quant_shared_init - SOURCES test_lorenzo_quant_shared_init.cu - LIBS fzgmod_modules fzgmod_mem fzgmod - LABELS stages gpu -) +# Includes cuda_runtime.h directly. +if(FZGMOD_BACKEND STREQUAL "CUDA") + fz_add_test( + NAME test_lorenzo_quant_shared_init + SOURCES test_lorenzo_quant_shared_init.cu + LIBS fzgmod_modules fzgmod_mem fzgmod + LABELS stages gpu + ) +endif() # ─── Per-tile adaptive multi-order Lorenzo (FSZ prediction stage) ──────────── fz_add_test( From 07919bfc68f93fd16b71c831a57639f40335da45 Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Wed, 4 Nov 2026 15:52:50 -0500 Subject: [PATCH 3/8] fz::pfpl: clear a new Plan on the calling stream; add tests and notes Plan::reserve cleared its ticket counter and look-back state with a stream-less cudaMemset, which is not ordered with a non-blocking stream's kernels. The first encode could read old page contents as its ticket, compute a wild chunk id and stop with a memory aperture violation. The clears now run on the calling stream, ahead of that call's kernel. See CN-PFPLSTREAM-2. tests/stages/test_pfpl_stream.cpp compares fz::pfpl with the staged pipeline. On HIP, CTest runs it with the AMD kernels and with the generic harness kernels. The device-side printf in the decode error paths is gone, and the measurements moved from code comments to CN-PFPLSTREAM-1. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 + docs/codebase_notes.md | 43 +++ modules/fused/chunk_fusion/pfpl_stream.cu | 18 +- modules/fused/chunk_fusion/pfpl_stream.h | 5 +- .../fused/chunk_fusion/pfpl_stream_amd.inl | 9 +- tests/stages/CMakeLists.txt | 16 + tests/stages/test_pfpl_stream.cpp | 363 ++++++++++++++++++ 7 files changed, 448 insertions(+), 12 deletions(-) create mode 100644 tests/stages/test_pfpl_stream.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 96286d98..0eec188b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,9 +10,15 @@ Version numbers follow [Semantic Versioning](https://semver.org/). ## [Unreleased] ### Added +- `fz::pfpl` (`fused/chunk_fusion/pfpl_stream.h`): a float field compresses from one device buffer into another and back in one stream-ordered kernel each, with no host synchronization and the same archive bytes as the PFPL preset in ABS mode. +- Wave64 AMD kernels for `fz::pfpl`, used by default on HIP; `FZ_PFPL_HARNESS=1` selects the generic harness kernels. See CN-PFPLSTREAM-1. - `-DFZGMOD_BACKEND=HIP` configures again, with `hip` CMake presets. NVRTC specialization is not built on HIP, so every stage runs staged, and `test_fusion_f64`, `test_fusion_fixed`, `test_fusion_planner` and `test_lorenzo_quant_shared_init` are CUDA-only. +### Fixed +- On HIP, the ABS bound check could accept values that decode up to half an ulp past the bound, because clang fused its multiply and subtract into one FMA. `fz::backend::fmulRounded` keeps the product rounded in the staged quantizer and the chunk harness. + ### Changed +- The chunk harness decode copies payloads with `loadChunkBytes` instead of byte by byte, on CUDA as well. - CMake now locates the build-time Python interpreter with `find_package(Python3 3.7 REQUIRED COMPONENTS Interpreter)` and runs the header generators through `${Python3_EXECUTABLE}`, so a missing or too-old Python fails at configure time instead of mid-build. Override with `-DPython3_EXECUTABLE=...`. ## [2.0.1] — 2026-10-15 diff --git a/docs/codebase_notes.md b/docs/codebase_notes.md index f25c9a05..4211775a 100644 --- a/docs/codebase_notes.md +++ b/docs/codebase_notes.md @@ -1816,3 +1816,46 @@ it in `stage_registry.h`, and add FZGM integration tests. Use the existing SDK capability guards for optional families. Registration order retains the previous first-match behavior; individual handlers own their token aliases. The main plugin requires no stage-specific edits. + +## CN-PFPLSTREAM-1 — Why `fz::pfpl` has AMD kernels, and what they gain + +**Source:** `modules/fused/chunk_fusion/pfpl_stream_amd.inl`. + +The generic harness kernels keep a chunk in shared memory with block-wide scans +and about 20 barriers per chunk. Profiled on MI250X, they wait about 65% of the +time at 8 waves per CU. Their 36 KiB of shared memory also limits a CU to one +block. The AMD kernels produce the same stream, byte for byte, with 6-7 barriers +and 19 KiB of shared memory, which fits 3 blocks per CU. + +Measured on one MI250X GCD with ROCm 7.1.1, on a polar tile with `abs_eb` 1e-4 +and a compression ratio of about 1.7. Throughput is in GB/s: + +| Field | Encode: two-array | Encode: staged | Decode: two-array | Decode: staged | +|---|---:|---:|---:|---:| +| 1 MiB | 57.3 | 6.7 | 46.1 | 6.8 | +| 8 MiB | 151.6 | 25.7 | 129.2 | 29.2 | +| 32 MiB | 191.7 | 37.1 | 156.9 | 45.4 | +| 128 MiB | 203.3 | 41.8 | 165.6 | 50.3 | +| 512 MiB | 202.5 | 43.0 | 167.7 | 52.0 | + +## CN-PFPLSTREAM-2 — A new `Plan` clears its memory on the calling stream + +**Source:** `Plan::reserve` and `Plan::nextEpoch` in +`modules/fused/chunk_fusion/pfpl_stream.cu`. + +`Plan::reserve` used to clear its new ticket counter and look-back state with a +stream-less `cudaMemset`. On MI250X with ROCm 7.1.1 that clear could land after +the first encode kernel on a non-blocking stream had started. A read on the +calling stream right after the call returned the page's previous contents, and +the same read after a 50 ms wait, with no GPU synchronization, returned zero. + +A stale counter gives `atomicAdd(ticket, 1) - ticket_base` a wild chunk id, so the +kernel reads from `in + cid * kElems`. Both the AMD and the harness encode kernels +stopped with `HSA_STATUS_ERROR_MEMORY_APERTURE_VIOLATION`. Stale look-back words +tagged with the same epoch would give a wrong offset with no fault. + +The clears now run in `Plan::nextEpoch` with `cudaMemsetAsync` on the calling +stream, ahead of that call's kernel, so they need no host synchronization. +`PfplStream.FreshPlanIgnoresStaleMemory` keeps the default stream busy so that a +clear queued there always loses the race. Without this fix it faults in its first +round on both kernel sets. diff --git a/modules/fused/chunk_fusion/pfpl_stream.cu b/modules/fused/chunk_fusion/pfpl_stream.cu index 2e424447..cf70a5c3 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.cu +++ b/modules/fused/chunk_fusion/pfpl_stream.cu @@ -200,9 +200,6 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t const uint32_t chunks = gridDim.x; if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { - if (blockIdx.x == 0 && threadIdx.x == 0) - printf("[fzgmod] pfpl::decompress: not a PFPL stream of %u floats " - "(magic %08x, n %u, chunks %u, archive %u B)\n", n, w[0], w[3], w[5], w[1]); streamTrap(); } const uint32_t cid = blockIdx.x; @@ -274,7 +271,7 @@ void Plan::reserve(size_t n) { const size_t chunks = std::max(chunksOf(n), 1); if (ticket_ == nullptr) { FZ_CUDA_CHECK(cudaMalloc(&ticket_, sizeof(unsigned long long))); - FZ_CUDA_CHECK(cudaMemset(ticket_, 0, sizeof(unsigned long long))); + clear_ticket_ = true; tickets_ = 0; } if (chunks <= chunks_) return; @@ -284,12 +281,23 @@ void Plan::reserve(size_t n) { state_ = nullptr; } FZ_CUDA_CHECK(cudaMalloc(&state_, chunks * sizeof(unsigned long long))); - FZ_CUDA_CHECK(cudaMemset(state_, 0, chunks * sizeof(unsigned long long))); + clear_state_ = true; chunks_ = chunks; epoch_ = 0; // a cleared state: every epoch is new again } uint32_t Plan::nextEpoch(fz::stream_t stream) { + // New memory can hold earlier bytes: a stale ticket is a wild chunk id, and stale look-back words are a wrong + // offset. Clear it on the calling stream, because a stream-less cudaMemset is not ordered with a non-blocking + // stream's kernels. See CN-PFPLSTREAM-2. + if (clear_ticket_) { + FZ_CUDA_CHECK(cudaMemsetAsync(ticket_, 0, sizeof(unsigned long long), stream)); + clear_ticket_ = false; + } + if (clear_state_) { + FZ_CUDA_CHECK(cudaMemsetAsync(state_, 0, chunks_ * sizeof(unsigned long long), stream)); + clear_state_ = false; + } if (++epoch_ > kEpochMask) { // after 2^30 calls, start over on a cleared state FZ_CUDA_CHECK(cudaMemsetAsync(state_, 0, chunks_ * sizeof(unsigned long long), stream)); epoch_ = 1; diff --git a/modules/fused/chunk_fusion/pfpl_stream.h b/modules/fused/chunk_fusion/pfpl_stream.h index a8afe931..7d62067e 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.h +++ b/modules/fused/chunk_fusion/pfpl_stream.h @@ -61,7 +61,8 @@ class Plan { Plan(const Plan&) = delete; Plan& operator=(const Plan&) = delete; - /** Make room for fields of up to n floats. Allocating synchronizes the device. */ + /** Make room for fields of up to n floats. Replacing a smaller state synchronizes the device; new memory is + * cleared on the stream of the next compress or decompress. */ void reserve(size_t n); private: @@ -74,6 +75,8 @@ class Plan { size_t chunks_ = 0; ///< chunks state_ holds unsigned long long tickets_ = 0; ///< the ticket counter's value between calls uint32_t epoch_ = 0; ///< last epoch used (0 = never: fresh state reads as empty) + bool clear_ticket_ = false; ///< ticket_ is new: zero it ahead of the next call's kernel + bool clear_state_ = false; ///< state_ is new: zero it ahead of the next call's kernel }; /** diff --git a/modules/fused/chunk_fusion/pfpl_stream_amd.inl b/modules/fused/chunk_fusion/pfpl_stream_amd.inl index 8d6b1327..ce2b0dfe 100644 --- a/modules/fused/chunk_fusion/pfpl_stream_amd.inl +++ b/modules/fused/chunk_fusion/pfpl_stream_amd.inl @@ -1,9 +1,9 @@ // Included by pfpl_stream.cu (inside fz::pfpl::): the AMD (wave64) encode/decode kernels. // // The generic harness kernels keep a chunk in two 16 KiB shared buffers plus 4 KiB of coder scratch (36 KiB: one -// 512-thread block per CU on MI200's 64 KiB LDS) and run d_RZE's bitmap recursion with block-wide scans (~20 -// barriers per chunk); profiled on MI250X they wait ~65% of the time at 8 waves per CU. These produce the same -// stream, byte for byte, from one padded 16.1 KiB buffer and 2 KiB of bitmap (3 blocks per CU), with 6-7 barriers: +// 512-thread block per CU on MI200's 64 KiB LDS) and run d_RZE's bitmap recursion with block-wide scans, which +// spend most of their time waiting at barriers (CN-PFPLSTREAM-1). These produce the same stream, byte for byte, from +// one padded 16.1 KiB buffer and 2 KiB of bitmap (3 blocks per CU), with 6-7 barriers: // - quantize, difference and bitshuffle in registers (the difference's cross-wavefront neighbours through a // 64-word exchange), scattering the bitshuffled planes into a buffer padded to 129 words per plane so the // plane scatter hits 32 different banks (the unpadded layout puts all 32 lanes of a group in one bank); @@ -354,9 +354,6 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t const uint32_t chunks = gridDim.x; if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { - if (blockIdx.x == 0 && t == 0) - printf("[fzgmod] pfpl::decompress: not a PFPL stream of %u floats " - "(magic %08x, n %u, chunks %u, archive %u B)\n", n, w[0], w[3], w[5], w[1]); streamTrap(); } const uint32_t cid = blockIdx.x; diff --git a/tests/stages/CMakeLists.txt b/tests/stages/CMakeLists.txt index 7357a5bb..4033dc51 100644 --- a/tests/stages/CMakeLists.txt +++ b/tests/stages/CMakeLists.txt @@ -84,6 +84,22 @@ fz_add_test( LABELS stages gpu ) +# ─── Two-array PFPL against the staged PFPL pipeline ───────────────────────── +fz_add_test( + NAME test_pfpl_stream + SOURCES test_pfpl_stream.cpp + LIBS fzgmod_modules fzgmod_mem fzgmod + LABELS stages gpu +) +# HIP defaults to AMD kernels: also run the generic harness kernels, which are the ones CUDA uses. +if(FZGMOD_BACKEND STREQUAL "HIP") + add_test(NAME test_pfpl_stream_harness COMMAND test_pfpl_stream) + set_tests_properties(test_pfpl_stream_harness PROPERTIES + ENVIRONMENT "HIP_VISIBLE_DEVICES=0;FZ_PFPL_HARNESS=1" + LABELS "stages;gpu" + ) +endif() + # ─── GolombRiceStage GPU round-trip ────────────────────────────────────────── fz_add_test( NAME test_golomb_rice diff --git a/tests/stages/test_pfpl_stream.cpp b/tests/stages/test_pfpl_stream.cpp new file mode 100644 index 00000000..e793fe80 --- /dev/null +++ b/tests/stages/test_pfpl_stream.cpp @@ -0,0 +1,363 @@ +/** + * tests/stages/test_pfpl_stream.cpp + * + * Two-array PFPL (fz::pfpl) against the staged PFPL pipeline in ABS mode. + * + * PS1 PfplStreamCase/MatchesStagedPipeline — archive bytes and decoded values equal the staged + * pipeline's, within the bound, non-finite values restored + * PS2 PfplStream/OverflowStaysInCapacity — a short capacity reports the full length and writes + * nothing past itself + * PS3 PfplStream/ConcurrentStreamsMatchSerial — streams run together produce the serial streams + * PS4 PfplStream/RepeatedCallsOnOnePlan — look-back epochs and tickets survive many calls + * PS5 PfplStream/StreamBytesReadsTheHeader — streamBytes() on valid, foreign and truncated data + * PS6 PfplStream/FreshPlanIgnoresStaleMemory — a new Plan on dirty pages produces the same stream + * + * On HIP, CTest runs this twice: with the AMD kernels (the default) and with the generic harness + * kernels (FZ_PFPL_HARNESS=1), which are the kernels CUDA uses. + */ + +#include +#include "helpers/fz_test_utils.h" +#include "coders/rze/rze_stage.h" +#include "fused/chunk_fusion/pfpl_stream.h" +#include "pipeline/compressor.h" +#include "predictors/diff/diff.h" +#include "quantizers/quantizer/quantizer.h" +#include "shufflers/bitshuffle/bitshuffle_stage.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using fz_test::CudaBuffer; + +namespace { + +constexpr size_t kChunkBytes = 16384; + +// The PFPL preset in ABS mode (examples/presets/pfpl_specialization.toml), staged. +std::unique_ptr stagedPfpl(size_t rawBytes, double eb) { + auto p = std::make_unique(rawBytes, fz::MemoryStrategy::PREALLOCATE); + p->setSpecializationPolicy(fz::SpecializationPolicy::Off); + auto* quant = p->addStage>(); + quant->setErrorBound(eb); + quant->setErrorBoundMode(fz::ErrorBoundMode::ABS); + quant->setQuantRadius(fz::pfpl::kQuantRadius); + quant->setZigzagCodes(true); + quant->setInplaceOutliers(true); + quant->setVerifyReconstruction(true); + auto* diff = p->addStage>(); + diff->setChunkSize(kChunkBytes); + auto* shuffle = p->addStage(); + shuffle->setElementWidth(4); + shuffle->setBlockSize(kChunkBytes); + auto* rze = p->addStage(); + rze->setWordSize(1); + rze->setChunkSize(kChunkBytes); + p->connect(diff, quant, "codes"); + p->connect(shuffle, diff); + p->connect(rze, shuffle); + p->finalize(); + return p; +} + +std::vector makeField(const std::string& kind, size_t n, uint64_t seed) { + std::mt19937_64 rng(seed); + std::normal_distribution g(0.0, 1.0); + std::uniform_real_distribution u(0.0, 1.0); + std::vector f(n); + if (kind == "polar") { + // log2 magnitude and phase of a random state vector, with 20% exact-zero sentinels. + const size_t half = n / 2; + const double scale = 1.0 / std::sqrt(2.0 * static_cast(std::max(half, 1))); + for (size_t i = 0; i < half; ++i) { + if (u(rng) < 0.2) continue; + const double re = g(rng) * scale, im = g(rng) * scale; + f[i] = static_cast(std::log2(std::hypot(re, im)) - 1.0); + f[half + i] = static_cast(std::atan2(im, re) / std::log(2.0)); + } + if (n % 2) f[n - 1] = 0.5f; + } else if (kind == "smooth") { + for (size_t i = 0; i < n; ++i) + f[i] = static_cast(std::sin(i * 1e-3) + 0.5 * std::cos(i * 7e-5)); + } else if (kind == "noise") { + for (size_t i = 0; i < n; ++i) f[i] = static_cast(2.0 * u(rng) - 1.0); + } else if (kind == "spiky") { + // Outliers past the quantization radius, and non-finite values. + const float inf = std::numeric_limits::infinity(); + for (size_t i = 0; i < n; ++i) { + const double r = u(rng); + f[i] = r < 0.001 ? static_cast(1e9 * g(rng)) + : r < 0.0015 ? std::numeric_limits::quiet_NaN() + : r < 0.002 ? inf + : r < 0.0025 ? -inf + : static_cast(1e-3 * g(rng)); + } + } // "zeros" stays all zero + return f; +} + +bool sameBits(float a, float b) { + uint32_t x, y; + std::memcpy(&x, &a, 4); + std::memcpy(&y, &b, 4); + return x == y; +} + +// A stream that does not block host threads, as a concurrent caller would use it. +struct Stream { + cudaStream_t s = nullptr; + Stream() { cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking); } + ~Stream() { cudaStreamDestroy(s); } + Stream(const Stream&) = delete; + Stream& operator=(const Stream&) = delete; + operator cudaStream_t() const { return s; } +}; + +size_t readDeviceSize(const CudaBuffer& d, cudaStream_t s) { + return d.download(s)[0]; +} + +void expectMatchesStaged(const std::string& kind, size_t n, double eb, cudaStream_t s, fz::pfpl::Plan& plan) { + SCOPED_TRACE(kind + " n=" + std::to_string(n) + " eb=" + std::to_string(eb)); + const size_t raw = n * sizeof(float); + const size_t padded = (raw + kChunkBytes - 1) / kChunkBytes * kChunkBytes; + const uint64_t seed = n * 31 + static_cast(eb * 1e9); + std::vector h = makeField(kind, n, seed); + std::vector hPadded = h; + hPadded.resize(padded / sizeof(float), 0.0f); + + CudaBuffer d_in(hPadded.size()), d_out(hPadded.size()); + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + CudaBuffer d_len(1); + d_in.upload(hPadded, s); + + fz::pfpl::compress(plan, d_in.get(), n, static_cast(eb), d_stream.get(), d_stream.bytes(), d_len.get(), s); + const size_t len = readDeviceSize(d_len, s); + ASSERT_LE(len, d_stream.bytes()); + ASSERT_GE(len, sizeof(fz::pfpl::StreamHeader)); + const std::vector stream = d_stream.download_bytes(len, s); + EXPECT_EQ(fz::pfpl::streamBytes(stream.data(), len), len); + fz::pfpl::StreamHeader hdr; + std::memcpy(&hdr, stream.data(), sizeof(hdr)); + EXPECT_EQ(hdr.n, n); + EXPECT_EQ(hdr.abs_eb, static_cast(eb)); + + // The staged pipeline codes the zero-padded input. + auto pipe = stagedPfpl(padded, eb); + const size_t cap = pipe->getMaxCompressedSize(padded); + CudaBuffer d_ref(cap); + CudaBuffer d_refOut(hPadded.size()); + size_t refLen = 0; + pipe->compress(d_in.get(), raw, d_ref.get(), cap, &refLen, s); + cudaStreamSynchronize(s); + const std::vector ref = d_ref.download_bytes(refLen, s); + + const size_t arc = len - sizeof(fz::pfpl::StreamHeader); + ASSERT_EQ(arc, refLen) << "archive length differs from the staged pipeline's"; + EXPECT_EQ(std::memcmp(ref.data(), stream.data() + sizeof(fz::pfpl::StreamHeader), arc), 0) + << "archive bytes differ from the staged pipeline's"; + + size_t refDecoded = 0; + pipe->decompressInto(d_ref.get(), refLen, d_refOut.get(), padded, &refDecoded, s); + fz::pfpl::decompress(plan, d_stream.get(), d_out.get(), n, s); + cudaStreamSynchronize(s); + const std::vector a = d_out.download_bytes(raw, s); + const std::vector b = d_refOut.download_bytes(raw, s); + + size_t differ = 0, overBound = 0, nonFinite = 0; + double maxErr = 0.0; + for (size_t i = 0; i < n; ++i) { + if (!sameBits(a[i], b[i])) ++differ; + if (!std::isfinite(h[i])) { + if (!sameBits(a[i], h[i])) ++nonFinite; + continue; + } + const double e = std::fabs(static_cast(a[i]) - static_cast(h[i])); + maxErr = std::max(maxErr, e); + if (e > eb) ++overBound; + } + EXPECT_EQ(differ, 0u) << "values differ from the staged decode"; + EXPECT_EQ(overBound, 0u) << "values past the bound; largest error " << maxErr; + EXPECT_EQ(nonFinite, 0u) << "non-finite values not restored"; +} + +using CaseParam = std::tuple; +class PfplStreamCase : public ::testing::TestWithParam {}; + +} // namespace + +// PS1 +TEST_P(PfplStreamCase, MatchesStagedPipeline) { + Stream s; + fz::pfpl::Plan plan; + expectMatchesStaged(std::get<0>(GetParam()), std::get<1>(GetParam()), std::get<2>(GetParam()), s, plan); +} + +INSTANTIATE_TEST_SUITE_P( + Fields, PfplStreamCase, + ::testing::Combine(::testing::Values("polar", "smooth", "noise", "zeros", "spiky"), + // under a chunk, one chunk, not a whole number of chunks, and many chunks + ::testing::Values(size_t{512}, size_t{4096}, size_t{3 * 4096 + 2}, size_t{1} << 18, + size_t{1} << 23), + ::testing::Values(1e-2, 1e-4, 1e-6)), + [](const ::testing::TestParamInfo& info) { + return std::string(std::get<0>(info.param)) + "_n" + std::to_string(std::get<1>(info.param)) + "_eb" + + std::to_string(static_cast(std::lround(-std::log10(std::get<2>(info.param))))); + }); + +// PS2 +TEST(PfplStream, OverflowStaysInCapacity) { + Stream s; + fz::pfpl::Plan plan; + const size_t n = 3 * 4096 + 77; + const std::vector h = makeField("noise", n, 5); + CudaBuffer d_in(n); + CudaBuffer d_len(1); + d_in.upload(h, s); + + size_t full = 0; + { + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + fz::pfpl::compress(plan, d_in.get(), n, 1e-7f, d_stream.get(), d_stream.bytes(), d_len.get(), s); + full = readDeviceSize(d_len, s); + } + ASSERT_GT(full, size_t{100}); + + for (const size_t cap : {full - 1, full / 2, size_t{100}, size_t{24}, size_t{8}}) { + SCOPED_TRACE("capacity " + std::to_string(cap)); + const size_t guard = 4096; + CudaBuffer d_stream(cap + guard); + ASSERT_EQ(cudaMemset(d_stream.get(), 0xA5, cap + guard), cudaSuccess); + fz::pfpl::compress(plan, d_in.get(), n, 1e-7f, d_stream.get(), cap, d_len.get(), s); + EXPECT_EQ(readDeviceSize(d_len, s), full); + const std::vector all = d_stream.download(s); + const size_t touched = static_cast( + std::count_if(all.begin() + static_cast(cap), all.end(), [](uint8_t b) { return b != 0xA5; })); + EXPECT_EQ(touched, 0u) << "bytes written past the capacity"; + } +} + +// PS3 +TEST(PfplStream, ConcurrentStreamsMatchSerial) { + constexpr int kStreams = 4, kRounds = 20; + constexpr size_t n = size_t{1} << 22; + constexpr float eb = 1e-4f; + + std::vector> fields; + std::vector> streams; + std::vector> plans; + std::vector>> in, out; + std::vector>> st; + std::vector>> len; + std::vector> serial(kStreams); + for (int i = 0; i < kStreams; ++i) { + fields.push_back(makeField(i % 2 ? "polar" : "smooth", n, 100 + i)); + streams.emplace_back(new Stream); + plans.emplace_back(new fz::pfpl::Plan); + in.emplace_back(new CudaBuffer(n)); + out.emplace_back(new CudaBuffer(n)); + st.emplace_back(new CudaBuffer(fz::pfpl::maxStreamBytes(n))); + len.emplace_back(new CudaBuffer(1)); + in[i]->upload(fields[i], *streams[i]); + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::compress(*plans[i], in[i]->get(), n, eb, st[i]->get(), st[i]->bytes(), len[i]->get(), *streams[i]); + serial[i] = st[i]->download_bytes(readDeviceSize(*len[i], *streams[i]), *streams[i]); + } + + for (int round = 0; round < kRounds; ++round) { + SCOPED_TRACE("round " + std::to_string(round)); + for (int i = 0; i < kStreams; ++i) { + ASSERT_EQ(cudaMemsetAsync(st[i]->get(), 0, st[i]->bytes(), *streams[i]), cudaSuccess); + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::compress(*plans[i], in[i]->get(), n, eb, st[i]->get(), st[i]->bytes(), len[i]->get(), *streams[i]); + } + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + for (int i = 0; i < kStreams; ++i) { + EXPECT_EQ(st[i]->download_bytes(serial[i].size(), *streams[i]), serial[i]) << "stream " << i; + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::decompress(*plans[i], st[i]->get(), out[i]->get(), n, *streams[i]); + } + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + for (int i = 0; i < kStreams; ++i) { + const std::vector d = out[i]->download(*streams[i]); + size_t overBound = 0; + for (size_t j = 0; j < n; ++j) { + overBound += std::fabs(static_cast(d[j]) - static_cast(fields[i][j])) > eb; + } + EXPECT_EQ(overBound, 0u) << "stream " << i; + } + } +} + +// PS4 +TEST(PfplStream, RepeatedCallsOnOnePlan) { + Stream s; + fz::pfpl::Plan plan; + for (int r = 0; r < 300 && !::testing::Test::HasFatalFailure(); ++r) { + expectMatchesStaged(r % 3 ? "polar" : "smooth", (size_t{1} << 20) + (r % 5) * 1000, 1e-3, s, plan); + } +} + +// PS5 +TEST(PfplStream, StreamBytesReadsTheHeader) { + const fz::pfpl::StreamHeader good{fz::pfpl::kStreamMagic, 1000, 1e-3f, 4096}; + const size_t length = sizeof(good) + 1000; + EXPECT_EQ(fz::pfpl::streamBytes(&good), length); + EXPECT_EQ(fz::pfpl::streamBytes(&good, length), length); + EXPECT_EQ(fz::pfpl::streamBytes(&good, length - 1), 0u) << "claims more than the capacity"; + + fz::pfpl::StreamHeader foreign = good; + foreign.magic ^= 1u; + EXPECT_EQ(fz::pfpl::streamBytes(&foreign), 0u); +} + +// PS6: a Plan's counter and look-back state land on recycled pages, so their contents must not reach its first +// kernel. The pages are dirty, and the default stream is busy: a clear queued there would run after the fill, while +// the first kernel on this non-blocking stream goes ahead of it. A stale counter is a wild chunk id, and stale +// look-back words are a wrong offset. +TEST(PfplStream, FreshPlanIgnoresStaleMemory) { + Stream s; + const size_t n = 8 * 4096; + const std::vector h = makeField("smooth", n, 1); + CudaBuffer d_in(n); + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + CudaBuffer d_busy(size_t{1} << 30); + CudaBuffer d_len(1); + d_in.upload(h, s); + + std::vector expected; + { + fz::pfpl::Plan plan; + fz::pfpl::compress(plan, d_in.get(), n, 1e-4f, d_stream.get(), d_stream.bytes(), d_len.get(), s); + expected = d_stream.download_bytes(readDeviceSize(d_len, s), s); + } + + for (int round = 0; round < 40 && !::testing::Test::HasFailure(); ++round) { + SCOPED_TRACE("round " + std::to_string(round)); + { + std::vector>> dirty; + for (int i = 0; i < 64; ++i) { + dirty.emplace_back(new CudaBuffer(4096)); + ASSERT_EQ(cudaMemsetAsync(dirty.back()->get(), 0xFF, 4096, s), cudaSuccess); + } + ASSERT_EQ(cudaStreamSynchronize(s), cudaSuccess); + } + ASSERT_EQ(cudaMemsetAsync(d_busy.get(), 0, d_busy.bytes(), nullptr), cudaSuccess); + fz::pfpl::Plan plan; + fz::pfpl::compress(plan, d_in.get(), n, 1e-4f, d_stream.get(), d_stream.bytes(), d_len.get(), s); + ASSERT_EQ(readDeviceSize(d_len, s), expected.size()); + EXPECT_EQ(d_stream.download_bytes(expected.size(), s), expected); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + } +} From 6ce75f51f24e195754a50e235ee7b91b4263a818 Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Thu, 5 Nov 2026 16:19:50 +0000 Subject: [PATCH 4/8] Link fzgmod_modules against fzgmod_mem fzgmod_modules calls fz::MemoryPool but never declared the dependency, so libfzgmod_modules.so had no DT_NEEDED entry for libfzgmod_mem.so. A consumer that links the library without calling fzgmod_mem itself failed with undefined references on toolchains that default to --as-needed, such as Ubuntu's GCC. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + CMakeLists.txt | 1 + 2 files changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0eec188b..56baa538 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ Version numbers follow [Semantic Versioning](https://semver.org/). ### Fixed - On HIP, the ABS bound check could accept values that decode up to half an ulp past the bound, because clang fused its multiply and subtract into one FMA. `fz::backend::fmulRounded` keeps the product rounded in the staged quantizer and the chunk harness. +- `fzgmod_modules` now links `fzgmod_mem`, which it calls into, so a consumer of the shared library links on toolchains that default to `--as-needed` (such as Ubuntu's GCC). ### Changed - The chunk harness decode copies payloads with `loadChunkBytes` instead of byte by byte, on CUDA as well. diff --git a/CMakeLists.txt b/CMakeLists.txt index 93112a60..bfe30451 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -320,6 +320,7 @@ target_link_libraries(fzgmod_modules fzgmod_compile_settings ${FZGMOD_GPU_RUNTIME} PRIVATE + fzgmod_mem fzgmod_utils ) add_subdirectory(modules) From eef094a2775b293ad4c6f873437f89f56b7fb430 Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Thu, 5 Nov 2026 12:26:09 -0500 Subject: [PATCH 5/8] adds FP64 for ROCm --- include/backend/fp.h | 10 + modules/fused/chunk_fusion/pfpl_stream.cu | 434 ++++++++++++++---- modules/fused/chunk_fusion/pfpl_stream.h | 54 ++- .../fused/chunk_fusion/pfpl_stream_amd.inl | 197 +++++--- tests/stages/test_pfpl_stream.cpp | 333 ++++++++++++++ 5 files changed, 852 insertions(+), 176 deletions(-) diff --git a/include/backend/fp.h b/include/backend/fp.h index 96f96c8e..e04ddb77 100644 --- a/include/backend/fp.h +++ b/include/backend/fp.h @@ -29,5 +29,15 @@ __device__ __forceinline__ float fmulRounded(float a, float b) { #endif } +/** The same for a double: the product the decoder stores, never fused into the operation that consumes it. */ +__device__ __forceinline__ double fmulRounded(double a, double b) { +#if defined(FZGMOD_BACKEND_HIP) + #pragma clang fp contract(off) + return a * b; +#else + return __dmul_rn(a, b); +#endif +} + } // namespace backend } // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream.cu b/modules/fused/chunk_fusion/pfpl_stream.cu index cf70a5c3..61f68f46 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.cu +++ b/modules/fused/chunk_fusion/pfpl_stream.cu @@ -35,10 +35,64 @@ using fused::chunk::Bitshuffle32; using fused::chunk::RZECoder; constexpr uint32_t kElems = Geom::NELEM; // floats per chunk +constexpr uint32_t kElems64 = kChunkBytes / 8; // doubles per chunk +constexpr uint32_t kNPP64 = kElems64 / 32; // 32-bit words per bit-plane of a double chunk constexpr uint32_t kLcHeader = 8; // LC archive: [u32 coded bytes][u32 chunks] constexpr uint32_t kStreamHeader = sizeof(StreamHeader); constexpr int kWave = fz::backend::kWaveSize; +// The double quantizer's parameters: the float ones' rule, with no threshold (the radius check also turns away +// non-finite values). +struct QuantParams64 { + double ebx2_r; + double ebx2; + double abs_eb; +}; + +template struct Elem; +template <> struct Elem { + using Word = uint32_t; + using SWord = int32_t; + using Params = QuantInplaceZigzagParams; + static constexpr uint32_t kMagic = kStreamMagic; +}; +template <> struct Elem { + using Word = uint64_t; + using SWord = int64_t; + using Params = QuantParams64; + static constexpr uint32_t kMagic = kStreamMagic64; +}; + +__device__ __forceinline__ uint64_t doubleBits(double x) { + union { double d; uint64_t u; } b; + b.d = x; + return b.u; +} + +__device__ __forceinline__ double bitsDouble(uint64_t u) { + union { double d; uint64_t u; } b; + b.u = u; + return b.d; +} + +// x's code: its zigzagged quantization index, or its raw bits when that index would not reconstruct within the bound. +// Raw bits read as raw because a raw value is never subnormal (one that small quantizes to index 0 at any allowed +// bound), so they lie at or above 2 * kQuantRadius64. The range check comes first: converting an index out of +// long long's range, or a non-finite one, is undefined. +__device__ __forceinline__ uint64_t quantCode64(double x, const QuantParams64& p) { + const double r = rint(x * p.ebx2_r); + const bool repr = fabs(r) < static_cast(kQuantRadius64) && + fabs(x - fz::backend::fmulRounded(r, p.ebx2)) <= p.abs_eb; + const long long q = repr ? static_cast(r) : 0ll; + return repr ? (static_cast(q) << 1) ^ static_cast(q >> 63) : doubleBits(x); +} + +__device__ __forceinline__ double valueOfCode64(uint64_t code, double ebx2) { + if ((code >> 1) >= kQuantRadius64) return bitsDouble(code); + const long long q = static_cast((code >> 1) ^ (0ull - (code & 1ull))); + return static_cast(q) * ebx2; +} + // Look-back word: [63:62] kind (0 empty, 1 the chunk's own size, 2 its inclusive // prefix) | [61:32] epoch of the call that wrote it | [31:0] bytes. A word left by // another call reads as empty, so the state is never cleared between calls. @@ -122,29 +176,12 @@ __device__ uint32_t lookbackOffset(uint32_t cid, uint32_t own, unsigned long lon return s_offset; } -// One block per chunk: code it in shared memory, find its offset, write its -// size-table entry and payload into the stream; the last chunk adds the headers -// and the length. Layout: [StreamHeader][LC header][u32 sizes x chunks][payloads]. -__global__ void __launch_bounds__(TPB) -encodeKernel(const float* __restrict__ in, uint32_t n, uint32_t coded_n, QuantInplaceZigzagParams qp, - byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, - unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, - unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { - __shared__ uint32_t s_cid; - __shared__ QuantInplaceZigzagParams s_qp; - if (threadIdx.x == 0) { - s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); - s_qp = qp; - } - __syncthreads(); - const uint32_t cid = s_cid; - - const ChunkEncoded e = - fused::chunk::chunk_encode_at, - DiffNegabinary, Bitshuffle32>( - cid, in, coded_n, reinterpret_cast(&s_qp), ChunkSideCtx{}, n); - const uint32_t offset = lookbackOffset(cid, e.nbytes, state, epoch, nullptr); - +// A coded chunk's size-table entry and payload go into the stream; the last chunk adds the headers and the length. +// Layout: [StreamHeader][LC header][u32 sizes x chunks][payloads]. +__device__ __forceinline__ void storeEncoded(const ChunkEncoded& e, uint32_t cid, uint32_t offset, uint32_t magic, + float abs_eb, uint32_t n, byte* __restrict__ out, + unsigned long long capacity, uint32_t chunks, + size_t* __restrict__ d_stream_bytes) { const unsigned long long table = kStreamHeader + kLcHeader; const unsigned long long data = table + 4ull * chunks; if (threadIdx.x == 0 && table + 4ull * (cid + 1u) <= capacity) @@ -159,28 +196,107 @@ encodeKernel(const float* __restrict__ in, uint32_t n, uint32_t coded_n, QuantIn uint32_t* w = reinterpret_cast(out); if (total <= capacity) { for (unsigned long long b = kStreamHeader + archive; b < total; ++b) out[b] = 0; - w[4] = coded_n * 4u; + w[4] = chunks * static_cast(kChunkBytes); w[5] = chunks; } // The stream header goes in whenever it fits, so a stream too long for the // buffer still records the length it needed. if (kStreamHeader <= capacity) { - w[0] = kStreamMagic; + w[0] = magic; w[1] = static_cast(padded); - w[2] = __float_as_uint(qp.abs_eb); + w[2] = __float_as_uint(abs_eb); w[3] = n; } if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); } } +// One block per chunk: code it in shared memory, find its offset, store it. +__global__ void __launch_bounds__(TPB) +encodeKernel(const float* __restrict__ in, uint32_t n, uint32_t coded_n, QuantInplaceZigzagParams qp, + byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, + unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, + unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { + __shared__ uint32_t s_cid; + __shared__ QuantInplaceZigzagParams s_qp; + if (threadIdx.x == 0) { + s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); + s_qp = qp; + } + __syncthreads(); + const uint32_t cid = s_cid; + + const ChunkEncoded e = + fused::chunk::chunk_encode_at, + DiffNegabinary, Bitshuffle32>( + cid, in, coded_n, reinterpret_cast(&s_qp), ChunkSideCtx{}, n); + const uint32_t offset = lookbackOffset(cid, e.nbytes, state, epoch, nullptr); + storeEncoded(e, cid, offset, kStreamMagic, qp.abs_eb, n, out, capacity, chunks, d_stream_bytes); +} + +// One double chunk in the block's shared memory: codes, difference and negabinary in 64 bits, then the bitshuffle and +// RZE. A 64-plane bitshuffle is two 32-plane ones: each 32-element group's high words give planes 0..31 (MSB first, as +// the float planes) and its low words planes 32..63. Each plane's 64-bit words hold an even group's bits above the next +// odd group's, so in memory the odd group's word comes first (group ^ 1), as in the staged 8-byte bitshuffle. +__device__ __forceinline__ ChunkEncoded encodeChunk64(uint32_t cid, const double* __restrict__ in, uint32_t n, + const QuantParams64& qp) { + __shared__ __align__(16) uint64_t sA[kElems64]; + __shared__ __align__(16) uint64_t sB[kElems64]; + __shared__ __align__(16) byte sTemp[fused::chunk::TEMP_BYTES]; + + const size_t base = static_cast(cid) * kElems64; + for (uint32_t i = threadIdx.x; i < kElems64; i += TPB) { + const size_t e = base + i; + sA[i] = quantCode64(e < n ? in[e] : 0.0, qp); + } + __syncthreads(); + + for (uint32_t i = threadIdx.x; i < kElems64; i += TPB) { + const uint64_t d = i == 0 ? sA[i] : sA[i] - sA[i - 1]; + sB[i] = Negabinary::encode(static_cast(d)); + } + __syncthreads(); + + uint32_t* planes = reinterpret_cast(sA); + const int lane = threadIdx.x & 31; + for (uint32_t i = threadIdx.x; i < kElems64; i += TPB) { + const uint32_t g = (i / 32) ^ 1u; + planes[g + lane * kNPP64] = fused::chunk::butterfly32(static_cast(sB[i] >> 32), lane); + planes[g + (32 + lane) * kNPP64] = fused::chunk::butterfly32(static_cast(sB[i]), lane); + } + __syncthreads(); + + int csize = kChunkBytes; + const bool good = RZECoder::encode(csize, reinterpret_cast(sA), reinterpret_cast(sB), + sTemp, nullptr); + __syncthreads(); + if (good && csize < kChunkBytes) + return ChunkEncoded{reinterpret_cast(sB), static_cast(csize), static_cast(csize)}; + return ChunkEncoded{reinterpret_cast(sA), static_cast(kChunkBytes), + (1u << 31) | static_cast(kChunkBytes)}; +} + +__global__ void __launch_bounds__(TPB) +encodeKernel64(const double* __restrict__ in, uint32_t n, QuantParams64 qp, float abs_eb, + byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, + unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, + unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { + __shared__ uint32_t s_cid; + if (threadIdx.x == 0) s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); + __syncthreads(); + const uint32_t cid = s_cid; + const ChunkEncoded e = encodeChunk64(cid, in, n, qp); + const uint32_t offset = lookbackOffset(cid, e.nbytes, state, epoch, nullptr); + storeEncoded(e, cid, offset, kStreamMagic64, abs_eb, n, out, capacity, chunks, d_stream_bytes); +} + // A field of no elements: the headers alone. -__global__ void emptyKernel(byte* __restrict__ out, unsigned long long capacity, float abs_eb, +__global__ void emptyKernel(byte* __restrict__ out, unsigned long long capacity, uint32_t magic, float abs_eb, size_t* __restrict__ d_stream_bytes) { const unsigned long long total = kStreamHeader + kLcHeader; uint32_t* w = reinterpret_cast(out); if (kStreamHeader <= capacity) { - w[0] = kStreamMagic; + w[0] = magic; w[1] = kLcHeader; w[2] = __float_as_uint(abs_eb); w[3] = 0u; @@ -192,16 +308,21 @@ __global__ void emptyKernel(byte* __restrict__ out, unsigned long long capacity, if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); } +// The headers must describe this field. A stream that does not is a caller's error: stop rather than decode garbage. +__device__ __forceinline__ void checkStream(const uint32_t* w, uint32_t magic, uint32_t n, uint32_t chunks) { + if (w[0] != magic || w[3] != n || w[4] != chunks * static_cast(kChunkBytes) || w[5] != chunks || + w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { + streamTrap(); + } +} + // One block per chunk (gridDim.x == chunks): find the chunk's payload, decode it. __global__ void __launch_bounds__(TPB) decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t n, unsigned long long* __restrict__ state, uint32_t epoch) { const uint32_t* w = reinterpret_cast(stream); const uint32_t chunks = gridDim.x; - if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || - w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { - streamTrap(); - } + checkStream(w, kStreamMagic, n, chunks); const uint32_t cid = blockIdx.x; const uint32_t* table = w + (kStreamHeader + kLcHeader) / 4u; const uint32_t entry = table[cid]; @@ -213,49 +334,140 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t 2.0f * __uint_as_float(w[2]), /*inplace_outliers=*/true, kQuantRadius, out + base); } +// One double chunk from the archive into out[0, cnt): RZE, the inverse bitshuffle, then the inverse negabinary and a +// chunk-local inclusive scan, as chunk_inverse_pfpl_at does for floats: each 512-value slab is scanned warp-wise, with +// 16 warp totals and one carry in the coder scratch. +__device__ __forceinline__ void decodeChunk64(const byte* __restrict__ payload, uint32_t entry, int cnt, double ebx2, + double* __restrict__ out) { + __shared__ __align__(16) uint64_t sA[kElems64]; + __shared__ __align__(16) uint64_t sB[kElems64]; + __shared__ __align__(16) byte sTemp[fused::chunk::TEMP_BYTES]; + + const uint32_t stored = entry & 0x7fffffffu; + const bool raw = (entry & 0x80000000u) != 0; + fused::chunk::loadChunkBytes(reinterpret_cast(sA), payload, stored); + __syncthreads(); + + uint64_t* cur = sA; + uint64_t* alt = sB; + if (!raw) { + int csize = static_cast(stored); + RZECoder::decode(csize, reinterpret_cast(sA), reinterpret_cast(sB), sTemp); + __syncthreads(); + cur = sB; + alt = sA; + } + + const int lane = threadIdx.x & 31; + const int warp = threadIdx.x >> 5; + const uint32_t* planes = reinterpret_cast(cur); + for (uint32_t i = threadIdx.x; i < kElems64; i += TPB) { + const uint32_t g = (i / 32) ^ 1u; + const uint32_t hi = fused::chunk::butterfly32(planes[g + lane * kNPP64], lane); + const uint32_t lo = fused::chunk::butterfly32(planes[g + (32 + lane) * kNPP64], lane); + alt[i] = (static_cast(hi) << 32) | lo; + } + __syncthreads(); + + const int64_t* words = reinterpret_cast(alt); + int64_t* qout = reinterpret_cast(cur); // the bitshuffled bytes are done with + int64_t* scan = reinterpret_cast(sTemp); + if (threadIdx.x == 0) scan[16] = 0; + __syncthreads(); + for (int slab = 0; slab < cnt; slab += TPB) { + const int idx = slab + static_cast(threadIdx.x); + const bool valid = idx < cnt; + int64_t v = valid ? Negabinary::decode(static_cast(words[idx])) : 0; + for (int delta = 1; delta < 32; delta <<= 1) { + const int64_t up = fz::backend::shflUp(v, delta, 32); + if (lane >= delta) v += up; + } + if (lane == 31) scan[warp] = v; + __syncthreads(); + if (warp == 0) { + int64_t t = lane < 16 ? scan[lane] : 0; + for (int delta = 1; delta < 32; delta <<= 1) { + const int64_t up = fz::backend::shflUp(t, delta, 32); + if (lane >= delta) t += up; + } + if (lane < 16) scan[lane] = t; + } + __syncthreads(); + const int64_t q = v + (warp ? scan[warp - 1] : 0) + scan[16]; + if (valid) qout[idx] = q; + __syncthreads(); + if (threadIdx.x == 0) scan[16] += scan[15]; + __syncthreads(); + } + for (int i = threadIdx.x; i < cnt; i += TPB) + out[i] = valueOfCode64(static_cast(qout[i]), ebx2); +} + +__global__ void __launch_bounds__(TPB) +decodeKernel64(const byte* __restrict__ stream, double* __restrict__ out, uint32_t n, + unsigned long long* __restrict__ state, uint32_t epoch) { + const uint32_t* w = reinterpret_cast(stream); + const uint32_t chunks = gridDim.x; + checkStream(w, kStreamMagic64, n, chunks); + const uint32_t cid = blockIdx.x; + const uint32_t* table = w + (kStreamHeader + kLcHeader) / 4u; + const uint32_t entry = table[cid]; + const uint32_t offset = lookbackOffset(cid, entry & 0x7fffffffu, state, epoch, table); + const byte* payload = stream + kStreamHeader + kLcHeader + 4ull * chunks + offset; + const uint32_t base = cid * kElems64; + decodeChunk64(payload, entry, static_cast(min(kElems64, n - base)), + 2.0 * static_cast(__uint_as_float(w[2])), out + base); +} + #if defined(FZGMOD_BACKEND_HIP) #include "fused/chunk_fusion/pfpl_stream_amd.inl" -#endif // FZ_PFPL_HARNESS=1 runs the generic harness kernels instead of the AMD ones (same stream; for comparison). bool useHarness() { -#if defined(FZGMOD_BACKEND_HIP) static const bool harness = [] { const char* e = std::getenv("FZ_PFPL_HARNESS"); return e != nullptr && e[0] == '1'; }(); return harness; -#else - return true; -#endif } +#endif -uint32_t chunksOf(size_t n) { return static_cast((n + kElems - 1) / kElems); } +uint32_t chunksOfBytes(size_t bytes) { return static_cast((bytes + kChunkBytes - 1) / kChunkBytes); } + +template +uint32_t chunksOf(size_t n) { return chunksOfBytes(n * sizeof(Real)); } // The archive's size fields are 32-bit: even an all-raw stream must stay below 4 GiB. +template bool fitsArchive(size_t n) { - return kStreamHeader + kLcHeader + static_cast(chunksOf(n)) * (4u + kChunkBytes) + return n < (static_cast(1) << 32) / sizeof(Real) && + kStreamHeader + kLcHeader + static_cast(chunksOf(n)) * (4u + kChunkBytes) < (static_cast(1) << 32); } -void requireAligned(const void* p, const char* what) { - if (reinterpret_cast(p) % 4u != 0) - throw std::invalid_argument(std::string("fz::pfpl: ") + what + " must be 4-byte aligned"); +void requireAligned(const void* p, size_t alignment, const char* what) { + if (reinterpret_cast(p) % alignment != 0) + throw std::invalid_argument(std::string("fz::pfpl: ") + what + " must be " + std::to_string(alignment) + + "-byte aligned"); } } // namespace +template size_t maxStreamBytes(size_t n) { - const size_t chunks = chunksOf(n); - return kStreamHeader + kLcHeader + chunks * (4u + static_cast(kChunkBytes)); + return kStreamHeader + kLcHeader + chunksOf(n) * (4u + static_cast(kChunkBytes)); } +template size_t maxStreamBytes(size_t); +template size_t maxStreamBytes(size_t); size_t streamBytes(const void* h_stream, size_t capacity) { if (h_stream == nullptr || capacity < kStreamHeader) return 0; StreamHeader h; std::memcpy(&h, h_stream, sizeof(h)); - if (h.magic != kStreamMagic || (h.archive_bytes & 3u) != 0 || - h.archive_bytes < kLcHeader + 4ull * chunksOf(h.n)) + if ((h.magic != kStreamMagic && h.magic != kStreamMagic64) || (h.archive_bytes & 3u) != 0) + return 0; + const uint32_t chunks = h.magic == kStreamMagic ? chunksOf(h.n) : chunksOf(h.n); + if (h.archive_bytes < kLcHeader + 4ull * chunks) return 0; const size_t total = kStreamHeader + static_cast(h.archive_bytes); return total <= capacity ? total : 0; @@ -267,8 +479,10 @@ Plan::~Plan() { if (ticket_ != nullptr) (void)cudaFree(ticket_); } -void Plan::reserve(size_t n) { - const size_t chunks = std::max(chunksOf(n), 1); +void Plan::reserve(size_t n) { reserveBytes(n * sizeof(float)); } + +void Plan::reserveBytes(size_t fieldBytes) { + const size_t chunks = std::max(chunksOfBytes(fieldBytes), 1); if (ticket_ == nullptr) { FZ_CUDA_CHECK(cudaMalloc(&ticket_, sizeof(unsigned long long))); clear_ticket_ = true; @@ -305,17 +519,72 @@ uint32_t Plan::nextEpoch(fz::stream_t stream) { return epoch_; } +struct PlanAccess { + template + static void encode(Plan& plan, const Real* d_in, size_t n, const typename Elem::Params& qp, float abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { + requireAligned(d_in, sizeof(Real), "the input"); + requireAligned(d_out, 4, "the stream buffer"); + plan.reserveBytes(n * sizeof(Real)); + const uint32_t epoch = plan.nextEpoch(stream); + const uint32_t chunks = chunksOf(n); + byte* out = static_cast(d_out); + if (chunks == 0) { + emptyKernel<<<1, 1, 0, stream>>>(out, capacity, Elem::kMagic, abs_eb, d_stream_bytes); + FZ_CUDA_CHECK(cudaGetLastError()); + return; + } +#if defined(FZGMOD_BACKEND_HIP) + if (!useHarness()) + amd::encodeKernel<<>>(d_in, static_cast(n), qp, abs_eb, out, + capacity, chunks, plan.state_, plan.ticket_, + plan.tickets_, epoch, d_stream_bytes); + else +#endif + if constexpr (sizeof(Real) == 4) + encodeKernel<<>>(d_in, static_cast(n), chunks * kElems, qp, out, capacity, + chunks, plan.state_, plan.ticket_, plan.tickets_, epoch, + d_stream_bytes); + else + encodeKernel64<<>>(d_in, static_cast(n), qp, abs_eb, out, capacity, chunks, + plan.state_, plan.ticket_, plan.tickets_, epoch, d_stream_bytes); + FZ_CUDA_CHECK(cudaGetLastError()); + plan.tickets_ += chunks; + } + + template + static void decode(Plan& plan, const void* d_stream, Real* d_out, size_t n, fz::stream_t stream) { + if (n == 0) return; + if (!fitsArchive(n)) + throw std::invalid_argument("fz::pfpl::decompress: no stream holds a field of 4 GiB or more"); + requireAligned(d_stream, 4, "the stream"); + requireAligned(d_out, sizeof(Real), "the output"); + plan.reserveBytes(n * sizeof(Real)); + const uint32_t epoch = plan.nextEpoch(stream); + const uint32_t chunks = chunksOf(n); + const byte* in = static_cast(d_stream); +#if defined(FZGMOD_BACKEND_HIP) + if (!useHarness()) + amd::decodeKernel<<>>(in, d_out, static_cast(n), plan.state_, + epoch); + else +#endif + if constexpr (sizeof(Real) == 4) + decodeKernel<<>>(in, d_out, static_cast(n), plan.state_, epoch); + else + decodeKernel64<<>>(in, d_out, static_cast(n), plan.state_, epoch); + FZ_CUDA_CHECK(cudaGetLastError()); + } +}; + void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { if (!(abs_eb > 0.0f) || !std::isfinite(abs_eb) || !std::isfinite(0.5f / abs_eb)) throw std::invalid_argument("fz::pfpl::compress: the error bound must be positive and finite, " "with a finite quantizer scale 0.5/abs_eb"); - if (!fitsArchive(n)) + if (!fitsArchive(n)) throw std::invalid_argument("fz::pfpl::compress: a field of 4 GiB or more (in whole chunks) exceeds " "the LC archive's 32-bit sizes"); - requireAligned(d_out, "the stream buffer"); - plan.reserve(n); - const uint32_t epoch = plan.nextEpoch(stream); // examples/presets/pfpl_specialization.toml in ABS mode, as QuantizerStage::getFusedOp() packs it. QuantInplaceZigzagParams qp{}; @@ -325,45 +594,36 @@ void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, qp.radius = kQuantRadius; qp.threshold = std::numeric_limits::infinity(); qp.verify_reconstruction = 1; + PlanAccess::encode(plan, d_in, n, qp, abs_eb, d_out, capacity, d_stream_bytes, stream); +} - const uint32_t chunks = chunksOf(n); - if (chunks == 0) { - emptyKernel<<<1, 1, 0, stream>>>(static_cast(d_out), capacity, abs_eb, d_stream_bytes); - FZ_CUDA_CHECK(cudaGetLastError()); - return; - } -#if defined(FZGMOD_BACKEND_HIP) - if (!useHarness()) - amd::encodeKernel<<>>(d_in, static_cast(n), qp, - static_cast(d_out), capacity, chunks, - plan.state_, plan.ticket_, plan.tickets_, epoch, - d_stream_bytes); - else -#endif - encodeKernel<<>>(d_in, static_cast(n), chunks * kElems, qp, - static_cast(d_out), capacity, chunks, plan.state_, - plan.ticket_, plan.tickets_, epoch, d_stream_bytes); - FZ_CUDA_CHECK(cudaGetLastError()); - plan.tickets_ += chunks; +void compress(Plan& plan, const double* d_in, size_t n, double abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { + if (!(abs_eb > 0.0) || !std::isfinite(abs_eb)) + throw std::invalid_argument("fz::pfpl::compress: the error bound must be positive and finite"); + // The bound the stream records: never above the one asked for. + float recorded = static_cast(abs_eb); + if (static_cast(recorded) > abs_eb) recorded = std::nextafter(recorded, 0.0f); + if (recorded < std::numeric_limits::min()) + throw std::invalid_argument("fz::pfpl::compress: a double stream records its bound as a float, so the bound " + "must be at least the smallest normal float"); + if (!fitsArchive(n)) + throw std::invalid_argument("fz::pfpl::compress: a field of 4 GiB or more (in whole chunks) exceeds " + "the LC archive's 32-bit sizes"); + const double eb = static_cast(recorded); + QuantParams64 qp{}; + qp.ebx2_r = 0.5 / eb; + qp.ebx2 = 2.0 * eb; + qp.abs_eb = eb; + PlanAccess::encode(plan, d_in, n, qp, recorded, d_out, capacity, d_stream_bytes, stream); } void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream) { - if (n == 0) return; - if (!fitsArchive(n)) - throw std::invalid_argument("fz::pfpl::decompress: no stream holds a field of 4 GiB or more"); - requireAligned(d_stream, "the stream"); - requireAligned(d_out, "the output"); - plan.reserve(n); - const uint32_t epoch = plan.nextEpoch(stream); -#if defined(FZGMOD_BACKEND_HIP) - if (!useHarness()) - amd::decodeKernel<<>>(static_cast(d_stream), d_out, - static_cast(n), plan.state_, epoch); - else -#endif - decodeKernel<<>>(static_cast(d_stream), d_out, - static_cast(n), plan.state_, epoch); - FZ_CUDA_CHECK(cudaGetLastError()); + PlanAccess::decode(plan, d_stream, d_out, n, stream); +} + +void decompress(Plan& plan, const void* d_stream, double* d_out, size_t n, fz::stream_t stream) { + PlanAccess::decode(plan, d_stream, d_out, n, stream); } } // namespace pfpl diff --git a/modules/fused/chunk_fusion/pfpl_stream.h b/modules/fused/chunk_fusion/pfpl_stream.h index 7d62067e..79029525 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.h +++ b/modules/fused/chunk_fusion/pfpl_stream.h @@ -2,9 +2,9 @@ /** * @file modules/fused/chunk_fusion/pfpl_stream.h - * @brief Two-array PFPL: a float field compresses straight into one caller - * buffer and decompresses straight back, stream-ordered, with no host - * synchronization and no device temporaries. + * @brief Two-array PFPL: a float or double field compresses straight into one + * caller buffer and decompresses straight back, stream-ordered, with no + * host synchronization and no device temporaries. * * The stream is self-describing: a 16-byte StreamHeader, then the LC/RZE * archive that the PFPL preset (examples/presets/pfpl_specialization.toml with @@ -16,9 +16,17 @@ * decoupled look-back, so neither direction touches device memory beyond the * two arrays and a Plan's 8 bytes of look-back state per chunk. * + * A double field takes the same steps in 64 bits: a 64-bit zigzag code, the + * difference and negabinary in 64 bits, a 64-plane bitshuffle, and RZE over the + * same 16 KiB chunks. The preset's quantizer has no double form with in-place + * outliers, so there is no staged twin of the whole chain, but the difference, + * bitshuffle and RZE match the staged stages byte for byte. kStreamMagic64 marks + * these streams, and they are read back only as doubles. + * * Archive bytes are a multiple of 4; a field whose length is not a whole number * of chunks is coded padded with zeros (as the staged pipeline pads it) and - * decodes to exactly its own length. Both arrays must be 4-byte aligned. + * decodes to exactly its own length. The stream must be 4-byte aligned and the + * field aligned for its type. */ #include "backend/types.h" @@ -34,15 +42,18 @@ struct StreamHeader { uint32_t magic; ///< kStreamMagic uint32_t archive_bytes; ///< LC archive bytes after the header (a multiple of 4) float abs_eb; ///< absolute error bound every value was coded to - uint32_t n; ///< float elements of the field + uint32_t n; ///< elements of the field }; static_assert(sizeof(StreamHeader) == 16, "the PFPL stream header is 16 bytes"); -constexpr uint32_t kStreamMagic = 0x4C504650u; ///< "PFPL" in little-endian byte order -constexpr int kChunkBytes = 16384; ///< the preset's Difference/Bitshuffle/RZE chunk -constexpr uint32_t kQuantRadius = 4194304; ///< the preset's quant_radius +constexpr uint32_t kStreamMagic = 0x4C504650u; ///< "PFPL" in little-endian byte order: a float stream +constexpr uint32_t kStreamMagic64 = 0x44504650u; ///< "PFPD": a double stream +constexpr int kChunkBytes = 16384; ///< the preset's Difference/Bitshuffle/RZE chunk +constexpr uint32_t kQuantRadius = 4194304; ///< the preset's quant_radius: 2^(float's mantissa bits - 1) +constexpr uint64_t kQuantRadius64 = 1ull << 51; ///< the same for a double: raw bits of a normal value lie above it -/** Longest stream n floats can produce (every chunk stored raw). */ +/** Longest stream n elements of T (float or double) can produce (every chunk stored raw). */ +template size_t maxStreamBytes(size_t n); /** Length of the stream whose first bytes are at `h_stream` (host memory), or 0 @@ -54,6 +65,8 @@ size_t streamBytes(const void* h_stream, size_t capacity = SIZE_MAX); * stream, or synchronized between streams): use one Plan per concurrently used * stream. It grows to the largest field it has seen. */ +struct PlanAccess; + class Plan { public: Plan() = default; @@ -65,9 +78,11 @@ class Plan { * cleared on the stream of the next compress or decompress. */ void reserve(size_t n); + /** The same, for a field of `fieldBytes` of either type: the state depends only on its chunks. */ + void reserveBytes(size_t fieldBytes); + private: - friend void compress(Plan&, const float*, size_t, float, void*, size_t, size_t*, fz::stream_t); - friend void decompress(Plan&, const void*, float*, size_t, fz::stream_t); + friend struct PlanAccess; uint32_t nextEpoch(fz::stream_t stream); unsigned long long* state_ = nullptr; ///< one look-back word per chunk, tagged with its call's epoch @@ -93,12 +108,21 @@ void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream); /** - * Decompress the compress() stream at d_stream, which must describe n floats, - * into d_out[0, n). Stream-ordered, no host synchronization. A stream that is - * not a PFPL stream of n floats stops the kernel (device trap) rather than - * decoding garbage. + * compress() for doubles. The stream records its bound as a float, so abs_eb is + * rounded down to the nearest float and that is the bound the values are coded + * to; it must be at least the smallest normal float (about 1.2e-38). + */ +void compress(Plan& plan, const double* d_in, size_t n, double abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream); + +/** + * Decompress the compress() stream at d_stream, which must describe n elements + * of this type, into d_out[0, n). Stream-ordered, no host synchronization. A + * stream that is not a PFPL stream of n elements of this type stops the kernel + * (device trap) rather than decoding garbage. */ void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream); +void decompress(Plan& plan, const void* d_stream, double* d_out, size_t n, fz::stream_t stream); } // namespace pfpl } // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream_amd.inl b/modules/fused/chunk_fusion/pfpl_stream_amd.inl index ce2b0dfe..0dded014 100644 --- a/modules/fused/chunk_fusion/pfpl_stream_amd.inl +++ b/modules/fused/chunk_fusion/pfpl_stream_amd.inl @@ -13,19 +13,56 @@ // Decoding mirrors it: an element of a repetition-eliminated level is data[k - 1], k the set bitmap bits up to and // including it, so each level decodes in parallel from popcounts; the chunk-local prefix sum runs as eight // wave scans at once with one barrier, in the element order the inverse bitshuffle leaves in registers. +// A double chunk runs the same kernels with 64-bit codes in registers (4 elements per thread, not 8) and twice the +// planes at half the length, each 32-element group's high and low words going through butterfly32 separately (see +// encodeChunk64); from zero elimination on, a chunk is bytes and nothing differs. // Tuned for, and only used on, 64-lane wavefronts. namespace amd { constexpr int kTPB = 512; // 8 wavefronts; must equal the harness's TPB (storeChunkBytes) -constexpr int kPerT = kElems / kTPB; // 8 elements per thread, and 8 bitshuffled words -constexpr int kNPP = kElems / 32; // words per bitshuffle plane -constexpr int kPS = kNPP + 1; // padded plane stride (words) -constexpr int kBufWords = 32 * kPS; +constexpr int kWordsT = kChunkBytes / 4 / kTPB; // 8 bitshuffled words per thread constexpr int kL1 = kChunkBytes / 8; // zero-elimination bitmap bytes constexpr uint32_t kAvail = kChunkBytes - 2; // d_RZE: room before its 2-byte trailer -static_assert(kTPB == TPB && kTPB == 8 * 64 && kPerT == 8, "one 512-thread block per 16 KiB chunk"); -static_assert(kBufWords * 4 >= kChunkBytes + 4, "the packed chunk is compacted in place in the plane buffer"); +static_assert(kTPB == TPB && kTPB == 8 * 64 && kWordsT == 8, "one 512-thread block per 16 KiB chunk"); + +// A chunk's geometry in registers and in the plane buffer. +template +struct Layout { + using Word = typename Elem::Word; + using SWord = typename Elem::SWord; + static constexpr int kPerChunk = kChunkBytes / sizeof(Real); + static constexpr int kPerT = kPerChunk / kTPB; // elements per thread: 8 floats or 4 doubles + static constexpr int kPlanes = 8 * sizeof(Real); + static constexpr int kNPP = kPerChunk / 32; // 32-bit words per bitshuffle plane + static constexpr int kPS = kNPP + 1; // padded plane stride (words) + static constexpr int kBufWords = kPlanes * kPS; + static constexpr int kRowLog2 = sizeof(Real) == 4 ? 4 : 3; // log2 threads per plane + static_assert(kBufWords * 4 >= kChunkBytes + 4, "the packed chunk is compacted in place in the plane buffer"); + static_assert(kPlanes * kNPP * 4 == kChunkBytes && (kNPP / kWordsT) == (1 << kRowLog2), "planes fill the chunk"); +}; + +// Thread t's 8 words of the bitshuffled chunk: bytes 32t..32t+31. +template +__device__ __forceinline__ const uint32_t* planeRow(const uint32_t* buf, int t) { + using L = Layout; + return buf + (t >> L::kRowLog2) * L::kPS + (t & ((1 << L::kRowLog2) - 1)) * kWordsT; +} + +// Quantized code of x (QuantInplaceZigzag, verify_reconstruction): zigzag code, raw bits for outliers. +__device__ __forceinline__ uint32_t quantCode(float x, const QuantInplaceZigzagParams& qp) { + const int q = static_cast(roundf(x * qp.ebx2_r)); + bool repr = q > -static_cast(qp.radius) && q < static_cast(qp.radius) && fabsf(x) < qp.threshold; + if (repr) repr = fabsf(x - fz::backend::fmulRounded(static_cast(q), qp.ebx2)) <= qp.abs_eb; + return repr ? static_cast((q << 1) ^ (q >> 31)) : __float_as_uint(x); +} +__device__ __forceinline__ uint64_t quantCode(double x, const QuantParams64& qp) { return quantCode64(x, qp); } + +__device__ __forceinline__ float valueOfCode(uint32_t code, float ebx2) { + if ((code >> 1) >= kQuantRadius) return __uint_as_float(code); + return static_cast(static_cast((code >> 1) ^ (0u - (code & 1u)))) * ebx2; +} +__device__ __forceinline__ double valueOfCode(uint64_t code, double ebx2) { return valueOfCode64(code, ebx2); } __device__ __forceinline__ uint32_t inclMask(int j) { return j == 31 ? 0xffffffffu : (2u << j) - 1u; } @@ -90,16 +127,22 @@ __device__ __forceinline__ void store32(byte* __restrict__ dst, const uint32_t v } } +template __global__ void __launch_bounds__(kTPB) -encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams qp, +encodeKernel(const Real* __restrict__ in, uint32_t n, typename Elem::Params qp, float abs_eb, byte* __restrict__ out, unsigned long long capacity, uint32_t chunks, unsigned long long* __restrict__ state, unsigned long long* __restrict__ ticket, unsigned long long ticket_base, uint32_t epoch, size_t* __restrict__ d_stream_bytes) { - __shared__ __align__(16) uint32_t buf[kBufWords]; // bitshuffled planes (padded), then the packed chunk + using L = Layout; + using Word = typename L::Word; + using SWord = typename L::SWord; + constexpr int kPerT = L::kPerT; + constexpr int kPS = L::kPS; + __shared__ __align__(16) uint32_t buf[L::kBufWords]; // bitshuffled planes (padded), then the packed chunk __shared__ __align__(16) uint32_t l1w[kL1 / 4]; // L1 bitmap, then its repetition-eliminated bytes __shared__ uint8_t re2[kL1 / 8]; // L2's repetition-eliminated bytes __shared__ uint8_t re3[32]; // L3's - __shared__ uint32_t scratch[64]; + __shared__ __align__(8) uint32_t scratch[64]; __shared__ uint32_t s_cid; __shared__ uint32_t s_info[8]; // raw, csize, offset, n1, n2, n3, L4, wpos @@ -107,46 +150,48 @@ encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams if (t == 0) s_cid = static_cast(atomicAdd(ticket, 1ull) - ticket_base); __syncthreads(); const uint32_t cid = s_cid; - const size_t base = static_cast(cid) * kElems; + const size_t base = static_cast(cid) * L::kPerChunk; - // Quantize elements i = t + 512k (QuantInplaceZigzag, verify_reconstruction): zigzag codes, raw bits for outliers. - uint32_t c[kPerT]; + // Quantize elements i = t + 512k. + Word c[kPerT]; #pragma unroll for (int k = 0; k < kPerT; ++k) { const size_t i = base + t + kTPB * k; - const float x = i < n ? in[i] : 0.0f; - const int q = static_cast(roundf(x * qp.ebx2_r)); - bool repr = q > -static_cast(qp.radius) && q < static_cast(qp.radius) && fabsf(x) < qp.threshold; - if (repr) repr = fabsf(x - fz::backend::fmulRounded(static_cast(q), qp.ebx2)) <= qp.abs_eb; - c[k] = repr ? static_cast((q << 1) ^ (q >> 31)) : __float_as_uint(x); + c[k] = quantCode(i < n ? in[i] : static_cast(0), qp); } // Difference with element i-1 (DiffNegabinary): lane-1 of this wavefront, or for lane 0 the previous // wavefront's lane 63 (wavefront 0: wavefront 7's of slab k-1; element 0 keeps its code). + Word* lastLane = reinterpret_cast(scratch); if (lane == 63) { #pragma unroll - for (int k = 0; k < kPerT; ++k) scratch[wave * kPerT + k] = c[k]; + for (int k = 0; k < kPerT; ++k) lastLane[wave * kPerT + k] = c[k]; } __syncthreads(); #pragma unroll for (int k = 0; k < kPerT; ++k) { - uint32_t prev = fz::backend::shflUp(c[k], 1, 64); - if (lane == 0) prev = wave > 0 ? scratch[(wave - 1) * kPerT + k] : (k > 0 ? scratch[7 * kPerT + k - 1] : 0u); - const int32_t d = (t == 0 && k == 0) ? static_cast(c[k]) : static_cast(c[k] - prev); - // Bitshuffle32: the butterfly leaves plane (t & 31) of group (t >> 5) + 16k in this lane. - buf[(t & 31) * kPS + (t >> 5) + 16 * k] = - fused::chunk::butterfly32(Negabinary::encode(d), t & 31); + Word prev = fz::backend::shflUp(c[k], 1, 64); + if (lane == 0) prev = wave > 0 ? lastLane[(wave - 1) * kPerT + k] : (k > 0 ? lastLane[7 * kPerT + k - 1] : Word(0)); + const Word d = (t == 0 && k == 0) ? c[k] : c[k] - prev; + const Word nb = Negabinary::encode(static_cast(d)); + // The butterfly leaves plane (t & 31) of group (t >> 5) + 16k in this lane. + const int group = (t >> 5) + 16 * k; + if constexpr (sizeof(Real) == 4) { + buf[(t & 31) * kPS + group] = fused::chunk::butterfly32(nb, t & 31); + } else { + buf[(t & 31) * kPS + (group ^ 1)] = fused::chunk::butterfly32(static_cast(nb >> 32), t & 31); + buf[(32 + (t & 31)) * kPS + (group ^ 1)] = fused::chunk::butterfly32(static_cast(nb), t & 31); + } } __syncthreads(); - // Zero elimination (d_ZEencodebyteshort): thread t owns bytes 32t..32t+31 of the bitshuffled chunk, i.e. plane - // t >> 4, words 8(t & 15)..+7. - uint32_t v[kPerT]; + // Zero elimination (d_ZEencodebyteshort): thread t owns bytes 32t..32t+31 of the bitshuffled chunk. + uint32_t v[kWordsT]; uint32_t bmp = 0; { - const uint32_t* row = buf + (t >> 4) * kPS + (t & 15) * 8; + const uint32_t* row = planeRow(buf, t); #pragma unroll - for (int r = 0; r < kPerT; ++r) { + for (int r = 0; r < kWordsT; ++r) { v[r] = row[r]; #pragma unroll for (int b = 0; b < 4; ++b) bmp |= static_cast(((v[r] >> (8 * b)) & 0xffu) != 0u) << (4 * r + b); @@ -163,14 +208,14 @@ encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams raw = 1; } else if (wpos != 0) { // L1 -> RE data + L2: lane owns L1 bytes 32*lane..+31; L2 bit e = (byte e != byte e-1), byte -1 = 0. - uint32_t b1[kPerT]; + uint32_t b1[kWordsT]; #pragma unroll - for (int r = 0; r < kPerT; ++r) b1[r] = l1w[kPerT * lane + r]; + for (int r = 0; r < kWordsT; ++r) b1[r] = l1w[kWordsT * lane + r]; uint32_t prevByte = fz::backend::shflUp(b1[7], 1, 64) >> 24; if (lane == 0) prevByte = 0; uint32_t m1 = 0; #pragma unroll - for (int r = 0; r < kPerT; ++r) + for (int r = 0; r < kWordsT; ++r) #pragma unroll for (int b = 0; b < 4; ++b) { const uint32_t cur = (b1[r] >> (8 * b)) & 0xffu; @@ -186,7 +231,7 @@ encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams byte* l1b = reinterpret_cast(l1w); // compacted in place: positions only move down uint32_t q = i1 - c1; #pragma unroll - for (int r = 0; r < kPerT; ++r) + for (int r = 0; r < kWordsT; ++r) #pragma unroll for (int b = 0; b < 4; ++b) if ((m1 >> (4 * r + b)) & 1u) l1b[q++] = static_cast(b1[r] >> (8 * b)); @@ -253,7 +298,7 @@ encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams byte* bb = reinterpret_cast(buf); uint32_t q = pos; #pragma unroll - for (int r = 0; r < kPerT; ++r) + for (int r = 0; r < kWordsT; ++r) #pragma unroll for (int b = 0; b < 4; ++b) if ((bmp >> (4 * r + b)) & 1u) bb[q++] = static_cast(v[r] >> (8 * b)); @@ -280,13 +325,13 @@ encodeKernel(const float* __restrict__ in, uint32_t n, QuantInplaceZigzagParams uint32_t* w = reinterpret_cast(out); if (total <= capacity) { for (unsigned long long b = kStreamHeader + archive; b < total; ++b) out[b] = 0; - w[4] = static_cast(chunks) * kElems * 4u; + w[4] = chunks * static_cast(kChunkBytes); w[5] = chunks; } if (kStreamHeader <= capacity) { - w[0] = kStreamMagic; + w[0] = Elem::kMagic; w[1] = static_cast(padded); - w[2] = __float_as_uint(qp.abs_eb); + w[2] = __float_as_uint(abs_eb); w[3] = n; } if (d_stream_bytes != nullptr) *d_stream_bytes = static_cast(total); @@ -329,33 +374,36 @@ __device__ __forceinline__ uint32_t decodeBitmaps(const byte* payload, uint32_t const uint32_t c1 = __popc(l2); const uint32_t e1 = waveInclusiveScan(c1) - c1; #pragma unroll - for (int r = 0; r < kPerT; ++r) { + for (int r = 0; r < kWordsT; ++r) { uint32_t word = 0; #pragma unroll for (int b = 0; b < 4; ++b) { const uint32_t k = e1 + __popc(l2 & inclMask(4 * r + b)); word |= (k ? static_cast(payload[rpos + k - 1]) : 0u) << (8 * b); } - l1w[kPerT * lane + r] = word; + l1w[kWordsT * lane + r] = word; } return rpos; } +template __global__ void __launch_bounds__(kTPB) -decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t n, +decodeKernel(const byte* __restrict__ stream, Real* __restrict__ out, uint32_t n, unsigned long long* __restrict__ state, uint32_t epoch) { - __shared__ __align__(16) uint32_t buf[kBufWords]; // the packed chunk, then the bitshuffled planes (padded) + using L = Layout; + using Word = typename L::Word; + using SWord = typename L::SWord; + constexpr int kPerT = L::kPerT; + constexpr int kPS = L::kPS; + __shared__ __align__(16) uint32_t buf[L::kBufWords]; // the packed chunk, then the bitshuffled planes (padded) __shared__ __align__(16) uint32_t l1w[kL1 / 4]; - __shared__ uint32_t scratch[64]; + __shared__ __align__(8) uint32_t scratch[64]; __shared__ uint32_t s_offset; const int t = threadIdx.x, lane = t & 63, wave = t >> 6; const uint32_t* w = reinterpret_cast(stream); const uint32_t chunks = gridDim.x; - if (w[0] != kStreamMagic || w[3] != n || w[4] != chunks * kElems * 4u || w[5] != chunks || - w[1] < kLcHeader + 4u * chunks || !(__uint_as_float(w[2]) > 0.0f)) { - streamTrap(); - } + checkStream(w, Elem::kMagic, n, chunks); const uint32_t cid = blockIdx.x; const uint32_t* table = w + (kStreamHeader + kLcHeader) / 4u; const uint32_t entry = table[cid]; @@ -368,7 +416,7 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t __syncthreads(); const byte* payload = stream + kStreamHeader + kLcHeader + 4ull * chunks + s_offset; - uint32_t v[kPerT]; + uint32_t v[kWordsT]; if (raw) { load32(payload + 32u * t, v); } else { @@ -381,7 +429,7 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t uint32_t q = blockExclusiveScan(__popc(bmp), scratch, total); const byte* bb = reinterpret_cast(buf); #pragma unroll - for (int r = 0; r < kPerT; ++r) { + for (int r = 0; r < kWordsT; ++r) { uint32_t word = 0; #pragma unroll for (int b = 0; b < 4; ++b) @@ -391,56 +439,57 @@ decodeKernel(const byte* __restrict__ stream, float* __restrict__ out, uint32_t __syncthreads(); // every gather done before the planes overwrite the packed bytes } { - uint32_t* row = buf + (t >> 4) * kPS + (t & 15) * 8; + uint32_t* row = const_cast(planeRow(buf, t)); #pragma unroll - for (int r = 0; r < kPerT; ++r) row[r] = v[r]; + for (int r = 0; r < kWordsT; ++r) row[r] = v[r]; } __syncthreads(); // Inverse bitshuffle and negabinary: element i = t + 512k; then its chunk-local inclusive prefix sum. - uint32_t d[kPerT]; + Word d[kPerT]; #pragma unroll - for (int k = 0; k < kPerT; ++k) - d[k] = static_cast(Negabinary::decode( - fused::chunk::butterfly32(buf[(t & 31) * kPS + (t >> 5) + 16 * k], t & 31))); + for (int k = 0; k < kPerT; ++k) { + const int group = (t >> 5) + 16 * k; + Word nb; + if constexpr (sizeof(Real) == 4) { + nb = fused::chunk::butterfly32(buf[(t & 31) * kPS + group], t & 31); + } else { + const uint32_t hi = fused::chunk::butterfly32(buf[(t & 31) * kPS + (group ^ 1)], t & 31); + const uint32_t lo = fused::chunk::butterfly32(buf[(32 + (t & 31)) * kPS + (group ^ 1)], t & 31); + nb = (static_cast(hi) << 32) | lo; + } + d[k] = static_cast(Negabinary::decode(nb)); + } #pragma unroll for (int dd = 1; dd < 64; dd <<= 1) #pragma unroll for (int k = 0; k < kPerT; ++k) { - const uint32_t up = fz::backend::shflUp(d[k], dd, 64); + const Word up = fz::backend::shflUp(d[k], dd, 64); if (lane >= dd) d[k] += up; } + Word* waveTotal = reinterpret_cast(scratch); if (lane == 63) { #pragma unroll - for (int k = 0; k < kPerT; ++k) scratch[k * 8 + wave] = d[k]; + for (int k = 0; k < kPerT; ++k) waveTotal[k * 8 + wave] = d[k]; } __syncthreads(); - const size_t base = static_cast(cid) * kElems; - const uint32_t cnt = min(static_cast(kElems), n - cid * kElems); - const float ebx2 = 2.0f * __uint_as_float(w[2]); - uint32_t carry = 0; + const size_t base = static_cast(cid) * L::kPerChunk; + const uint32_t cnt = min(static_cast(L::kPerChunk), n - cid * static_cast(L::kPerChunk)); + const Real ebx2 = static_cast(2) * static_cast(__uint_as_float(w[2])); + Word carry = 0; #pragma unroll for (int k = 0; k < kPerT; ++k) { - uint32_t before = 0, slab = 0; + Word before = 0, slab = 0; #pragma unroll for (int ww = 0; ww < 8; ++ww) { - const uint32_t x = scratch[k * 8 + ww]; - before += ww < wave ? x : 0u; + const Word x = waveTotal[k * 8 + ww]; + before += ww < wave ? x : Word(0); slab += x; } - const uint32_t code = d[k] + before + carry; + const Word code = d[k] + before + carry; carry += slab; const uint32_t e = t + kTPB * k; - if (e < cnt) { - float y; - if ((code >> 1) >= kQuantRadius) { - y = __uint_as_float(code); - } else { - const int32_t qq = static_cast((code >> 1) ^ (0u - (code & 1u))); - y = static_cast(qq) * ebx2; - } - out[base + e] = y; - } + if (e < cnt) out[base + e] = valueOfCode(code, ebx2); } } diff --git a/tests/stages/test_pfpl_stream.cpp b/tests/stages/test_pfpl_stream.cpp index e793fe80..01a8bbc0 100644 --- a/tests/stages/test_pfpl_stream.cpp +++ b/tests/stages/test_pfpl_stream.cpp @@ -12,6 +12,18 @@ * PS5 PfplStream/StreamBytesReadsTheHeader — streamBytes() on valid, foreign and truncated data * PS6 PfplStream/FreshPlanIgnoresStaleMemory — a new Plan on dirty pages produces the same stream * + * Doubles have no staged quantizer with in-place outliers, so the quantizer's rule is computed here and the rest of + * the chain is compared with the staged stages: + * + * PD1 PfplStream64Case/MatchesStagedStages — archive bytes equal the staged difference, 8-byte bitshuffle and RZE + * of the expected codes; decoded values within the bound, non-finite + * values restored + * PD2 PfplStream64/OverflowStaysInCapacity — as PS2 + * PD3 PfplStream64/ConcurrentStreamsMatchSerial — as PS3 + * PD4 PfplStream64/SharesAPlanWithFloats — float and double fields alternate through one Plan + * PD5 PfplStream64/StreamBytesReadsTheHeader — as PS5 + * PD6 PfplStream64/BoundIsRecordedAsAFloat — never rounded up; too small or not finite is refused + * * On HIP, CTest runs this twice: with the AMD kernels (the default) and with the generic harness * kernels (FZ_PFPL_HARNESS=1), which are the kernels CUDA uses. */ @@ -361,3 +373,324 @@ TEST(PfplStream, FreshPlanIgnoresStaleMemory) { ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); } } + +// ─── Doubles ──────────────────────────────────────────────────────────────────────────────────────────────────────── + +namespace { + +constexpr size_t kElems64 = kChunkBytes / sizeof(double); + +std::vector makeField64(const std::string& kind, size_t n, uint64_t seed) { + std::mt19937_64 rng(seed); + std::normal_distribution g(0.0, 1.0); + std::uniform_real_distribution u(0.0, 1.0); + std::vector f(n); + if (kind == "polar") { + const size_t half = n / 2; + const double scale = 1.0 / std::sqrt(2.0 * static_cast(std::max(half, 1))); + for (size_t i = 0; i < half; ++i) { + if (u(rng) < 0.2) continue; + const double re = g(rng) * scale, im = g(rng) * scale; + f[i] = std::log2(std::hypot(re, im)) - 1.0; + f[half + i] = std::atan2(im, re) / std::log(2.0); + } + if (n % 2) f[n - 1] = 0.5; + } else if (kind == "smooth") { + for (size_t i = 0; i < n; ++i) f[i] = std::sin(i * 1e-3) + 0.5 * std::cos(i * 7e-5); + } else if (kind == "noise") { + for (size_t i = 0; i < n; ++i) f[i] = 2.0 * u(rng) - 1.0; + } else if (kind == "spiky") { + // Past the radius at small bounds, non-finite, subnormal, and values no index reconstructs within the bound. + const double inf = std::numeric_limits::infinity(); + for (size_t i = 0; i < n; ++i) { + const double r = u(rng); + f[i] = r < 0.001 ? 1e300 * g(rng) + : r < 0.002 ? 1e9 * g(rng) + : r < 0.0025 ? std::numeric_limits::quiet_NaN() + : r < 0.003 ? inf + : r < 0.0035 ? -inf + : r < 0.004 ? std::numeric_limits::denorm_min() * (1.0 + 1000.0 * u(rng)) + : r < 0.0045 ? -std::numeric_limits::min() * u(rng) + : 1e-3 * g(rng); + } + } // "zeros" stays all zero + return f; +} + +uint64_t bitsOf(double x) { + uint64_t u; + std::memcpy(&u, &x, 8); + return u; +} + +double doubleOf(uint64_t u) { + double x; + std::memcpy(&x, &u, 8); + return x; +} + +// The code fz::pfpl gives x at this bound. volatile keeps the host compiler from fusing the product into the +// subtraction: the decoder stores the rounded product, so that is what the bound applies to. +uint64_t expectedCode(double x, double eb) { + const double r = std::rint(x * (0.5 / eb)); + if (!(std::fabs(r) < static_cast(fz::pfpl::kQuantRadius64))) return bitsOf(x); + const volatile double product = r * (2.0 * eb); + if (!(std::fabs(x - product) <= eb)) return bitsOf(x); + const int64_t q = static_cast(r); + return (static_cast(q) << 1) ^ static_cast(q >> 63); +} + +double expectedValue(uint64_t code, double eb) { + if ((code >> 1) >= fz::pfpl::kQuantRadius64) return doubleOf(code); + const int64_t q = static_cast((code >> 1) ^ (0ull - (code & 1ull))); + return static_cast(q) * (2.0 * eb); +} + +// Everything after the quantizer, staged: difference, negabinary, 8-byte bitshuffle, RZE. +std::unique_ptr stagedTail64(size_t rawBytes) { + auto p = std::make_unique(rawBytes, fz::MemoryStrategy::PREALLOCATE); + p->setSpecializationPolicy(fz::SpecializationPolicy::Off); + auto* diff = p->addStage>(); + diff->setChunkSize(kChunkBytes); + auto* shuffle = p->addStage(); + shuffle->setElementWidth(8); + shuffle->setBlockSize(kChunkBytes); + auto* rze = p->addStage(); + rze->setWordSize(1); + rze->setChunkSize(kChunkBytes); + p->connect(shuffle, diff); + p->connect(rze, shuffle); + p->finalize(); + return p; +} + +void expectMatches64(const std::string& kind, size_t n, double eb, cudaStream_t s, fz::pfpl::Plan& plan) { + SCOPED_TRACE(kind + " n=" + std::to_string(n) + " eb=" + std::to_string(eb)); + const size_t raw = n * sizeof(double); + const size_t padded = (raw + kChunkBytes - 1) / kChunkBytes * kChunkBytes; + const std::vector h = makeField64(kind, n, n * 31 + static_cast(-std::log10(eb))); + + CudaBuffer d_in(n), d_out(n); + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + CudaBuffer d_len(1); + d_in.upload(h, s); + + fz::pfpl::compress(plan, d_in.get(), n, eb, d_stream.get(), d_stream.bytes(), d_len.get(), s); + const size_t len = readDeviceSize(d_len, s); + ASSERT_LE(len, d_stream.bytes()); + ASSERT_GE(len, sizeof(fz::pfpl::StreamHeader)); + const std::vector stream = d_stream.download_bytes(len, s); + EXPECT_EQ(fz::pfpl::streamBytes(stream.data(), len), len); + fz::pfpl::StreamHeader hdr; + std::memcpy(&hdr, stream.data(), sizeof(hdr)); + EXPECT_EQ(hdr.magic, fz::pfpl::kStreamMagic64); + EXPECT_EQ(hdr.n, n); + EXPECT_LE(static_cast(hdr.abs_eb), eb); + const double used = static_cast(hdr.abs_eb); + + // The expected codes, zero-padded to whole chunks as the field is coded. + std::vector codes(padded / sizeof(uint64_t), 0); + for (size_t i = 0; i < n; ++i) codes[i] = expectedCode(h[i], used); + CudaBuffer d_codes(codes.size()); + d_codes.upload(codes, s); + cudaStreamSynchronize(s); + + auto pipe = stagedTail64(padded); + const size_t cap = pipe->getMaxCompressedSize(padded); + CudaBuffer d_ref(cap); + size_t refLen = 0; + pipe->compress(d_codes.get(), padded, d_ref.get(), cap, &refLen, s); + cudaStreamSynchronize(s); + const std::vector ref = d_ref.download_bytes(refLen, s); + + const size_t arc = len - sizeof(fz::pfpl::StreamHeader); + ASSERT_EQ(arc, refLen) << "archive length differs from the staged stages'"; + EXPECT_EQ(std::memcmp(ref.data(), stream.data() + sizeof(fz::pfpl::StreamHeader), arc), 0) + << "archive bytes differ from the staged stages'"; + + fz::pfpl::decompress(plan, d_stream.get(), d_out.get(), n, s); + cudaStreamSynchronize(s); + const std::vector a = d_out.download_bytes(raw, s); + + size_t differ = 0, overBound = 0, nonFinite = 0; + double maxErr = 0.0; + for (size_t i = 0; i < n; ++i) { + if (bitsOf(a[i]) != bitsOf(expectedValue(codes[i], used))) ++differ; + if (!std::isfinite(h[i])) { + if (bitsOf(a[i]) != bitsOf(h[i])) ++nonFinite; + continue; + } + const double e = std::fabs(a[i] - h[i]); + maxErr = std::max(maxErr, e); + if (e > used) ++overBound; + } + EXPECT_EQ(differ, 0u) << "values differ from the codes' reconstruction"; + EXPECT_EQ(overBound, 0u) << "values past the bound; largest error " << maxErr; + EXPECT_EQ(nonFinite, 0u) << "non-finite values not restored"; +} + +class PfplStream64Case : public ::testing::TestWithParam {}; + +} // namespace + +// PD1 +TEST_P(PfplStream64Case, MatchesStagedStages) { + Stream s; + fz::pfpl::Plan plan; + expectMatches64(std::get<0>(GetParam()), std::get<1>(GetParam()), std::get<2>(GetParam()), s, plan); +} + +INSTANTIATE_TEST_SUITE_P( + Fields, PfplStream64Case, + ::testing::Combine(::testing::Values("polar", "smooth", "noise", "zeros", "spiky"), + ::testing::Values(size_t{256}, kElems64, 3 * kElems64 + 2, size_t{1} << 17, size_t{1} << 22), + ::testing::Values(1e-3, 1e-8, 1e-13)), + [](const ::testing::TestParamInfo& info) { + return std::string(std::get<0>(info.param)) + "_n" + std::to_string(std::get<1>(info.param)) + "_eb" + + std::to_string(static_cast(std::lround(-std::log10(std::get<2>(info.param))))); + }); + +// PD2 +TEST(PfplStream64, OverflowStaysInCapacity) { + Stream s; + fz::pfpl::Plan plan; + const size_t n = 3 * kElems64 + 77; + const std::vector h = makeField64("noise", n, 5); + CudaBuffer d_in(n); + CudaBuffer d_len(1); + d_in.upload(h, s); + + size_t full = 0; + { + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + fz::pfpl::compress(plan, d_in.get(), n, 1e-12, d_stream.get(), d_stream.bytes(), d_len.get(), s); + full = readDeviceSize(d_len, s); + } + ASSERT_GT(full, size_t{100}); + + for (const size_t cap : {full - 1, full / 2, size_t{100}, size_t{24}, size_t{8}}) { + SCOPED_TRACE("capacity " + std::to_string(cap)); + const size_t guard = 4096; + CudaBuffer d_stream(cap + guard); + ASSERT_EQ(cudaMemset(d_stream.get(), 0xA5, cap + guard), cudaSuccess); + fz::pfpl::compress(plan, d_in.get(), n, 1e-12, d_stream.get(), cap, d_len.get(), s); + EXPECT_EQ(readDeviceSize(d_len, s), full); + const std::vector all = d_stream.download(s); + const size_t touched = static_cast( + std::count_if(all.begin() + static_cast(cap), all.end(), [](uint8_t b) { return b != 0xA5; })); + EXPECT_EQ(touched, 0u) << "bytes written past the capacity"; + } +} + +// PD3 +TEST(PfplStream64, ConcurrentStreamsMatchSerial) { + constexpr int kStreams = 4, kRounds = 20; + constexpr size_t n = size_t{1} << 21; + constexpr double eb = 1e-7; + + std::vector> fields; + std::vector> streams; + std::vector> plans; + std::vector>> in, out; + std::vector>> st; + std::vector>> len; + std::vector> serial(kStreams); + for (int i = 0; i < kStreams; ++i) { + fields.push_back(makeField64(i % 2 ? "polar" : "smooth", n, 100 + i)); + streams.emplace_back(new Stream); + plans.emplace_back(new fz::pfpl::Plan); + in.emplace_back(new CudaBuffer(n)); + out.emplace_back(new CudaBuffer(n)); + st.emplace_back(new CudaBuffer(fz::pfpl::maxStreamBytes(n))); + len.emplace_back(new CudaBuffer(1)); + in[i]->upload(fields[i], *streams[i]); + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::compress(*plans[i], in[i]->get(), n, eb, st[i]->get(), st[i]->bytes(), len[i]->get(), *streams[i]); + serial[i] = st[i]->download_bytes(readDeviceSize(*len[i], *streams[i]), *streams[i]); + } + + for (int round = 0; round < kRounds; ++round) { + SCOPED_TRACE("round " + std::to_string(round)); + for (int i = 0; i < kStreams; ++i) { + ASSERT_EQ(cudaMemsetAsync(st[i]->get(), 0, st[i]->bytes(), *streams[i]), cudaSuccess); + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::compress(*plans[i], in[i]->get(), n, eb, st[i]->get(), st[i]->bytes(), len[i]->get(), *streams[i]); + } + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + for (int i = 0; i < kStreams; ++i) { + EXPECT_EQ(st[i]->download_bytes(serial[i].size(), *streams[i]), serial[i]) << "stream " << i; + } + for (int i = 0; i < kStreams; ++i) { + fz::pfpl::decompress(*plans[i], st[i]->get(), out[i]->get(), n, *streams[i]); + } + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + for (int i = 0; i < kStreams; ++i) { + const std::vector d = out[i]->download(*streams[i]); + size_t overBound = 0; + for (size_t j = 0; j < n; ++j) overBound += std::fabs(d[j] - fields[i][j]) > eb; + EXPECT_EQ(overBound, 0u) << "stream " << i; + } + } +} + +// PD4: a stream workspace's plan is sized in bytes, and one run could hold either type. +TEST(PfplStream64, SharesAPlanWithFloats) { + Stream s; + fz::pfpl::Plan plan; + plan.reserveBytes(size_t{1} << 23); + for (int r = 0; r < 30 && !::testing::Test::HasFatalFailure(); ++r) { + if (r % 2) expectMatchesStaged("polar", (size_t{1} << 20) + (r % 5) * 1000, 1e-3, s, plan); + else expectMatches64("polar", (size_t{1} << 19) + (r % 5) * 500, 1e-8, s, plan); + } +} + +// PD5 +TEST(PfplStream64, StreamBytesReadsTheHeader) { + const fz::pfpl::StreamHeader good{fz::pfpl::kStreamMagic64, 1000, 1e-3f, 2048}; + const size_t length = sizeof(good) + 1000; + EXPECT_EQ(fz::pfpl::streamBytes(&good), length); + EXPECT_EQ(fz::pfpl::streamBytes(&good, length), length); + EXPECT_EQ(fz::pfpl::streamBytes(&good, length - 1), 0u) << "claims more than the capacity"; + + // 8 chunks of doubles need a size table of 32 bytes, but as floats the same n is only 4 chunks. + const fz::pfpl::StreamHeader doubles{fz::pfpl::kStreamMagic64, 8 + 4 * 4, 1e-3f, 8 * 2048}; + const fz::pfpl::StreamHeader floats{fz::pfpl::kStreamMagic, 8 + 4 * 4, 1e-3f, 8 * 2048}; + EXPECT_EQ(fz::pfpl::streamBytes(&doubles), 0u) << "shorter than its size table"; + EXPECT_NE(fz::pfpl::streamBytes(&floats), 0u); + + EXPECT_EQ(fz::pfpl::maxStreamBytes(2048), fz::pfpl::maxStreamBytes(4096)); +} + +// PD6 +TEST(PfplStream64, BoundIsRecordedAsAFloat) { + Stream s; + fz::pfpl::Plan plan; + const size_t n = 2048; + const std::vector h = makeField64("smooth", n, 3); + CudaBuffer d_in(n); + CudaBuffer d_stream(fz::pfpl::maxStreamBytes(n)); + CudaBuffer d_len(1); + d_in.upload(h, s); + + for (const double eb : {1e-5, 0.1, 1.0 / 3.0, 3e-30}) { + SCOPED_TRACE(eb); + fz::pfpl::compress(plan, d_in.get(), n, eb, d_stream.get(), d_stream.bytes(), d_len.get(), s); + const size_t len = readDeviceSize(d_len, s); + fz::pfpl::StreamHeader hdr; + const std::vector head = d_stream.download_bytes(sizeof(hdr), s); + std::memcpy(&hdr, head.data(), sizeof(hdr)); + EXPECT_GT(len, sizeof(hdr)); + EXPECT_LE(static_cast(hdr.abs_eb), eb); + EXPECT_GE(static_cast(hdr.abs_eb), eb * (1.0 - 1.2e-7)); + } + const auto attempt = [&](double eb) { + fz::pfpl::compress(plan, d_in.get(), n, eb, d_stream.get(), d_stream.bytes(), d_len.get(), s); + }; + EXPECT_THROW(attempt(1e-40), std::invalid_argument); + EXPECT_THROW(attempt(0.0), std::invalid_argument); + EXPECT_THROW(attempt(-1.0), std::invalid_argument); + EXPECT_THROW(attempt(std::numeric_limits::infinity()), std::invalid_argument); + EXPECT_THROW(attempt(std::numeric_limits::quiet_NaN()), std::invalid_argument); +} From 1811d4515b0976f49fa006b05249ece92a619e6b Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Thu, 5 Nov 2026 15:59:33 -0500 Subject: [PATCH 6/8] fz::pfpl: record doubles in the changelog and notes CN-PFPLSTREAM-3 has the 64-bit chain, the radius, how the float-sized header records the bound, and the throughput on one MI250X GCD. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + docs/codebase_notes.md | 28 ++++++++++++++++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 56baa538..918c3553 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ Version numbers follow [Semantic Versioning](https://semver.org/). ### Added - `fz::pfpl` (`fused/chunk_fusion/pfpl_stream.h`): a float field compresses from one device buffer into another and back in one stream-ordered kernel each, with no host synchronization and the same archive bytes as the PFPL preset in ABS mode. - Wave64 AMD kernels for `fz::pfpl`, used by default on HIP; `FZ_PFPL_HARNESS=1` selects the generic harness kernels. See CN-PFPLSTREAM-1. +- `fz::pfpl` takes double fields as well: a 64-bit zigzag code, difference and negabinary, a 64-plane bitshuffle, and the same RZE over 16 KiB chunks. Double streams carry their own magic and read back only as doubles. See CN-PFPLSTREAM-3. - `-DFZGMOD_BACKEND=HIP` configures again, with `hip` CMake presets. NVRTC specialization is not built on HIP, so every stage runs staged, and `test_fusion_f64`, `test_fusion_fixed`, `test_fusion_planner` and `test_lorenzo_quant_shared_init` are CUDA-only. ### Fixed diff --git a/docs/codebase_notes.md b/docs/codebase_notes.md index 4211775a..8aeae55e 100644 --- a/docs/codebase_notes.md +++ b/docs/codebase_notes.md @@ -1859,3 +1859,31 @@ stream, ahead of that call's kernel, so they need no host synchronization. `PfplStream.FreshPlanIgnoresStaleMemory` keeps the default stream busy so that a clear queued there always loses the race. Without this fix it faults in its first round on both kernel sets. + +## CN-PFPLSTREAM-3 — Doubles in `fz::pfpl` + +**Source:** `quantCode64` and `encodeChunk64` in `modules/fused/chunk_fusion/pfpl_stream.cu`, `Layout` in `pfpl_stream_amd.inl`. + +A double takes the float chain's steps in 64 bits. The chunk stays 16 KiB, so it holds 2048 elements in 64 planes of 256 bytes, +and everything after the bitshuffle is unchanged. A 64-plane bitshuffle is two 32-plane transposes: a group's high words give +planes 0..31 and its low words planes 32..63, with the groups of a plane in swapped pairs. That is the staged 8-byte +`BitshuffleStage`'s order. + +The staged quantizer has no double form with in-place outliers, so no staged pipeline matches the whole chain. +`PfplStream64Case.MatchesStagedStages` computes the expected codes on the host and compares the archive with staged +`DifferenceStage`, `BitshuffleStage` and `RZEStage`, byte for byte. + +- **Quantization radius:** 2^51, the analog of the float preset's 2^22. A raw value's bits have to read as raw on decode, and + the bits of a normal double lie at or above 2 * radius. A subnormal always quantizes to index 0 at an allowed bound. +- **Bound:** the 16-byte header is unchanged, so the stream records the bound as a float, rounded down, and that rounded value + is the bound the field is coded to. It has to be at least the smallest normal float. A power-of-two bound is exact. +- **Rounding:** `rint`, with a range check before the conversion, which is undefined out of range. + +Measured on one MI250X GCD with ROCm 7.1.1, on a 128 MiB polar tile (throughput in GB/s of raw data): + +| Type, bound | Ratio | AMD encode | AMD decode | Generic encode | Generic decode | +|---|---:|---:|---:|---:|---:| +| float, 1e-4 | 1.74 | 186.4 | 164.5 | 62.5 | 55.1 | +| double, 1e-4 | 3.56 | 219.2 | 182.5 | 70.0 | 65.1 | +| double, 1e-9 | 1.81 | 214.2 | 177.2 | | | +| double, 1e-12 | 1.40 | 207.9 | 173.0 | 68.1 | 61.6 | From c7b048a94701bc2c1a23c9483534df34cdaa095c Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Thu, 5 Nov 2026 15:59:33 -0500 Subject: [PATCH 7/8] fz::pfpl: a host coder that writes the kernels' stream compressHost and decompressHost read and write the same stream as the kernels, byte for byte, in plain C++, with no device, Plan or GPU. A caller off the GPU path no longer has to stage a stream through the device: in QCLAB++ that ran out of memory reading back double tiles of 512 MiB, and a rank with no GPU could not code at all. PfplHostCase holds them to the kernels on both kernel sets, both types and every size class, including a hostile field for conversions out of range. compress() and compressHost() now share their checks and quantizer parameters through quantParams(). See CN-PFPLSTREAM-4. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + docs/codebase_notes.md | 22 ++ modules/fused/chunk_fusion/pfpl_stream.cu | 54 ++- modules/fused/chunk_fusion/pfpl_stream.h | 17 + .../fused/chunk_fusion/pfpl_stream_host.inl | 326 ++++++++++++++++++ tests/stages/test_pfpl_stream.cpp | 196 +++++++++++ 6 files changed, 607 insertions(+), 9 deletions(-) create mode 100644 modules/fused/chunk_fusion/pfpl_stream_host.inl diff --git a/CHANGELOG.md b/CHANGELOG.md index 918c3553..a74289db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ Version numbers follow [Semantic Versioning](https://semver.org/). - `fz::pfpl` (`fused/chunk_fusion/pfpl_stream.h`): a float field compresses from one device buffer into another and back in one stream-ordered kernel each, with no host synchronization and the same archive bytes as the PFPL preset in ABS mode. - Wave64 AMD kernels for `fz::pfpl`, used by default on HIP; `FZ_PFPL_HARNESS=1` selects the generic harness kernels. See CN-PFPLSTREAM-1. - `fz::pfpl` takes double fields as well: a 64-bit zigzag code, difference and negabinary, a 64-plane bitshuffle, and the same RZE over 16 KiB chunks. Double streams carry their own magic and read back only as doubles. See CN-PFPLSTREAM-3. +- `fz::pfpl::compressHost` and `decompressHost`: the same stream, byte for byte, in host memory, with no device, `Plan` or GPU. See CN-PFPLSTREAM-4. - `-DFZGMOD_BACKEND=HIP` configures again, with `hip` CMake presets. NVRTC specialization is not built on HIP, so every stage runs staged, and `test_fusion_f64`, `test_fusion_fixed`, `test_fusion_planner` and `test_lorenzo_quant_shared_init` are CUDA-only. ### Fixed diff --git a/docs/codebase_notes.md b/docs/codebase_notes.md index 8aeae55e..0909d17a 100644 --- a/docs/codebase_notes.md +++ b/docs/codebase_notes.md @@ -1887,3 +1887,25 @@ Measured on one MI250X GCD with ROCm 7.1.1, on a 128 MiB polar tile (throughput | double, 1e-4 | 3.56 | 219.2 | 182.5 | 70.0 | 65.1 | | double, 1e-9 | 1.81 | 214.2 | 177.2 | | | | double, 1e-12 | 1.40 | 207.9 | 173.0 | 68.1 | 61.6 | + +## CN-PFPLSTREAM-4 — The host coder + +**Source:** `modules/fused/chunk_fusion/pfpl_stream_host.inl`. + +A caller that has to read or write a stream off the GPU path used to stage it through the device. QCLAB++ did, and it failed +two ways. Reading a state back, 14 threads per rank each allocated a tile's worth of device memory beside what PaRSEC already +held, and double tiles of 512 MiB stopped with out of memory. And a rank that could see no GPU could not code at all. + +`compressHost` and `decompressHost` write and read the same stream as the kernels, in plain C++, one thread per call, with stack +buffers only. `PfplHostCase` holds them to the kernels byte for byte on both kernel sets, both types and every size class. Its +`hostile` field exists for the one place they can differ, a float-to-int conversion out of range: the GPU saturates and the host +conversion is undefined, so the host code tests the range first, and both then store the raw bits. `-ffp-contract` matters as +it does in the kernels (see `fmulRounded`), so the rebuilt value goes through a `volatile`. + +Measured on one Frontier login core, with no GPU in the process, on a 128 MiB polar tile: + +| Type, bound | Ratio | Encode | Decode | +|---|---:|---:|---:| +| float, 1e-3 | 2.43 | 350 MB/s | 550 MB/s | +| double, 1e-3 | 4.87 | 565 MB/s | 656 MB/s | +| double, 1e-9 | 1.92 | 460 MB/s | 638 MB/s | diff --git a/modules/fused/chunk_fusion/pfpl_stream.cu b/modules/fused/chunk_fusion/pfpl_stream.cu index 61f68f46..f2dfb5db 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.cu +++ b/modules/fused/chunk_fusion/pfpl_stream.cu @@ -577,8 +577,10 @@ struct PlanAccess { } }; -void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, - void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { +namespace { + +// What compress() and compressHost() share: the checks, and the quantizer's parameters. +QuantInplaceZigzagParams quantParams(float abs_eb, size_t n) { if (!(abs_eb > 0.0f) || !std::isfinite(abs_eb) || !std::isfinite(0.5f / abs_eb)) throw std::invalid_argument("fz::pfpl::compress: the error bound must be positive and finite, " "with a finite quantizer scale 0.5/abs_eb"); @@ -594,27 +596,41 @@ void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, qp.radius = kQuantRadius; qp.threshold = std::numeric_limits::infinity(); qp.verify_reconstruction = 1; - PlanAccess::encode(plan, d_in, n, qp, abs_eb, d_out, capacity, d_stream_bytes, stream); + return qp; } -void compress(Plan& plan, const double* d_in, size_t n, double abs_eb, - void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { +QuantParams64 quantParams(double abs_eb, size_t n, float* recorded) { if (!(abs_eb > 0.0) || !std::isfinite(abs_eb)) throw std::invalid_argument("fz::pfpl::compress: the error bound must be positive and finite"); // The bound the stream records: never above the one asked for. - float recorded = static_cast(abs_eb); - if (static_cast(recorded) > abs_eb) recorded = std::nextafter(recorded, 0.0f); - if (recorded < std::numeric_limits::min()) + *recorded = static_cast(abs_eb); + if (static_cast(*recorded) > abs_eb) *recorded = std::nextafter(*recorded, 0.0f); + if (*recorded < std::numeric_limits::min()) throw std::invalid_argument("fz::pfpl::compress: a double stream records its bound as a float, so the bound " "must be at least the smallest normal float"); if (!fitsArchive(n)) throw std::invalid_argument("fz::pfpl::compress: a field of 4 GiB or more (in whole chunks) exceeds " "the LC archive's 32-bit sizes"); - const double eb = static_cast(recorded); + const double eb = static_cast(*recorded); QuantParams64 qp{}; qp.ebx2_r = 0.5 / eb; qp.ebx2 = 2.0 * eb; qp.abs_eb = eb; + return qp; +} + +} // namespace + +void compress(Plan& plan, const float* d_in, size_t n, float abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { + const QuantInplaceZigzagParams qp = quantParams(abs_eb, n); + PlanAccess::encode(plan, d_in, n, qp, abs_eb, d_out, capacity, d_stream_bytes, stream); +} + +void compress(Plan& plan, const double* d_in, size_t n, double abs_eb, + void* d_out, size_t capacity, size_t* d_stream_bytes, fz::stream_t stream) { + float recorded = 0.0f; + const QuantParams64 qp = quantParams(abs_eb, n, &recorded); PlanAccess::encode(plan, d_in, n, qp, recorded, d_out, capacity, d_stream_bytes, stream); } @@ -626,5 +642,25 @@ void decompress(Plan& plan, const void* d_stream, double* d_out, size_t n, fz::s PlanAccess::decode(plan, d_stream, d_out, n, stream); } +#include "fused/chunk_fusion/pfpl_stream_host.inl" + +size_t compressHost(const float* in, size_t n, float abs_eb, void* out, size_t capacity) { + return host::encode(in, n, quantParams(abs_eb, n), abs_eb, out, capacity); +} + +size_t compressHost(const double* in, size_t n, double abs_eb, void* out, size_t capacity) { + float recorded = 0.0f; + const QuantParams64 qp = quantParams(abs_eb, n, &recorded); + return host::encode(in, n, qp, recorded, out, capacity); +} + +bool decompressHost(const void* stream, float* out, size_t n, size_t outBytes) { + return host::decode(stream, out, n, outBytes); +} + +bool decompressHost(const void* stream, double* out, size_t n, size_t outBytes) { + return host::decode(stream, out, n, outBytes); +} + } // namespace pfpl } // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream.h b/modules/fused/chunk_fusion/pfpl_stream.h index 79029525..3a26791b 100644 --- a/modules/fused/chunk_fusion/pfpl_stream.h +++ b/modules/fused/chunk_fusion/pfpl_stream.h @@ -124,5 +124,22 @@ void compress(Plan& plan, const double* d_in, size_t n, double abs_eb, void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream); void decompress(Plan& plan, const void* d_stream, double* d_out, size_t n, fz::stream_t stream); +/** + * compress() on the host: the same stream, byte for byte, in host memory, from any thread, with no device, Plan or + * GPU. One thread, so far slower than the kernels. Returns the stream's length, greater than `capacity` when it does + * not fit; nothing at or past out + capacity is written, and a stream that does not fit leaves its 16-byte header + * (when that fits), as compress() does. The bound is taken as compress() takes it. + */ +size_t compressHost(const float* in, size_t n, float abs_eb, void* out, size_t capacity); +size_t compressHost(const double* in, size_t n, double abs_eb, void* out, size_t capacity); + +/** + * decompress() on the host, under the same terms: the stream at `stream` (host memory; compress() or compressHost() + * wrote it) must describe n elements of this type. Writes the first `outBytes` bytes of the field and stops; false, with + * nothing more written, if the stream is not such a stream or its sizes do not add up. + */ +bool decompressHost(const void* stream, float* out, size_t n, size_t outBytes = SIZE_MAX); +bool decompressHost(const void* stream, double* out, size_t n, size_t outBytes = SIZE_MAX); + } // namespace pfpl } // namespace fz diff --git a/modules/fused/chunk_fusion/pfpl_stream_host.inl b/modules/fused/chunk_fusion/pfpl_stream_host.inl new file mode 100644 index 00000000..1966a8d0 --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream_host.inl @@ -0,0 +1,326 @@ +// The coder on the host: the stream the kernels write, byte for byte, from plain C++. No device, stream or Plan, so any +// thread may call it, and it needs no GPU. One thread, one chunk at a time. Part of pfpl_stream.cu. +namespace host { + +constexpr uint32_t kAvail = kChunkBytes - 2; // d_RZE: room before its 2-byte trailer +constexpr int kL1 = kChunkBytes / 8; // zero-elimination bitmap bytes +constexpr int kL2 = kL1 / 8; +constexpr int kL3 = kL2 / 8; +constexpr int kL4 = kL3 / 8; +constexpr int kWords = kChunkBytes / 4; // 32-bit words of a bitshuffled chunk + +// What the kernels' butterfly32 computes over 32 lanes: word p's bit (31 - r) is word r's bit (31 - p). Its own inverse. +inline void transpose32(uint32_t* a) { + uint32_t m = 0x0000ffffu; + for (int j = 16; j != 0; j >>= 1, m ^= (m << j)) + for (int k = 0; k < 32; k = (k + j + 1) & ~j) { + const uint32_t t = (a[k] ^ (a[k + j] >> j)) & m; + a[k] ^= t; + a[k + j] ^= (t << j); + } +} + +// The range check comes before the conversion: out of int's range it is undefined here, where the GPU saturates, and +// either way the value is stored raw. `volatile` keeps a compiler from fusing the product into x - q * ebx2: the bound +// check reads the rounded product, as fz::backend::fmulRounded makes the kernels do. +inline uint32_t codeOf(float x, const QuantInplaceZigzagParams& p) { + const float r = std::round(x * p.ebx2_r); + if (std::fabs(r) < static_cast(p.radius) && std::fabs(x) < p.threshold) { + const int q = static_cast(r); + volatile float rebuilt = static_cast(q) * p.ebx2; + const float value = rebuilt; + if (std::fabs(x - value) <= p.abs_eb) + return (static_cast(q) << 1) ^ static_cast(q >> 31); + } + uint32_t bits; + std::memcpy(&bits, &x, sizeof(bits)); + return bits; +} + +inline uint64_t codeOf(double x, const QuantParams64& p) { + const double r = std::rint(x * p.ebx2_r); + if (std::fabs(r) < static_cast(kQuantRadius64)) { + volatile double rebuilt = r * p.ebx2; + const double value = rebuilt; + if (std::fabs(x - value) <= p.abs_eb) { + const long long q = static_cast(r); + return (static_cast(q) << 1) ^ static_cast(q >> 63); + } + } + uint64_t bits; + std::memcpy(&bits, &x, sizeof(bits)); + return bits; +} + +inline float valueOf(uint32_t code, float ebx2) { + if ((code >> 1) >= kQuantRadius) { + float x; + std::memcpy(&x, &code, sizeof(x)); + return x; + } + return static_cast(static_cast((code >> 1) ^ (0u - (code & 1u)))) * ebx2; +} + +inline double valueOf(uint64_t code, double ebx2) { + if ((code >> 1) >= kQuantRadius64) { + double x; + std::memcpy(&x, &code, sizeof(x)); + return x; + } + return static_cast(static_cast((code >> 1) ^ (0ull - (code & 1ull)))) * ebx2; +} + +// Negabinary codes to the bitshuffled chunk's words, and back: plane by plane, most significant first. A float chunk's +// plane is one word per group of 32 elements. A double's high words make planes 0..31 and its low words planes 32..63, +// and a plane holds its groups' words in pairs swapped (the staged 8-byte bitshuffle's order). +inline void shuffle(const uint32_t* nb, uint32_t* words) { + for (int g = 0; g < kWords / 32; ++g) { + uint32_t a[32]; + std::memcpy(a, nb + 32 * g, sizeof(a)); + transpose32(a); + for (int p = 0; p < 32; ++p) words[p * (kWords / 32) + g] = a[p]; + } +} + +inline void unshuffle(const uint32_t* words, uint32_t* nb) { + for (int g = 0; g < kWords / 32; ++g) { + uint32_t a[32]; + for (int p = 0; p < 32; ++p) a[p] = words[p * (kWords / 32) + g]; + transpose32(a); + std::memcpy(nb + 32 * g, a, sizeof(a)); + } +} + +inline void shuffle(const uint64_t* nb, uint32_t* words) { + constexpr int kNpp = kWords / 64; + for (int g = 0; g < kNpp; ++g) { + uint32_t hi[32], lo[32]; + for (int r = 0; r < 32; ++r) { + hi[r] = static_cast(nb[32 * g + r] >> 32); + lo[r] = static_cast(nb[32 * g + r]); + } + transpose32(hi); + transpose32(lo); + for (int p = 0; p < 32; ++p) { + words[p * kNpp + (g ^ 1)] = hi[p]; + words[(32 + p) * kNpp + (g ^ 1)] = lo[p]; + } + } +} + +inline void unshuffle(const uint32_t* words, uint64_t* nb) { + constexpr int kNpp = kWords / 64; + for (int g = 0; g < kNpp; ++g) { + uint32_t hi[32], lo[32]; + for (int p = 0; p < 32; ++p) { + hi[p] = words[p * kNpp + (g ^ 1)]; + lo[p] = words[(32 + p) * kNpp + (g ^ 1)]; + } + transpose32(hi); + transpose32(lo); + for (int r = 0; r < 32; ++r) nb[32 * g + r] = (static_cast(hi[r]) << 32) | lo[r]; + } +} + +inline uint32_t popcountBytes(const uint8_t* b, int n) { + uint32_t c = 0; + for (int i = 0; i < n; ++i) c += static_cast(__builtin_popcount(b[i])); + return c; +} + +// One level of repetition elimination: bit e of `bitmap` says byte e differs from the one before it (byte -1 is 0), and +// the bytes so marked go to `kept`. Returns how many. +inline uint32_t repeatEliminate(const uint8_t* src, int n, uint8_t* bitmap, uint8_t* kept) { + std::memset(bitmap, 0, static_cast(n) / 8); + uint32_t count = 0; + uint8_t prev = 0; + for (int e = 0; e < n; ++e) { + if (src[e] != prev) { + bitmap[e >> 3] = static_cast(bitmap[e >> 3] | (1u << (e & 7))); + kept[count++] = src[e]; + prev = src[e]; + } + } + return count; +} + +// Its inverse. +inline void repeatRestore(const uint8_t* kept, const uint8_t* bitmap, int n, uint8_t* out) { + uint8_t prev = 0; + uint32_t k = 0; + for (int e = 0; e < n; ++e) { + if ((bitmap[e >> 3] >> (e & 7)) & 1u) prev = kept[k++]; + out[e] = prev; + } +} + +// d_RZE on one bitshuffled chunk: its packed form in `packed` and that length, with the kernels' fall-back checks in +// their order, or 0 if the chunk is stored raw. +inline uint32_t packChunk(const uint8_t* chunk, uint8_t* packed) { + uint8_t l1[kL1] = {}; + uint32_t wpos = 0; + for (int i = 0; i < kChunkBytes; ++i) { + if (chunk[i] != 0) { + l1[i >> 3] = static_cast(l1[i >> 3] | (1u << (i & 7))); + packed[wpos++] = chunk[i]; + } + } + if (wpos >= kAvail) return 0; + uint32_t csize = 2; + if (wpos != 0) { + uint8_t l2[kL2], l3[kL3], l4[kL4], re1[kL1], re2[kL2], re3[kL3]; + const uint32_t n1 = repeatEliminate(l1, kL1, l2, re1); + if (n1 > kAvail - wpos) return 0; + const uint32_t n2 = repeatEliminate(l2, kL2, l3, re2); + if (n2 > kAvail - wpos - n1) return 0; + const uint32_t n3 = repeatEliminate(l3, kL3, l4, re3); + if (n3 > kAvail - wpos - n1 - n2) return 0; + const uint32_t w3 = wpos + n1 + n2 + n3; + if (w3 >= kAvail - 4) return 0; + std::memcpy(packed + wpos, re1, n1); + std::memcpy(packed + wpos + n1, re2, n2); + std::memcpy(packed + wpos + n1 + n2, re3, n3); + std::memcpy(packed + w3, l4, kL4); + csize = w3 + kL4 + 2; + } + packed[csize - 2] = static_cast(kChunkBytes & 0xff); + packed[csize - 1] = static_cast(kChunkBytes >> 8); + return csize; +} + +// The bitshuffled chunk from its packed form. False if the lengths in it do not add up to `stored`. +inline bool unpackChunk(const uint8_t* payload, uint32_t stored, uint8_t* chunk) { + if (stored < 2 || stored >= kChunkBytes) return false; + uint32_t rpos = stored - 2; + if (rpos == 0) { + std::memset(chunk, 0, kChunkBytes); + return true; + } + if (rpos < kL4) return false; + rpos -= kL4; + uint8_t l1[kL1], l2[kL2], l3[kL3]; + const uint8_t* l4 = payload + rpos; + const uint32_t n3 = popcountBytes(l4, kL4); + if (rpos < n3) return false; + rpos -= n3; + repeatRestore(payload + rpos, l4, kL3, l3); + const uint32_t n2 = popcountBytes(l3, kL3); + if (rpos < n2) return false; + rpos -= n2; + repeatRestore(payload + rpos, l3, kL2, l2); + const uint32_t n1 = popcountBytes(l2, kL2); + if (rpos < n1) return false; + rpos -= n1; + repeatRestore(payload + rpos, l2, kL1, l1); + if (popcountBytes(l1, kL1) != rpos) return false; + uint32_t j = 0; + for (int i = 0; i < kChunkBytes; ++i) chunk[i] = ((l1[i >> 3] >> (i & 7)) & 1u) ? payload[j++] : uint8_t(0); + return true; +} + +template +size_t encode(const Real* in, size_t n, const typename Elem::Params& params, float abs_eb, void* out_, size_t capacity) { + using Word = typename Elem::Word; + using SWord = typename Elem::SWord; + constexpr size_t kPer = kChunkBytes / sizeof(Real); + uint8_t* out = static_cast(out_); + const uint32_t chunks = chunksOf(n); + const uint32_t magic = Elem::kMagic; + uint32_t boundBits; + std::memcpy(&boundBits, &abs_eb, sizeof(boundBits)); + if (chunks == 0) { + const uint32_t w[6] = {magic, kLcHeader, boundBits, 0u, 0u, 0u}; + if (kStreamHeader <= capacity) std::memcpy(out, w, kStreamHeader); + if (kStreamHeader + kLcHeader <= capacity) std::memcpy(out + kStreamHeader, w + 4, kLcHeader); + return kStreamHeader + kLcHeader; + } + const unsigned long long table = kStreamHeader + kLcHeader; + const unsigned long long data = table + 4ull * chunks; + unsigned long long offset = 0; + Word nb[kPer]; + uint32_t words[kWords]; + uint8_t packed[kChunkBytes]; + for (uint32_t cid = 0; cid < chunks; ++cid) { + const size_t base = static_cast(cid) * kPer; + Word prev = 0; + for (size_t i = 0; i < kPer; ++i) { + const Word code = codeOf(base + i < n ? in[base + i] : static_cast(0), params); + nb[i] = Negabinary::encode(static_cast(code - prev)); + prev = code; + } + shuffle(nb, words); + const uint8_t* chunk = reinterpret_cast(words); + const uint32_t csize = packChunk(chunk, packed); + const uint32_t entry = csize != 0 ? csize : (0x80000000u | static_cast(kChunkBytes)); + const uint32_t nbytes = entry & 0x7fffffffu; + if (table + 4ull * (cid + 1u) <= capacity) std::memcpy(out + table + 4ull * cid, &entry, sizeof(entry)); + if (data + offset + nbytes <= capacity) std::memcpy(out + data + offset, csize != 0 ? packed : chunk, nbytes); + offset += nbytes; + } + const unsigned long long archive = kLcHeader + 4ull * chunks + offset; + const unsigned long long padded = (archive + 3ull) & ~3ull; + const unsigned long long total = kStreamHeader + padded; + if (total <= capacity) { + for (unsigned long long b = kStreamHeader + archive; b < total; ++b) out[b] = 0; + const uint32_t lc[2] = {chunks * static_cast(kChunkBytes), chunks}; + std::memcpy(out + kStreamHeader, lc, sizeof(lc)); + } + if (kStreamHeader <= capacity) { + const uint32_t w[4] = {magic, static_cast(padded), boundBits, static_cast(n)}; + std::memcpy(out, w, sizeof(w)); + } + return static_cast(total); +} + +template +bool decode(const void* stream_, Real* out, size_t n, size_t outBytes) { + using Word = typename Elem::Word; + using SWord = typename Elem::SWord; + constexpr size_t kPer = kChunkBytes / sizeof(Real); + const uint8_t* stream = static_cast(stream_); + if (stream == nullptr || streamBytes(stream) == 0 || !fitsArchive(n)) return false; + const uint32_t chunks = chunksOf(n); + uint32_t w[6]; + std::memcpy(w, stream, sizeof(w)); + float bound; + std::memcpy(&bound, &w[2], sizeof(bound)); + if (w[0] != Elem::kMagic || w[3] != n || w[4] != chunks * static_cast(kChunkBytes) || w[5] != chunks || + w[1] < kLcHeader + 4u * chunks || !(bound > 0.0f)) + return false; + const Real ebx2 = static_cast(2) * static_cast(bound); + const unsigned long long end = kLcHeader + 4ull * chunks; + const size_t want = std::min(outBytes, n * sizeof(Real)); + uint8_t* to = reinterpret_cast(out); + size_t written = 0; + unsigned long long offset = 0; + uint32_t words[kWords]; + Word nb[kPer]; + Real values[kPer]; + for (uint32_t cid = 0; cid < chunks && written < want; ++cid) { + uint32_t entry; + std::memcpy(&entry, stream + kStreamHeader + kLcHeader + 4ull * cid, sizeof(entry)); + const uint32_t stored = entry & 0x7fffffffu; + const bool raw = (entry & 0x80000000u) != 0u; + if (end + offset + stored > w[1]) return false; + const uint8_t* payload = stream + kStreamHeader + end + offset; + offset += stored; + if (raw) { + if (stored != kChunkBytes) return false; + std::memcpy(words, payload, kChunkBytes); + } else if (!unpackChunk(payload, stored, reinterpret_cast(words))) { + return false; + } + unshuffle(words, nb); + Word code = 0; + for (size_t i = 0; i < kPer; ++i) { + code += static_cast(Negabinary::decode(nb[i])); + values[i] = valueOf(code, ebx2); + } + const size_t bytes = std::min(want - written, kPer * sizeof(Real)); + std::memcpy(to + written, values, bytes); + written += bytes; + } + return written == want; +} + +} // namespace host diff --git a/tests/stages/test_pfpl_stream.cpp b/tests/stages/test_pfpl_stream.cpp index 01a8bbc0..c1b2df35 100644 --- a/tests/stages/test_pfpl_stream.cpp +++ b/tests/stages/test_pfpl_stream.cpp @@ -38,6 +38,7 @@ #include "shufflers/bitshuffle/bitshuffle_stage.h" #include +#include #include #include #include @@ -45,6 +46,7 @@ #include #include #include +#include #include #include @@ -694,3 +696,197 @@ TEST(PfplStream64, BoundIsRecordedAsAFloat) { EXPECT_THROW(attempt(std::numeric_limits::infinity()), std::invalid_argument); EXPECT_THROW(attempt(std::numeric_limits::quiet_NaN()), std::invalid_argument); } + +// ─── PH: the host coder against the kernels ───────────────────────────────────────────── + +namespace { + +// Everything a tile can hold that sends the quantizer down a different branch: raw values on either side of the radius, +// non-finite ones, subnormals, ties, and runs of zeros between them. +template +std::vector hostileField(size_t n, double eb, uint64_t seed) { + std::mt19937_64 rng(seed); + std::normal_distribution g(0.0, 1.0); + std::uniform_int_distribution pick(0, 15); + const double step = 2.0 * eb; + const double radius = sizeof(Real) == 4 ? static_cast(fz::pfpl::kQuantRadius) + : static_cast(fz::pfpl::kQuantRadius64); + std::vector f(n); + for (size_t i = 0; i < n; ++i) { + double v = 0.0; + switch (pick(rng)) { + case 0: v = radius * step * (1.0 + 1e-3 * g(rng)); break; // on the radius + case 1: v = -radius * step * (1.0 + 1e-3 * g(rng)); break; + case 2: v = std::numeric_limits::quiet_NaN(); break; + case 3: v = std::numeric_limits::infinity(); break; + case 4: v = -std::numeric_limits::infinity(); break; + case 5: v = std::numeric_limits::denorm_min(); break; + case 6: v = -std::numeric_limits::min() * 0.5; break; // subnormal + case 7: v = (std::floor(g(rng) * 1000.0) + 0.5) * step; break; // a rounding tie + case 8: v = std::numeric_limits::max() * (g(rng) > 0 ? 1.0 : -1.0); break; + case 9: v = g(rng) * 1e30; break; // far past int's range + case 10: case 11: break; // zero + default: v = g(rng); break; + } + f[i] = static_cast(v); + } + return f; +} + +template +std::vector fieldOf(const std::string& kind, size_t n, double eb, uint64_t seed); +template <> +std::vector fieldOf(const std::string& kind, size_t n, double eb, uint64_t seed) { + return kind == "hostile" ? hostileField(n, eb, seed) : makeField(kind, n, seed); +} +template <> +std::vector fieldOf(const std::string& kind, size_t n, double eb, uint64_t seed) { + return kind == "hostile" ? hostileField(n, eb, seed) : makeField64(kind, n, seed); +} + +template +void expectHostMatchesDevice(const std::string& kind, size_t n, double eb, cudaStream_t s, fz::pfpl::Plan& plan) { + SCOPED_TRACE(kind + " n=" + std::to_string(n) + " eb=" + std::to_string(eb)); + const std::vector h = fieldOf(kind, n, eb, n * 53 + static_cast(-std::log10(eb))); + const size_t cap = fz::pfpl::maxStreamBytes(n); + const size_t raw = n * sizeof(Real); + + CudaBuffer d_in(n), d_out(n); + CudaBuffer d_stream(cap); + CudaBuffer d_len(1); + d_in.upload(h, s); + fz::pfpl::compress(plan, d_in.get(), n, static_cast(eb), d_stream.get(), cap, d_len.get(), s); + const size_t len = readDeviceSize(d_len, s); + const std::vector device = d_stream.download_bytes(len, s); + + std::vector hostStream(cap, 0xcd); + const size_t hostLen = fz::pfpl::compressHost(h.data(), n, static_cast(eb), hostStream.data(), cap); + ASSERT_EQ(hostLen, len) << "host stream length differs from the device's"; + ASSERT_EQ(std::memcmp(hostStream.data(), device.data(), len), 0) << "host stream bytes differ from the device's"; + + // Decoded by the kernels and on the host, from the same stream: the same bits. + fz::pfpl::decompress(plan, d_stream.get(), d_out.get(), n, s); + cudaStreamSynchronize(s); + const std::vector a = d_out.download_bytes(raw, s); + std::vector b(n, static_cast(-7)); + ASSERT_TRUE(fz::pfpl::decompressHost(hostStream.data(), b.data(), n)); + ASSERT_EQ(std::memcmp(a.data(), b.data(), raw), 0) << "host decode differs from the device's"; + + // The first part of the field alone, ending inside an element, and nothing past it. + const size_t part = raw / 3 + 1; + std::vector partial(raw + 16, 0xab); + ASSERT_TRUE(fz::pfpl::decompressHost(hostStream.data(), reinterpret_cast(partial.data()), n, part)); + EXPECT_EQ(std::memcmp(partial.data(), b.data(), part), 0); + for (size_t i = part; i < partial.size(); ++i) ASSERT_EQ(partial[i], 0xab) << "wrote past the bytes asked for"; +} + +} // namespace + +class PfplHostCase : public ::testing::TestWithParam {}; + +TEST_P(PfplHostCase, FloatStreamEqualsTheKernels) { + const auto& [kind, n, eb] = GetParam(); + cudaStream_t s = nullptr; + ASSERT_EQ(cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking), cudaSuccess); + { + fz::pfpl::Plan plan; + expectHostMatchesDevice(kind, n, eb, s, plan); + } + cudaStreamDestroy(s); +} + +TEST_P(PfplHostCase, DoubleStreamEqualsTheKernels) { + const auto& [kind, n, eb] = GetParam(); + cudaStream_t s = nullptr; + ASSERT_EQ(cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking), cudaSuccess); + { + fz::pfpl::Plan plan; + expectHostMatchesDevice(kind, n, eb, s, plan); + } + cudaStreamDestroy(s); +} + +INSTANTIATE_TEST_SUITE_P( + Fields, PfplHostCase, + ::testing::Combine(::testing::Values("polar", "smooth", "noise", "zeros", "spiky", "hostile"), + ::testing::Values(1, 7, 256, 2048, 4096, 4097, 3 * 2048 + 2, 1u << 17, 1u << 20), + ::testing::Values(1e-1, 1e-3, 1e-5))); + +template +void hostOverflowStaysInCapacity() { + const size_t n = 1u << 17; + const std::vector h = fieldOf("noise", n, 1e-5, 91); + const size_t needed = fz::pfpl::compressHost(h.data(), n, static_cast(1e-5), nullptr, 0); + ASSERT_GT(needed, sizeof(fz::pfpl::StreamHeader)); + for (const size_t cap : {size_t(0), size_t(8), sizeof(fz::pfpl::StreamHeader), size_t(32), needed / 2, needed - 1, needed}) { + std::vector out(cap + 64, 0xee); + const size_t len = fz::pfpl::compressHost(h.data(), n, static_cast(1e-5), out.data(), cap); + EXPECT_EQ(len, needed) << "capacity " << cap; + for (size_t i = cap; i < out.size(); ++i) ASSERT_EQ(out[i], 0xee) << "wrote past capacity " << cap; + if (cap >= sizeof(fz::pfpl::StreamHeader)) { + fz::pfpl::StreamHeader hdr; + std::memcpy(&hdr, out.data(), sizeof(hdr)); + EXPECT_EQ(sizeof(hdr) + static_cast(hdr.archive_bytes), needed) << "the header must record the length"; + } + } +} + +TEST(PfplHost, OverflowStaysInCapacityFloat) { hostOverflowStaysInCapacity(); } +TEST(PfplHost, OverflowStaysInCapacityDouble) { hostOverflowStaysInCapacity(); } + +TEST(PfplHost, RefusesAStreamNotOfThisField) { + const size_t n = 1u << 14; + const std::vector f = fieldOf("polar", n, 1e-3, 5); + std::vector d = fieldOf("polar", n, 1e-3, 5); + std::vector sf(fz::pfpl::maxStreamBytes(n)), sd(fz::pfpl::maxStreamBytes(n)); + const size_t lf = fz::pfpl::compressHost(f.data(), n, 1e-3f, sf.data(), sf.size()); + const size_t ld = fz::pfpl::compressHost(d.data(), n, 1e-3, sd.data(), sd.size()); + ASSERT_LE(lf, sf.size()); + ASSERT_LE(ld, sd.size()); + std::vector of(n); + std::vector od(n); + EXPECT_TRUE(fz::pfpl::decompressHost(sf.data(), of.data(), n)); + EXPECT_TRUE(fz::pfpl::decompressHost(sd.data(), od.data(), n)); + EXPECT_FALSE(fz::pfpl::decompressHost(sd.data(), of.data(), n)) << "a double stream read as floats"; + EXPECT_FALSE(fz::pfpl::decompressHost(sf.data(), od.data(), n)) << "a float stream read as doubles"; + EXPECT_FALSE(fz::pfpl::decompressHost(sf.data(), of.data(), n - 1)) << "another field's length"; + EXPECT_FALSE(fz::pfpl::decompressHost(nullptr, of.data(), n)); + std::vector zero(64, 0); + EXPECT_FALSE(fz::pfpl::decompressHost(zero.data(), of.data(), n)); + + // A chunk whose packed form lies about its lengths, and a size table that runs past the archive. + std::vector bad(sf.begin(), sf.begin() + static_cast(lf)); + const size_t table = sizeof(fz::pfpl::StreamHeader) + 8; + uint32_t entry; + std::memcpy(&entry, bad.data() + table, sizeof(entry)); + if ((entry & 0x80000000u) == 0u) { + const uint32_t shorter = 5; + std::memcpy(bad.data() + table, &shorter, sizeof(shorter)); + EXPECT_FALSE(fz::pfpl::decompressHost(bad.data(), of.data(), n)); + } + const uint32_t tooLong = 0x7ffffff0u; + std::memcpy(bad.data() + table, &tooLong, sizeof(tooLong)); + EXPECT_FALSE(fz::pfpl::decompressHost(bad.data(), of.data(), n)); +} + +TEST(PfplHost, ThreadsShareNothing) { + const size_t n = 1u << 16; + const std::vector f = fieldOf("polar", n, 1e-3, 17); + std::vector serial(fz::pfpl::maxStreamBytes(n)); + const size_t len = fz::pfpl::compressHost(f.data(), n, 1e-3f, serial.data(), serial.size()); + std::vector threads; + std::atomic differ{0}; + for (int t = 0; t < 8; ++t) + threads.emplace_back([&] { + std::vector out(serial.size()); + std::vector back(n); + for (int r = 0; r < 4; ++r) { + if (fz::pfpl::compressHost(f.data(), n, 1e-3f, out.data(), out.size()) != len || + std::memcmp(out.data(), serial.data(), len) != 0 || + !fz::pfpl::decompressHost(out.data(), back.data(), n)) + differ.fetch_add(1); + } + }); + for (auto& t : threads) t.join(); + EXPECT_EQ(differ.load(), 0); +} From 74855a40c868860ac7cb77685ab9c2563d26dbd4 Mon Sep 17 00:00:00 2001 From: Seth Wolfgang Date: Fri, 6 Nov 2026 02:40:23 +0000 Subject: [PATCH 8/8] docs: state the host coder's reason in the library's terms CN-PFPLSTREAM-4 named the caller it came from. Say what any caller needs, and call the measured data a field, as the rest of the library does. Co-Authored-By: Claude Opus 5.5 --- docs/codebase_notes.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/codebase_notes.md b/docs/codebase_notes.md index 0909d17a..32ef0d8c 100644 --- a/docs/codebase_notes.md +++ b/docs/codebase_notes.md @@ -1827,7 +1827,7 @@ time at 8 waves per CU. Their 36 KiB of shared memory also limits a CU to one block. The AMD kernels produce the same stream, byte for byte, with 6-7 barriers and 19 KiB of shared memory, which fits 3 blocks per CU. -Measured on one MI250X GCD with ROCm 7.1.1, on a polar tile with `abs_eb` 1e-4 +Measured on one MI250X GCD with ROCm 7.1.1, on a polar field with `abs_eb` 1e-4 and a compression ratio of about 1.7. Throughput is in GB/s: | Field | Encode: two-array | Encode: staged | Decode: two-array | Decode: staged | @@ -1879,7 +1879,7 @@ The staged quantizer has no double form with in-place outliers, so no staged pip is the bound the field is coded to. It has to be at least the smallest normal float. A power-of-two bound is exact. - **Rounding:** `rint`, with a range check before the conversion, which is undefined out of range. -Measured on one MI250X GCD with ROCm 7.1.1, on a 128 MiB polar tile (throughput in GB/s of raw data): +Measured on one MI250X GCD with ROCm 7.1.1, on a 128 MiB polar field (throughput in GB/s of raw data): | Type, bound | Ratio | AMD encode | AMD decode | Generic encode | Generic decode | |---|---:|---:|---:|---:|---:| @@ -1892,9 +1892,9 @@ Measured on one MI250X GCD with ROCm 7.1.1, on a 128 MiB polar tile (throughput **Source:** `modules/fused/chunk_fusion/pfpl_stream_host.inl`. -A caller that has to read or write a stream off the GPU path used to stage it through the device. QCLAB++ did, and it failed -two ways. Reading a state back, 14 threads per rank each allocated a tile's worth of device memory beside what PaRSEC already -held, and double tiles of 512 MiB stopped with out of memory. And a rank that could see no GPU could not code at all. +A caller that has to read or write a stream off the GPU path used to stage it through the device, and that fails two ways. +Threads reading fields back each allocate a field's worth of device memory beside whatever the card already holds, and double +fields of 512 MiB stopped with out of memory. And a process that sees no GPU cannot code at all. `compressHost` and `decompressHost` write and read the same stream as the kernels, in plain C++, one thread per call, with stack buffers only. `PfplHostCase` holds them to the kernels byte for byte on both kernel sets, both types and every size class. Its @@ -1902,7 +1902,7 @@ buffers only. `PfplHostCase` holds them to the kernels byte for byte on both ker conversion is undefined, so the host code tests the range first, and both then store the raw bits. `-ffp-contract` matters as it does in the kernels (see `fmulRounded`), so the rebuilt value goes through a `volatile`. -Measured on one Frontier login core, with no GPU in the process, on a 128 MiB polar tile: +Measured on one Frontier login core, with no GPU in the process, on a 128 MiB polar field: | Type, bound | Ratio | Encode | Decode | |---|---:|---:|---:|