diff --git a/CHANGELOG.md b/CHANGELOG.md index 964936da..a74289db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,19 @@ 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. +- `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 +- 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. - 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/CMakeLists.txt b/CMakeLists.txt index 1f863647..bfe30451 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 @@ -228,6 +320,7 @@ target_link_libraries(fzgmod_modules fzgmod_compile_settings ${FZGMOD_GPU_RUNTIME} PRIVATE + fzgmod_mem fzgmod_utils ) add_subdirectory(modules) @@ -456,6 +549,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/docs/codebase_notes.md b/docs/codebase_notes.md index f25c9a05..32ef0d8c 100644 --- a/docs/codebase_notes.md +++ b/docs/codebase_notes.md @@ -1816,3 +1816,96 @@ 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 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 | +|---|---:|---:|---:|---:| +| 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. + +## 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 field (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 | + +## 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, 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 +`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 field: + +| 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/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..e04ddb77 --- /dev/null +++ b/include/backend/fp.h @@ -0,0 +1,43 @@ +#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 +} + +/** 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/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..f2dfb5db --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream.cu @@ -0,0 +1,666 @@ +/** + * @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 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. +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; +} + +// 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) + 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] = 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] = magic; + w[1] = static_cast(padded); + 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, 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] = magic; + 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); +} + +// 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; + 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]; + 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); +} + +// 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" + +// FZ_PFPL_HARNESS=1 runs the generic harness kernels instead of the AMD ones (same stream; for comparison). +bool useHarness() { + static const bool harness = [] { + const char* e = std::getenv("FZ_PFPL_HARNESS"); + return e != nullptr && e[0] == '1'; + }(); + return harness; +} +#endif + +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 n < (static_cast(1) << 32) / sizeof(Real) && + kStreamHeader + kLcHeader + static_cast(chunksOf(n)) * (4u + kChunkBytes) + < (static_cast(1) << 32); +} + +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) { + 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.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; +} + +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) { 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; + 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))); + 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; + } + 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()); + } +}; + +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"); + 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"); + + // 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; + return qp; +} + +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. + *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; + 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); +} + +void decompress(Plan& plan, const void* d_stream, float* d_out, size_t n, fz::stream_t stream) { + 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); +} + +#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 new file mode 100644 index 00000000..3a26791b --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream.h @@ -0,0 +1,145 @@ +#pragma once + +/** + * @file modules/fused/chunk_fusion/pfpl_stream.h + * @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 + * 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. + * + * 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. The stream must be 4-byte aligned and the + * field aligned for its type. + */ + +#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; ///< 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: 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 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 + * 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. + */ +struct PlanAccess; + +class Plan { +public: + Plan() = default; + ~Plan(); + Plan(const Plan&) = delete; + Plan& operator=(const Plan&) = delete; + + /** 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); + + /** The same, for a field of `fieldBytes` of either type: the state depends only on its chunks. */ + void reserveBytes(size_t fieldBytes); + +private: + 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 + 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) + 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 +}; + +/** + * 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); + +/** + * 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); + +/** + * 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_amd.inl b/modules/fused/chunk_fusion/pfpl_stream_amd.inl new file mode 100644 index 00000000..0dded014 --- /dev/null +++ b/modules/fused/chunk_fusion/pfpl_stream_amd.inl @@ -0,0 +1,496 @@ +// 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, 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); +// - 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. +// 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 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 && 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; } + +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)); + } +} + +template +__global__ void __launch_bounds__(kTPB) +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) { + 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__ __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 + + 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) * L::kPerChunk; + + // 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; + 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) lastLane[wave * kPerT + k] = c[k]; + } + __syncthreads(); + #pragma unroll + for (int k = 0; k < kPerT; ++k) { + 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. + uint32_t v[kWordsT]; + uint32_t bmp = 0; + { + const uint32_t* row = planeRow(buf, t); + #pragma unroll + 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); + } + } + 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[kWordsT]; + #pragma unroll + 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 < kWordsT; ++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 < 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)); + // 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 < 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)); + 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] = chunks * static_cast(kChunkBytes); + w[5] = chunks; + } + if (kStreamHeader <= capacity) { + w[0] = Elem::kMagic; + w[1] = static_cast(padded); + w[2] = __float_as_uint(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 < 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[kWordsT * lane + r] = word; + } + return rpos; +} + +template +__global__ void __launch_bounds__(kTPB) +decodeKernel(const byte* __restrict__ stream, Real* __restrict__ out, uint32_t n, + unsigned long long* __restrict__ state, uint32_t epoch) { + 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__ __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; + 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]; + 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[kWordsT]; + 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 < kWordsT; ++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 = const_cast(planeRow(buf, t)); + #pragma unroll + 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. + Word d[kPerT]; + #pragma unroll + 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 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) waveTotal[k * 8 + wave] = d[k]; + } + __syncthreads(); + 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) { + Word before = 0, slab = 0; + #pragma unroll + for (int ww = 0; ww < 8; ++ww) { + const Word x = waveTotal[k * 8 + ww]; + before += ww < wave ? x : Word(0); + slab += x; + } + const Word code = d[k] + before + carry; + carry += slab; + const uint32_t e = t + kTPB * k; + if (e < cnt) out[base + e] = valueOfCode(code, ebx2); + } +} + +} // namespace amd 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/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/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 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..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 @@ -189,12 +205,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( diff --git a/tests/stages/test_pfpl_stream.cpp b/tests/stages/test_pfpl_stream.cpp new file mode 100644 index 00000000..c1b2df35 --- /dev/null +++ b/tests/stages/test_pfpl_stream.cpp @@ -0,0 +1,892 @@ +/** + * 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 + * + * 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. + */ + +#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 +#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); + } +} + +// ─── 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); +} + +// ─── 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); +}