device_gzip.cuh

Go to the documentation of this file.

#pragma once
 
#include <cuda.h>
#include <cudaTypedefs.h>
#include <nvcomp/deflate.h>
#include <nvcomp/gzip.h>
 
#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <condition_variable>
#include <cstdint>
#include <cstring>
#include <filesystem>
#include <functional>
#include <limits>
#include <mutex>
#include <optional>
#include <span>
#include <stdexcept>
#include <string>
#include <string_view>
#include <thread>
#include <vector>
 
#include <fcntl.h>
#include <sys/stat.h>
#include <unistd.h>
 
#include <cub/block/block_reduce.cuh>
#include <cub/device/device_copy.cuh>
#include <cub/device/device_select.cuh>
#include <cuda/algorithm>
#include <cuda/buffer>
#include <cuda/iterator>
#include <cuda/memory_pool>
#include <cuda/stream>
 
#if defined(__x86_64__) || defined(_M_X64)
    #include <immintrin.h>
#elif defined(__aarch64__) || defined(_M_ARM64)
    #include <arm_neon.h>
#endif
 
#include <cuddl/detail/fastx_sequence_file.hpp>
#include <cuddl/detail/pinned_pool.cuh>
#include <cuddl/detail/worker_pool.hpp>
#include <cuddl/device_span.cuh>
#include <cuddl/error.hpp>
 
namespace cuddl::detail {
 
struct gzip_file_probe {
    uint64_t compressed = 0;  // file size
    // Inflated size and CRC32 of the whole file. For one member, its trailer's ISIZE and CRC32;
    // for BGZF, the sum and combination over its blocks. A plain file's size, without a CRC.
    uint32_t isize = 0;
    uint32_t crc = 0;
    uint32_t members = 0;  // gzip members the device inflates: 1, or BGZF's non-empty blocks
    bool device = false;   // a candidate the device may inflate, or plain FASTA
    bool plain = false;    // uncompressed FASTA, which the device normalises as read
    bool bgzf = false;     // blocked gzip, with each block's size in its header
 
    [[nodiscard]] uint64_t host_bytes() const noexcept {
        return plain ? uint64_t{isize} + 1 : compressed;
    }
};
 
[[nodiscard]] inline uint32_t bgzf_block_size(unsigned char const* header) noexcept {
    unsigned char const expected[16] = {0x1f, 0x8b, 8, 4, 0, 0, 0, 0, 0, 0, 6, 0, 'B', 'C', 2, 0};
    for (size_t i = 0; i < sizeof(expected); ++i) {
        // Bytes 4-9 are MTIME, XFL and OS, which may hold anything.
        if ((i < 4 || i > 9) && header[i] != expected[i]) return 0;
    }
    return (uint32_t{header[16]} | uint32_t{header[17]} << 8U) + 1U;
}
 
[[nodiscard]] inline uint32_t crc32_combine(uint32_t a, uint32_t b, uint64_t length) noexcept {
    constexpr uint32_t polynomial = 0xedb88320U;
    auto const multiply = [](uint32_t x, uint32_t y) {
        uint32_t product = 0;
        for (uint32_t bit = 1U << 31U; bit != 0; bit >>= 1U) {
            if ((x & bit) != 0) product ^= y;
            y = (y & 1U) != 0 ? (y >> 1U) ^ polynomial : y >> 1U;
        }
        return product;
    };
    // power = x^(2^k) for bit k of 8 * length; x^1 is 1 << 30 in reflected order.
    uint32_t shift = 1U << 31U, power = 1U << 30U;
    for (uint64_t bits = length * 8; bits != 0; bits >>= 1U) {
        if ((bits & 1U) != 0) shift = multiply(power, shift);
        power = multiply(power, power);
    }
    return multiply(shift, a) ^ b;
}
 
template <typename Block>
[[nodiscard]] std::optional<std::pair<uint32_t, uint64_t>>
walk_bgzf(unsigned char const* bytes, uint64_t size, Block&& block) {
    uint64_t at = 0, total = 0;
    uint32_t count = 0;
    while (at < size) {
        auto const block_size = size - at >= 18 ? bgzf_block_size(bytes + at) : 0;
        if (block_size < 26 || block_size > size - at) return std::nullopt;
        uint32_t inflated = 0;
        std::memcpy(&inflated, bytes + at + block_size - 4, 4);
        if (inflated != 0) {
            block(at + 18, uint64_t{block_size} - 26, inflated);
            total += inflated;
            ++count;
        }
        at += block_size;
    }
    return std::pair{count, total};
}
 
[[nodiscard]] inline gzip_file_probe probe_gzip_file(std::filesystem::path const& path) {
    gzip_file_probe probe;
    int const fd = ::open(path.c_str(), O_RDONLY);
    if (fd == -1) return probe;
    struct stat status{};
    unsigned char header[18] = {};
    unsigned char trailer[8 + 18] = {};  // a member's trailer, and the next BGZF header
    if (::fstat(fd, &status) == 0 && status.st_size >= 4 &&
        ::pread(fd, header, sizeof(header), 0) >= 4) {
        auto const size = static_cast<uint64_t>(status.st_size);
        probe.compressed = size;
        auto name = path.filename().string();
        std::transform(name.begin(), name.end(), name.begin(), [](unsigned char c) {
            return static_cast<char>(std::tolower(c));
        });
        auto const fastq =
            name.find(".fq") != std::string::npos || name.find(".fastq") != std::string::npos;
        if (header[0] == '>') {
            // Slot offsets count a plain file's bytes and its guard byte in 32 bits.
            probe.plain = size < std::numeric_limits<uint32_t>::max();
            probe.isize = static_cast<uint32_t>(probe.plain ? size : 0);
            probe.device = probe.plain && !fastq;
        } else if (size >= 26 && bgzf_block_size(header) != 0) {
            // One read per block: its trailer, and the header of the block after it.
            uint64_t at = 0, total = 0;
            uint32_t crc = 0, block_size = bgzf_block_size(header);
            bool tiled = false;
            while (block_size >= 26 && block_size <= size - at) {
                auto const next = at + block_size;
                auto const wanted = next == size ? 8 : std::min<uint64_t>(26, size - next + 8);
                if (::pread(fd, trailer, wanted, static_cast<off_t>(next - 8)) !=
                    static_cast<ssize_t>(wanted)) {
                    break;
                }
                uint32_t block_crc = 0, inflated = 0;
                std::memcpy(&block_crc, trailer, 4);
                std::memcpy(&inflated, trailer + 4, 4);
                if (inflated != 0) {
                    crc = crc32_combine(crc, block_crc, inflated);
                    total += inflated;
                    ++probe.members;
                }
                if (next == size) {
                    tiled = true;
                    break;
                }
                block_size = wanted == 26 ? bgzf_block_size(trailer + 8) : 0;
                at = next;
            }
            probe.bgzf = true;
            probe.crc = crc;
            probe.isize = static_cast<uint32_t>(total);
            probe.device =
                tiled && total != 0 && total < std::numeric_limits<uint32_t>::max() && !fastq;
        } else if (size >= 18 && ::pread(fd, trailer, 8, static_cast<off_t>(size - 8)) == 8) {
            std::memcpy(&probe.crc, trailer, 4);
            std::memcpy(&probe.isize, trailer + 4, 4);
            probe.members = 1;
            probe.device = header[0] == 0x1f && header[1] == 0x8b && header[2] == 8 &&
                           (header[3] & 0xe4U) == 0 && probe.isize != 0 && !fastq;
        }
    }
    ::close(fd);
    return probe;
}
 
template <typename Work>
void parallel_for(size_t count, size_t workers, Work&& work) {
    std::atomic<size_t> next{0};
    auto run = [&] {
        for (size_t index = next++; index < count; index = next++) work(index);
    };
    std::vector<std::jthread> threads;
    auto const helpers = std::min(workers, count) > 0 ? std::min(workers, count) - 1 : 0;
    threads.reserve(helpers);
    for (size_t thread = 0; thread < helpers; ++thread) threads.emplace_back(run);
    run();
}
 
[[nodiscard]] inline bool has_gzip_signature_portable(char const* data, size_t size) {
    for (size_t at = 0; size - at >= 4;) {
        auto const* const found =
            static_cast<char const*>(::memmem(data + at, size - at - 1, "\x1f\x8b\x08", 3));
        if (found == nullptr) return false;
        if ((static_cast<uint8_t>(found[3]) & 0xe0U) == 0) return true;
        at = static_cast<size_t>(found - data) + 1;
    }
    return false;
}
 
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
__attribute__((target("avx512f,avx512bw"))) inline bool
has_gzip_signature_avx512(char const* data, size_t size) {
    size_t at = 0;
    for (; at + 67 <= size; at += 64) {
        auto const first =
            _mm512_cmpeq_epi8_mask(_mm512_loadu_si512(data + at), _mm512_set1_epi8(0x1f));
        if (first == 0) continue;
        auto const second = _mm512_cmpeq_epi8_mask(
            _mm512_loadu_si512(data + at + 1), _mm512_set1_epi8(static_cast<char>(0x8b))
        );
        auto const third =
            _mm512_cmpeq_epi8_mask(_mm512_loadu_si512(data + at + 2), _mm512_set1_epi8(0x08));
        auto const flags = _mm512_testn_epi8_mask(
            _mm512_loadu_si512(data + at + 3), _mm512_set1_epi8(static_cast<char>(0xe0))
        );
        if ((first & second & third & flags) != 0) return true;
    }
    return has_gzip_signature_portable(data + at, size - at);
}
 
__attribute__((target("avx2"))) inline bool has_gzip_signature_avx2(char const* data, size_t size) {
    size_t at = 0;
    for (; at + 35 <= size; at += 32) {
        auto const* const block = reinterpret_cast<__m256i const*>(data + at);
        auto const first = _mm256_cmpeq_epi8(_mm256_loadu_si256(block), _mm256_set1_epi8(0x1f));
        if (_mm256_testz_si256(first, first)) continue;
        auto const second = _mm256_cmpeq_epi8(
            _mm256_loadu_si256(reinterpret_cast<__m256i const*>(data + at + 1)),
            _mm256_set1_epi8(static_cast<char>(0x8b))
        );
        auto const third = _mm256_cmpeq_epi8(
            _mm256_loadu_si256(reinterpret_cast<__m256i const*>(data + at + 2)),
            _mm256_set1_epi8(0x08)
        );
        auto const flags = _mm256_cmpeq_epi8(
            _mm256_and_si256(
                _mm256_loadu_si256(reinterpret_cast<__m256i const*>(data + at + 3)),
                _mm256_set1_epi8(static_cast<char>(0xe0))
            ),
            _mm256_setzero_si256()
        );
        auto const all =
            _mm256_and_si256(_mm256_and_si256(first, flags), _mm256_and_si256(second, third));
        if (!_mm256_testz_si256(all, all)) return true;
    }
    return has_gzip_signature_portable(data + at, size - at);
}
#endif
 
#if defined(__aarch64__) || defined(_M_ARM64)
inline bool has_gzip_signature_neon(char const* data, size_t size) {
    size_t at = 0;
    auto const* bytes = reinterpret_cast<uint8_t const*>(data);
    for (; at + 19 <= size; at += 16) {
        auto const first = vceqq_u8(vld1q_u8(bytes + at), vdupq_n_u8(0x1f));
        if (vmaxvq_u8(first) == 0) continue;
        auto const second = vceqq_u8(vld1q_u8(bytes + at + 1), vdupq_n_u8(0x8b));
        auto const third = vceqq_u8(vld1q_u8(bytes + at + 2), vdupq_n_u8(8));
        auto const flags =
            vceqq_u8(vandq_u8(vld1q_u8(bytes + at + 3), vdupq_n_u8(0xe0)), vdupq_n_u8(0));
        if (vmaxvq_u8(vandq_u8(vandq_u8(first, flags), vandq_u8(second, third))) != 0) {
            return true;
        }
    }
    return has_gzip_signature_portable(data + at, size - at);
}
#endif
 
[[nodiscard]] inline auto resolve_gzip_signature_scanner() -> bool (*)(char const*, size_t) {
    if (ablated("simd-parsing")) return has_gzip_signature_portable;
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
    if (__builtin_cpu_supports("avx512bw")) return has_gzip_signature_avx512;
    if (__builtin_cpu_supports("avx2")) return has_gzip_signature_avx2;
#endif
#if defined(__aarch64__) || defined(_M_ARM64)
    return has_gzip_signature_neon;
#endif
    return has_gzip_signature_portable;
}
 
static bool (*const has_gzip_signature)(char const*, size_t) = resolve_gzip_signature_scanner();
 
namespace device_gzip {
 
inline constexpr uint32_t crc_polynomial = 0xEDB88320U;
inline constexpr uint32_t block_size = 256;
 
__device__ __forceinline__ bool is_whitespace(char c) {
    return c == '\n' || c == '\r' || c == ' ' || c == '\t';
}
 
struct kept_byte {
    __device__ bool operator()(char c) const {
        return !is_whitespace(c);
    }
};
 
__device__ inline uint32_t crc_multiply(uint32_t a, uint32_t b) {
    uint32_t m = 1U << 31U, p = 0;
    while (true) {
        if ((a & m) != 0) {
            p ^= b;
            if ((a & (m - 1)) == 0) break;
        }
        m >>= 1U;
        b = (b & 1U) != 0 ? (b >> 1U) ^ crc_polynomial : b >> 1U;
    }
    return p;
}
 
__device__ inline uint32_t
crc_combine(uint32_t first, uint32_t second, uint64_t second_bytes, uint32_t const* x2n) {
    uint32_t shift = 1U << 31U;
    for (uint32_t k = 3; second_bytes != 0; second_bytes >>= 1U, ++k) {
        if ((second_bytes & 1U) != 0) shift = crc_multiply(x2n[k & 31U], shift);
    }
    return crc_multiply(shift, first) ^ second;
}
 
static __global__ void slot_init_kernel(
    char* buffer,
    uint64_t const* slots,
    uint32_t files,
    unsigned long long* first_header,
    unsigned long long* header_count
) {
    auto const file = blockIdx.x * blockDim.x + threadIdx.x;
    if (file == 0) *header_count = 0;
    if (file >= files) return;
    buffer[slots[file]] = '\n';
    first_header[file] = slots[file + 1];
}
 
static __global__ void slot_crc_kernel(char const* buffer, uint64_t const* slots, uint32_t* crcs) {
    __shared__ uint32_t table[4][256];
    __shared__ uint32_t x2n[32];
    __shared__ uint32_t crc[block_size];
    __shared__ uint64_t bytes[block_size];
    auto const t = threadIdx.x;
    uint32_t entry = t;
    _Pragma("unroll")
    for (int bit = 0; bit < 8; ++bit) {
        entry = (entry & 1U) != 0 ? (entry >> 1U) ^ crc_polynomial : entry >> 1U;
    }
    table[0][t] = entry;
    __syncthreads();
    for (int slice = 1; slice < 4; ++slice) {
        auto const previous = table[slice - 1][t];
        table[slice][t] = (previous >> 8U) ^ table[0][previous & 0xFFU];
        __syncthreads();
    }
    if (t == 0) {
        uint32_t power = 1U << 30U;
        x2n[0] = power;
        for (int n = 1; n < 32; ++n) x2n[n] = power = crc_multiply(power, power);
    }
    auto const* const data = reinterpret_cast<unsigned char const*>(buffer + slots[blockIdx.x] + 1);
    uint64_t const size = slots[blockIdx.x + 1] - slots[blockIdx.x] - 1;
    uint64_t const per = ((size + block_size - 1) / block_size + 15) & ~uint64_t{15};
    uint64_t const begin = cuda::std::min<uint64_t>(size, t * per);
    uint64_t const end = cuda::std::min<uint64_t>(size, begin + per);
    uint32_t c = 0xFFFFFFFFU;
    auto byte = [&](unsigned char value) {
        c = table[0][(c ^ value) & 0xFFU] ^ (c >> 8U);
    };
    uint64_t at = begin;
    for (; at < end && (reinterpret_cast<uintptr_t>(data + at) & 15U) != 0; ++at) byte(data[at]);
    for (; at + 16 <= end; at += 16) {
        auto const words = *reinterpret_cast<uint4 const*>(data + at);
        uint32_t const lanes[4] = {words.x, words.y, words.z, words.w};
        _Pragma("unroll")
        for (uint32_t word : lanes) {
            c ^= word;
            c = table[3][c & 0xFFU] ^ table[2][(c >> 8U) & 0xFFU] ^ table[1][(c >> 16U) & 0xFFU] ^
                table[0][c >> 24U];
        }
    }
    for (; at < end; ++at) byte(data[at]);
    crc[t] = ~c;
    bytes[t] = end - begin;
    __syncthreads();
    for (uint32_t stride = 1; stride < block_size; stride *= 2) {
        if (t % (2 * stride) == 0) {
            crc[t] = crc_combine(crc[t], crc[t + stride], bytes[t + stride], x2n);
            bytes[t] += bytes[t + stride];
        }
        __syncthreads();
    }
    if (t == 0) crcs[blockIdx.x] = crc[0];
}
 
static __global__ void header_kernel(
    char const* buffer,
    uint64_t size,
    uint64_t* headers,
    uint32_t capacity,
    unsigned long long* count
) {
    constexpr uint64_t ones = 0x0101010101010101ULL;
    auto const has_marker = [](uint64_t word) {
        auto const x = word ^ (ones * '>');
        return ((x - ones) & ~x & 0x8080808080808080ULL) != 0;
    };
    uint64_t const words = (size + 15) / 16;
    for (uint64_t w = blockIdx.x * uint64_t{blockDim.x} + threadIdx.x; w < words;
         w += uint64_t{gridDim.x} * blockDim.x) {
        auto const value = *reinterpret_cast<ulonglong2 const*>(buffer + w * 16);
        if (!has_marker(value.x) && !has_marker(value.y)) continue;
        for (uint32_t j = 0; j < 16; ++j) {
            auto const position = w * 16 + j;
            if (position == 0 || position >= size || buffer[position] != '>') continue;
            auto const previous = buffer[position - 1];
            if (previous != '\n' && previous != '\r') continue;
            auto const slot = atomicAdd(count, 1ULL);
            if (slot < capacity) headers[slot] = position;
        }
    }
}
 
__device__ __forceinline__ uint32_t
slot_containing(uint64_t const* slots, uint32_t files, uint64_t position) {
    uint32_t low = 0, high = files;
    while (high - low > 1) {
        auto const middle = (low + high) / 2;
        if (slots[middle] <= position) {
            low = middle;
        } else {
            high = middle;
        }
    }
    return low;
}
 
static __global__ void blank_header_kernel(
    char* buffer,
    uint64_t const* headers,
    unsigned long long const* count,
    uint32_t capacity,
    uint64_t const* slots,
    uint32_t files,
    unsigned long long* first_header
) {
    auto const lane = threadIdx.x % 32U;
    auto const total = cuda::std::min<unsigned long long>(*count, capacity);
    auto const warps = uint64_t{gridDim.x} * blockDim.x / 32U;
    for (uint64_t h = (blockIdx.x * uint64_t{blockDim.x} + threadIdx.x) / 32U; h < total;
         h += warps) {
        auto const position = headers[h];
        auto const low = slot_containing(slots, files, position);
        if (lane == 0) atomicMin(first_header + low, static_cast<unsigned long long>(position));
        auto const end = slots[low + 1];
        for (uint64_t at = position + 1; at < end; at += 32) {
            auto const index = at + lane;
            auto const c = index < end ? buffer[index] : '\n';
            auto const stops = __ballot_sync(0xFFFFFFFFU, c == '\n' || c == '\r');
            auto const stop = stops != 0 ? static_cast<uint32_t>(__ffs(stops) - 1) : 32U;
            if (lane < stop) buffer[index] = '\n';
            if (stops != 0) break;
        }
    }
}
 
__device__ __forceinline__ uint32_t kept_in_word(uint64_t word, uint32_t first, uint32_t last) {
    constexpr uint64_t low7 = 0x7F7F7F7F7F7F7F7FULL, ones = 0x0101010101010101ULL;
    // The high bit of every zero byte of x, exactly: the addition carries within bytes only.
    auto const zero_bytes = [](uint64_t x) {
        return ~(((x & low7) + low7) | x | low7);
    };
    auto const blank = zero_bytes(word ^ (ones * '\n')) | zero_bytes(word ^ (ones * '\r')) |
                       zero_bytes(word ^ (ones * ' ')) | zero_bytes(word ^ (ones * '\t'));
    auto const width = last - first;
    auto const span = (width == 8 ? ~uint64_t{0} : (uint64_t{1} << (8 * width)) - 1) << (8 * first);
    return static_cast<uint32_t>(__popcll(~blank & span & ~low7));
}
 
inline constexpr uint64_t kept_tile = uint64_t{block_size} * 16 * 8;
 
static __global__ void kept_count_kernel(
    char const* buffer,
    uint64_t const* slots,
    uint32_t files,
    unsigned long long* kept
) {
    using block_reduce = cub::BlockReduce<uint32_t, block_size>;
    __shared__ typename block_reduce::TempStorage reduce_storage;
    __shared__ uint32_t tile_slot;
    constexpr uint32_t straddles = ~uint32_t{0};
    auto const begin = slots[0], end = slots[files];
    auto const last_tile = (end + kept_tile - 1) / kept_tile;
    for (auto tile = begin / kept_tile + blockIdx.x; tile < last_tile; tile += gridDim.x) {
        auto const low = cuda::std::max(begin, tile * kept_tile);
        auto const high = cuda::std::min(end, (tile + 1) * kept_tile);
        if (threadIdx.x == 0) {
            auto const slot = slot_containing(slots, files, low);
            tile_slot = slots[slot + 1] >= high ? slot : straddles;
        }
        __syncthreads();
        auto const slot = tile_slot;
        uint32_t count = 0;
        for (auto at = tile * kept_tile + uint64_t{threadIdx.x} * 16; at < (tile + 1) * kept_tile;
             at += uint64_t{block_size} * 16) {
            if (at + 16 <= low || at >= high) continue;
            auto const words = *reinterpret_cast<ulonglong2 const*>(buffer + at);
            uint64_t const halves[2] = {words.x, words.y};
            _Pragma("unroll")
            for (uint32_t half = 0; half < 2; ++half) {
                auto const base = at + 8 * half;
                auto const first =
                    static_cast<uint32_t>(low > base ? cuda::std::min<uint64_t>(low - base, 8) : 0);
                auto const last = static_cast<uint32_t>(
                    high > base ? cuda::std::min<uint64_t>(high - base, 8) : 0
                );
                if (first >= last) continue;
                if (slot != straddles) {
                    count += kept_in_word(halves[half], first, last);
                    continue;
                }
                auto const head = slot_containing(slots, files, base + first);
                if (head == slot_containing(slots, files, base + last - 1)) {
                    atomicAdd(kept + head, kept_in_word(halves[half], first, last));
                    continue;
                }
                for (auto byte = first; byte < last; ++byte) {
                    atomicAdd(
                        kept + slot_containing(slots, files, base + byte),
                        kept_in_word(halves[half], byte, byte + 1)
                    );
                }
            }
        }
        if (slot != straddles) {
            auto const sum = block_reduce(reduce_storage).Sum(count);
            if (threadIdx.x == 0) atomicAdd(kept + slot, sum);
        }
        __syncthreads();
    }
}
 
static __global__ void slot_finish_kernel(
    char* buffer,
    uint64_t const* slots,
    unsigned long long const* first_header,
    char* first_char
) {
    auto const begin = slots[blockIdx.x] + 1;
    auto const end = slots[blockIdx.x + 1];
    if (threadIdx.x == 0) {
        auto at = begin;
        while (at < end && (buffer[at] == '\n' || buffer[at] == '\r')) ++at;
        first_char[blockIdx.x] = at < end ? buffer[at] : '\0';
    }
    __syncthreads();
    auto const head = cuda::std::min<uint64_t>(first_header[blockIdx.x], end);
    for (auto at = begin + threadIdx.x; at < head; at += block_size) buffer[at] = '\n';
}
 
inline size_t hardware_inflate_limit(int device) {
    void* address = nullptr;
    cudaDriverEntryPointQueryResult found{};
    if (cudaGetDriverEntryPointByVersion(
            "cuDeviceGetAttribute", &address, 12000U, cudaEnableDefault, &found
        ) != cudaSuccess ||
        found != cudaDriverEntryPointSuccess) {
        return 0;
    }
    int limit = 0;
    auto const attribute = reinterpret_cast<decltype(&cuDeviceGetAttribute)>(address);
    auto const status =
        attribute(&limit, CU_DEVICE_ATTRIBUTE_MEM_DECOMPRESS_MAXIMUM_LENGTH, device);
    return status == CUDA_SUCCESS && limit > 0 ? static_cast<size_t>(limit) : 0;
}
 
[[nodiscard]] inline cuda::device_memory_pool_ref hardware_inflate_pool(cuda::device_ref device) {
    static std::mutex mutex;
    static std::vector<cudaMemPool_t> pools(cuda::devices.size());
    std::lock_guard lock(mutex);
    auto& pool = pools.at(static_cast<size_t>(device.get()));
    if (pool == nullptr) {
        cudaMemPoolProps properties{};
        properties.allocType = cudaMemAllocationTypePinned;
        properties.location = {cudaMemLocationTypeDevice, device.get()};
        properties.usage = cudaMemPoolCreateUsageHwDecompress;
        CUDDL_CUDA_CALL(cudaMemPoolCreate(&pool, &properties));
    }
    return cuda::device_memory_pool_ref(pool);
}
 
inline PFN_cuMemBatchDecompressAsync_v12060 batch_decompress_entry() {
    void* address = nullptr;
    cudaDriverEntryPointQueryResult found{};
    if (cudaGetDriverEntryPointByVersion(
            "cuMemBatchDecompressAsync", &address, 12060U, cudaEnableDefault, &found
        ) != cudaSuccess ||
        found != cudaDriverEntryPointSuccess) {
        return nullptr;
    }
    return reinterpret_cast<PFN_cuMemBatchDecompressAsync_v12060>(address);
}
 
inline size_t deflate_offset(char const* bytes, size_t size) {
    auto const flags = static_cast<unsigned char>(bytes[3]);
    size_t at = 10;
    // FNAME, then FCOMMENT: zero-terminated strings.
    for (auto const flag : {8U, 16U}) {
        if ((flags & flag) == 0) continue;
        auto const* const end =
            static_cast<char const*>(std::memchr(bytes + at, '\0', size - std::min(at, size)));
        if (end == nullptr) return 0;
        at = static_cast<size_t>(end - bytes) + 1;
    }
    if ((flags & 2) != 0) at += 2;  // FHCRC
    return at + 8 < size ? at : 0;
}
 
}  // namespace device_gzip
 
struct inflated_genome {
    size_t id;
    char const* bases;
    size_t size;
};
 
struct decoded_fasta {
    size_t id;
    std::string_view bytes;
};
 
class device_fasta_pipeline {
    struct lane_state;
 
   public:
    static constexpr size_t lane_count = 6;
    static constexpr size_t max_files = 4096;
    static constexpr size_t max_members = 65536;
 
    device_fasta_pipeline(
        cuda::stream_ref stream,
        size_t slot_capacity,
        size_t compressed_capacity,
        size_t workers,
        bool hardware_inflate
    )
        : workers_(std::max<size_t>(1, workers)),
          compressed_capacity_(compressed_capacity),
          slot_capacity_(slot_capacity),
          header_capacity_(
              static_cast<uint32_t>(std::clamp<size_t>(
                  slot_capacity / 4096,
                  1 << 16,
                  std::numeric_limits<uint32_t>::max()
              ))
          ) {
        auto const options = nvcompBatchedDeflateDecompressDefaultOpts;
        nvcomp_temp_bytes_ = 0;
        if (compressed_capacity_ != 0 &&
            nvcompBatchedDeflateDecompressGetTempSizeAsync(
                max_members,
                std::min<size_t>(slot_capacity_, std::numeric_limits<uint32_t>::max()),
                options,
                &nvcomp_temp_bytes_,
                slot_capacity_
            ) != nvcompSuccess) {
            throw std::bad_alloc();
        }
        select_temp_bytes_ = 0;
        auto const selected = compressed_capacity_ == 0 ? cub::DeviceSelect::If(
                                                              nullptr,
                                                              select_temp_bytes_,
                                                              static_cast<char*>(nullptr),
                                                              static_cast<int64_t*>(nullptr),
                                                              static_cast<int64_t>(slot_capacity_),
                                                              device_gzip::kept_byte{},
                                                              stream.get()
                                                          )
                                                        : cub::DeviceSelect::If(
                                                              nullptr,
                                                              select_temp_bytes_,
                                                              static_cast<char*>(nullptr),
                                                              static_cast<char*>(nullptr),
                                                              static_cast<int64_t*>(nullptr),
                                                              static_cast<int64_t>(slot_capacity_),
                                                              device_gzip::kept_byte{},
                                                              stream.get()
                                                          );
        if (selected != cudaSuccess) throw std::runtime_error(cudaGetErrorString(selected));
        auto const device = stream.device();
        if (hardware_inflate) {
            batch_decompress_ = device_gzip::batch_decompress_entry();
            if (batch_decompress_ != nullptr) {
                hardware_limit_ = device_gzip::hardware_inflate_limit(device.get());
            }
        }
        streaming_window_ =
            compressed_capacity_ == 0
                ? std::min(
                      slot_capacity_,
                      static_cast<size_t>(
                          device.attribute(cuda::device_attributes::max_access_policy_window_size)
                      )
                  )
                : size_t{0};
        for (auto& lane : lanes_) {
            lane.stream.emplace(device);
            auto const s = cuda::stream_ref{*lane.stream};
            lane.uploaded.emplace(s);
            lane.classified.emplace(s);
            lane.released.emplace(stream);
        }
        // The other lanes allocate on first use, so the device starts on the first batch while
        // they page-lock their compressed buffers.
        allocate(lanes_[0]);
    }
 
    [[nodiscard]] size_t slot_capacity() const noexcept {
        return slot_capacity_;
    }
 
    [[nodiscard]] size_t compressed_capacity() const noexcept {
        return compressed_capacity_;
    }
 
    [[nodiscard]] size_t hardware_limit() const noexcept {
        return hardware_limit_;
    }
 
    [[nodiscard]] static size_t lane_overhead(size_t slot_capacity) noexcept {
        return slot_capacity / 4096 * sizeof(uint64_t) + (size_t{64} << 20);
    }
 
    [[nodiscard]] size_t batch_files(
        std::span<size_t const> ids,
        std::span<gzip_file_probe const> probes
    ) const noexcept {
        size_t files = 0, compressed = 0, slots = 0, members = 0;
        while (files < ids.size() && files < max_files) {
            auto const& probe = probes[ids[files]];
            if (compressed + probe.host_bytes() > compressed_capacity_ ||
                slots + probe.isize + 1 > slot_capacity_ || members + probe.members > max_members) {
                break;
            }
            compressed += probe.host_bytes();
            slots += probe.isize + 1;
            members += probe.members;
            ++files;
        }
        return files;
    }
 
    [[nodiscard]] Result<size_t> submit(
        size_t lane_index,
        std::span<size_t const> ids,
        std::span<std::filesystem::path const> paths,
        std::span<gzip_file_probe const> probes,
        bool direct_reads
    ) {
        auto& lane = lanes_[lane_index];
        CUDDL_TRY(cuda_try([&] { allocate(lane); }));
        lane.decoded = false;
        if (ids.empty()) return Err(Error::invalid_argument("an inflate batch needs a file"));
        auto const files = batch_files(ids, probes);
        // The pinned bytes are rewritten below, so the previous batch's upload must be done.
        CUDDL_CUDA_TRY(lane.uploaded->sync());
        auto table = host_table(lane);
        // Plain files lie in the page-locked buffer exactly as they will in the slots, which
        // they lead, so one copy places them all. Compressed bytes follow them.
        uint64_t plain_at = 0;
        for (size_t file = 0; file < files; ++file) {
            auto const& probe = probes[ids[file]];
            if (probe.plain) plain_at += probe.host_bytes();
        }
        uint64_t compressed_at = plain_at, slot_at = 0;
        for (size_t file = 0; file < files; ++file) {
            auto const& probe = probes[ids[file]];
            table.compressed_bytes[file] = probe.compressed;
            if (probe.plain) {
                table.compressed_offsets[file] = slot_at + 1;
                slot_at += probe.host_bytes();
            } else {
                table.compressed_offsets[file] = compressed_at;
                compressed_at += probe.compressed;
            }
        }
        std::atomic<bool> failed{false};
        readers_.for_each(files, [&](size_t file) {
            auto const* const path = paths[ids[file]].c_str();
            auto* const out = lane.pinned->data() + table.compressed_offsets[file];
            auto const bytes = table.compressed_bytes[file];
            int const fd = ::open(path, O_RDONLY);
            // Files lie back to back in the lane, so a direct read goes through a bounce buffer.
            bool const read = fd != -1 && read_whole_file(fd, out, bytes, direct_reads);
            auto const done = read ? bytes : 0;
            if (!read) failed = true;
            // A matching trailer cannot distinguish repeated first and last members.
            // A signature inside compressed data is ambiguous too, so use the host parser.
            // BGZF names every member in its block headers.
            auto const& probe = probes[ids[file]];
            lane.possible_members[file] =
                !probe.plain && !probe.bgzf && done > 10 && has_gzip_signature(out + 10, done - 10);
        });
        if (failed) return Err(Error::invalid_argument("cannot read reference inputs"));
        lane.files.clear();
        lane.fallback.clear();
        lane.member_first.clear();
        engine_params_.clear();
        slot_at = 0;
        uint32_t members = 0;
        auto const accept = [&](size_t file) {
            table.slots[lane.files.size()] = slot_at;
            lane.member_first.push_back(members);
            slot_at += probes[ids[file]].isize + 1;
            lane.files.push_back(ids[file]);
        };
        for (size_t file = 0; file < files; ++file) {
            if (probes[ids[file]].plain) accept(file);
        }
        auto const plain = lane.files.size();
        for (size_t file = 0; file < files; ++file) {
            if (!probes[ids[file]].plain && lane.possible_members[file]) {
                lane.fallback.push_back(ids[file]);
            }
        }
        // One member: raw DEFLATE at @p deflate in the page-locked copy, which the payload buffer
        // mirrors, inflating to @p size bytes at @p out.
        auto const add_member = [&](unsigned char const* deflate,
                                    uint64_t deflate_bytes,
                                    uint32_t size,
                                    char* out) {
            auto* const source = lane.payload->data() +
                                 (reinterpret_cast<char const*>(deflate) - lane.pinned->data());
            table.member_sources[members] = source;
            table.member_source_bytes[members] = deflate_bytes;
            table.member_outputs[members] = out;
            table.member_output_bytes[members] = size;
            if (hardware_limit_ != 0 && size <= hardware_limit_) {
                CUmemDecompressParams params{};
                params.src = source;
                params.srcNumBytes = deflate_bytes;
                params.dst = out;
                params.dstNumBytes = size;
                // The low half of the member's 64-bit inflated size, zeroed before the batch.
                params.dstActBytes =
                    reinterpret_cast<cuuint32_t*>(device_tables(lane).member_inflated + members);
                params.algo = CU_MEM_DECOMPRESS_ALGORITHM_DEFLATE;
                engine_params_.push_back(params);
            }
            ++members;
        };
        // Gzip files whose members all fit the decompression engine come first, then those
        // nvCOMP inflates on the SMs, so each group's members form one range. Both take raw
        // DEFLATE: the host steps over each member's header in its copy. Without an engine, every
        // gzip file is in the second group.
        auto const gzip_pass = [&](bool engine) {
            for (size_t file = 0; file < files; ++file) {
                auto const& probe = probes[ids[file]];
                if (probe.plain || lane.possible_members[file]) continue;
                bool const fits =
                    hardware_limit_ != 0 && (probe.bgzf || probe.isize <= hardware_limit_);
                if (fits != engine) continue;
                auto const* const bytes = reinterpret_cast<unsigned char const*>(
                    lane.pinned->data() + table.compressed_offsets[file]
                );
                auto* out = lane.slots->data() + slot_at + 1;
                if (probe.bgzf) {
                    // The blocks must still be the ones probed, or they would overrun the slot.
                    auto const blocks = walk_bgzf(bytes, probe.compressed, [](auto...) {});
                    if (!blocks || blocks->first != probe.members ||
                        blocks->second != probe.isize) {
                        lane.fallback.push_back(ids[file]);
                        continue;
                    }
                    accept(file);
                    (void)walk_bgzf(
                        bytes, probe.compressed, [&](uint64_t at, uint64_t length, uint32_t size) {
                            add_member(bytes + at, length, size, out);
                            out += size;
                        }
                    );
                    continue;
                }
                auto const deflate_at = device_gzip::deflate_offset(
                    reinterpret_cast<char const*>(bytes), probe.compressed
                );
                if (deflate_at == 0) {
                    lane.fallback.push_back(ids[file]);
                    continue;
                }
                accept(file);
                add_member(bytes + deflate_at, probe.compressed - deflate_at - 8, probe.isize, out);
            }
        };
        gzip_pass(true);
        auto const engine_members = members;
        gzip_pass(false);
        lane.member_first.push_back(members);
        table.slots[lane.files.size()] = slot_at;
        lane.slot_bytes = slot_at;
        if (lane.files.empty()) return files;
 
        auto const s = cuda::stream_ref{*lane.stream};
        auto* const stream = s.get();
        auto const device_table = device_tables(lane);
        // The slots are rewritten below, so the sketch kernel reading them must be done.
        CUDDL_CUDA_TRY(s.wait(*lane.released));
        // Plain bytes land in their slots; slot_init_kernel then writes the guard bytes between.
        if (plain_at != 0) {
            CUDDL_CUDA_TRY(
                cuda::copy_bytes(
                    s,
                    cuda::std::span{lane.pinned->data(), plain_at},
                    device_span<char>{lane.slots->data(), plain_at}
                )
            );
        }
        if (compressed_at != plain_at) {
            CUDDL_CUDA_TRY(
                cuda::copy_bytes(
                    s,
                    cuda::std::span{lane.pinned->data() + plain_at, compressed_at - plain_at},
                    device_span<char>{lane.payload->data() + plain_at, compressed_at - plain_at}
                )
            );
        }
        CUDDL_CUDA_TRY(
            cuda::copy_bytes(
                s,
                cuda::std::span{lane.host->data(), inputs_bytes()},
                device_span<char>{lane.tables->data(), inputs_bytes()}
            )
        );
        CUDDL_CUDA_TRY(lane.uploaded->record(s));
        auto const count = static_cast<uint32_t>(lane.files.size());
        device_gzip::slot_init_kernel<<<(count + 255) / 256, 256, 0, stream>>>(
            lane.slots->data(),
            device_table.slots,
            count,
            device_table.first_header,
            device_table.header_count
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        if (engine_members != 0) {
            // The driver takes the parameters from host memory. nvCOMP's own engine path builds
            // them on the device and synchronises the stream to read them back, which stalls this
            // thread and every lane it drives. The engine reports only inflated sizes, so the
            // statuses read as success; finish still checks each size and CRC32.
            static_assert(nvcompSuccess == 0);
            static_assert(std::endian::native == std::endian::little);
            CUDDL_CUDA_TRY(cudaMemsetAsync(
                device_table.member_inflated, 0, engine_members * sizeof(uint64_t), stream
            ));
            CUDDL_CUDA_TRY(cudaMemsetAsync(
                device_table.member_statuses, 0, engine_members * sizeof(nvcompStatus_t), stream
            ));
            size_t error_index = 0;
            auto const submitted = batch_decompress_(
                engine_params_.data(), engine_params_.size(), 0, &error_index, stream
            );
            if (submitted != CUDA_SUCCESS) {
                return Err(
                    Error::resource(
                        "hardware decompression launch failed: " + std::to_string(submitted)
                    )
                );
            }
        }
        if (engine_members != members) {
            auto options = nvcompBatchedDeflateDecompressDefaultOpts;
            options.backend = NVCOMP_DECOMPRESS_BACKEND_CUDA;
            auto const first = engine_members;
            auto const inflated = nvcompBatchedDeflateDecompressAsync(
                reinterpret_cast<void const* const*>(device_table.member_sources + first),
                device_table.member_source_bytes + first,
                device_table.member_output_bytes + first,
                device_table.member_inflated + first,
                members - first,
                lane.temp->data(),
                nvcomp_temp_bytes_,
                reinterpret_cast<void* const*>(device_table.member_outputs + first),
                options,
                device_table.member_statuses + first,
                stream
            );
            if (inflated != nvcompSuccess) {
                return Err(
                    Error::resource("nvCOMP DEFLATE launch failed: " + std::to_string(inflated))
                );
            }
        }
        if (auto const gzip = count - static_cast<uint32_t>(plain); gzip != 0) {
            device_gzip::slot_crc_kernel<<<gzip, device_gzip::block_size, 0, stream>>>(
                lane.slots->data(), device_table.slots + plain, device_table.crcs + plain
            );
            CUDDL_CUDA_TRY(cudaGetLastError());
        }
        CUDDL_TRY(normalise(lane, slot_at));
        return files;
    }
 
   private:
    [[nodiscard]] Result<void> normalise(lane_state& lane, uint64_t slot_at) {
        auto const s = cuda::stream_ref{*lane.stream};
        auto const stream = s.get();
        auto const device_table = device_tables(lane);
        auto const count = static_cast<uint32_t>(lane.files.size());
        auto const grid = static_cast<uint32_t>(
            s.device().attribute(cuda::device_attributes::multiprocessor_count) * 8
        );
        device_gzip::header_kernel<<<grid, 256, 0, stream>>>(
            lane.slots->data(),
            slot_at,
            lane.headers->data(),
            header_capacity_,
            device_table.header_count
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        device_gzip::blank_header_kernel<<<grid, 256, 0, stream>>>(
            lane.slots->data(),
            lane.headers->data(),
            device_table.header_count,
            header_capacity_,
            device_table.slots,
            count,
            device_table.first_header
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        device_gzip::slot_finish_kernel<<<count, device_gzip::block_size, 0, stream>>>(
            lane.slots->data(),
            device_table.slots,
            device_table.first_header,
            device_table.first_char
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        CUDDL_CUDA_TRY(
            cudaMemsetAsync(device_table.kept, 0, count * sizeof(unsigned long long), stream)
        );
        device_gzip::kept_count_kernel<<<grid, device_gzip::block_size, 0, stream>>>(
            lane.slots->data(), device_table.slots, count, device_table.kept
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        auto temp_bytes = select_temp_bytes_;
        if (lane.decoded) {
            CUDDL_CUDA_TRY(
                cub::DeviceSelect::If(
                    lane.temp->data(),
                    temp_bytes,
                    lane.slots->data(),
                    device_table.selected,
                    static_cast<int64_t>(slot_at),
                    device_gzip::kept_byte{},
                    stream
                )
            );
        } else {
            // Inflation has consumed the compressed bytes; reuse their storage for compacted FASTA.
            CUDDL_CUDA_TRY(
                cub::DeviceSelect::If(
                    lane.temp->data(),
                    temp_bytes,
                    lane.slots->data(),
                    lane.payload->data(),
                    device_table.selected,
                    static_cast<int64_t>(slot_at),
                    device_gzip::kept_byte{},
                    stream
                )
            );
        }
        CUDDL_CUDA_TRY(
            cuda::copy_bytes(
                s,
                cuda::std::span{lane.tables->data() + inputs_bytes(), outputs_bytes()},
                device_span<char>{lane.host->data() + inputs_bytes(), outputs_bytes()}
            )
        );
        CUDDL_CUDA_TRY(lane.classified->record(s));
        return Ok();
    }
 
   public:
    [[nodiscard]] Result<void>
    submit_decoded(size_t lane_index, std::span<decoded_fasta const> sources) {
        if (sources.empty() || sources.size() > max_files || compressed_capacity_ != 0) {
            return Err(Error::invalid_argument("invalid decoded FASTA batch"));
        }
        auto& lane = lanes_[lane_index];
        CUDDL_TRY(cuda_try([&] { allocate(lane); }));
        lane.uploaded->sync();
        lane.decoded = true;
        lane.files.clear();
        lane.fallback.clear();
        auto const table = host_table(lane);
        uint64_t bytes = 0;
        for (auto const& source : sources) {
            if (source.bytes.size() >= slot_capacity_ - bytes) {
                return Err(Error::invalid_argument("decoded FASTA batch exceeds its lane"));
            }
            auto const file = lane.files.size();
            table.slots[file] = bytes;
            table.compressed_ptrs[file] = const_cast<char*>(source.bytes.data());
            table.compressed_bytes[file] = source.bytes.size();
            table.output_ptrs[file] = lane.slots->data() + bytes + 1;
            bytes += source.bytes.size() + 1;
            lane.files.push_back(source.id);
        }
        table.slots[sources.size()] = bytes;
        auto const s = cuda::stream_ref{*lane.stream};
        auto const device_table = device_tables(lane);
        CUDDL_CUDA_TRY(s.wait(*lane.released));
        CUDDL_CUDA_TRY(
            cuda::copy_bytes(
                s,
                cuda::std::span{lane.host->data(), inputs_bytes()},
                device_span<char>{lane.tables->data(), inputs_bytes()}
            )
        );
        auto const count = static_cast<uint32_t>(sources.size());
        device_gzip::slot_init_kernel<<<(count + 255) / 256, 256, 0, s.get()>>>(
            lane.slots->data(),
            device_table.slots,
            count,
            device_table.first_header,
            device_table.header_count
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        static_assert(CUDART_VERSION >= 13000, "cuDDL needs CUDA 13.0 or newer");
        cudaMemcpyAttributes attributes{};
        attributes.srcAccessOrder = cudaMemcpySrcAccessOrderAny;
        attributes.srcLocHint = {cudaMemLocationTypeHost, 0};
        attributes.flags = cudaMemcpyFlagPreferOverlapWithCompute;
        size_t attributes_index = 0;
        CUDDL_CUDA_TRY(cudaMemcpyBatchAsync(
            reinterpret_cast<void* const*>(table.output_ptrs),
            reinterpret_cast<void const* const*>(table.compressed_ptrs),
            table.compressed_bytes,
            count,
            &attributes,
            &attributes_index,
            1,
            s.get()
        ));
        CUDDL_CUDA_TRY(lane.uploaded->record(s));
        return normalise(lane, bytes);
    }
 
    [[nodiscard]] Result<void> wait_lanes(cuda::stream_ref consumer) {
        for (auto& lane : lanes_) {
            CUDDL_CUDA_TRY(lane.stream->sync());
        }
        CUDDL_CUDA_TRY(consumer.sync());
        return Ok();
    }
 
    [[nodiscard]] Result<void> finish(
        size_t lane_index,
        std::span<gzip_file_probe const> probes,
        std::vector<inflated_genome>& genomes,
        std::vector<size_t>& fallback
    ) {
        auto& lane = lanes_[lane_index];
        fallback.insert(fallback.end(), lane.fallback.begin(), lane.fallback.end());
        if (lane.files.empty()) return Ok();
        CUDDL_CUDA_TRY(lane.classified->sync());
        auto const table = host_table(lane);
        auto const overflow = *table.header_count > header_capacity_;
        auto const* bases = lane.decoded ? lane.slots->data() : lane.payload->data();
        uint64_t offset = 0;
        for (size_t file = 0; file < lane.files.size(); ++file) {
            auto const id = lane.files[file];
            auto const size = table.kept[file];
            bool intact = lane.decoded || probes[id].plain;
            if (!intact) {
                uint64_t inflated = 0;
                bool succeeded = true;
                for (auto m = lane.member_first[file]; m < lane.member_first[file + 1]; ++m) {
                    succeeded = succeeded && table.member_statuses[m] == nvcompSuccess;
                    inflated += table.member_inflated[m];
                }
                intact =
                    succeeded && inflated == probes[id].isize && table.crcs[file] == probes[id].crc;
            }
            auto const fasta =
                table.first_char[file] != '@' && table.first_header[file] < table.slots[file + 1];
            // Each header contributes one separator, including the first header.
            // ponytail: CPU-validate shorter outputs; count per-file headers if this dominates.
            auto const has_sequence = size > *table.header_count;
            if (!overflow && intact && fasta && has_sequence) {
                genomes.push_back({id, bases + offset, size});
            } else {
                fallback.push_back(id);
            }
            offset += size;
        }
        if (static_cast<int64_t>(offset) != *table.selected) {
            return Err(Error::resource("device FASTA compaction lost track of its bytes"));
        }
        return Ok();
    }
 
    [[nodiscard]] Result<void> release(size_t lane_index, cuda::stream_ref consumer) {
        CUDDL_CUDA_TRY(lanes_[lane_index].released->record(consumer));
        return Ok();
    }
 
   private:
    struct table_view {
        uint64_t* compressed_offsets;  // host only
        char** compressed_ptrs;        // decoded input only
        uint64_t* compressed_bytes;
        char** output_ptrs;  // decoded input only
        uint64_t* slots;     // max_files + 1
        // Gzip members, max_members each: raw DEFLATE in, inflated bytes out.
        char const** member_sources;
        uint64_t* member_source_bytes;
        uint64_t* member_output_bytes;
        char** member_outputs;
        // Outputs.
        uint64_t* member_inflated;
        nvcompStatus_t* member_statuses;
        uint32_t* crcs;
        unsigned long long* first_header;
        unsigned long long* kept;
        unsigned long long* header_count;
        int64_t* selected;
        char* first_char;
    };
 
    static constexpr size_t inputs_bytes() noexcept {
        return 5 * max_files * 8 + 8 + 4 * max_members * 8;
    }
 
    static constexpr size_t outputs_bytes() noexcept {
        return max_members * (8 + sizeof(nvcompStatus_t)) + max_files * (4 + 8 + 8 + 1) + 16;
    }
 
    static constexpr size_t host_table_bytes() noexcept {
        return inputs_bytes() + outputs_bytes();
    }
 
    static table_view view(char* base) {
        table_view view{};
        auto* at = base;
        auto take = [&]<typename T>(T*& field, size_t count) {
            field = reinterpret_cast<T*>(at);
            at += count * sizeof(T);
        };
        take(view.compressed_offsets, max_files);
        take(view.compressed_ptrs, max_files);
        take(view.compressed_bytes, max_files);
        take(view.output_ptrs, max_files);
        take(view.slots, max_files + 1);
        take(view.member_sources, max_members);
        take(view.member_source_bytes, max_members);
        take(view.member_output_bytes, max_members);
        take(view.member_outputs, max_members);
        take(view.member_inflated, max_members);
        take(view.first_header, max_files);
        take(view.kept, max_files);
        take(view.header_count, 1);
        take(view.selected, 1);
        take(view.member_statuses, max_members);
        take(view.crcs, max_files);
        take(view.first_char, max_files);
        return view;
    }
 
    struct lane_state {
        bool decoded = false;
        std::optional<cuda::stream> stream;
        std::optional<cuda::event> uploaded, classified, released;
        std::optional<cuda::buffer<char, cuda::mr::host_accessible, cuda::mr::device_accessible>>
            pinned, host;
        std::optional<cuda::device_buffer<char>> payload, slots, tables, temp;
        std::optional<cuda::device_buffer<uint64_t>> headers;
        std::vector<size_t> files, fallback;
        std::vector<uint32_t> member_first;  // each file's first member, then the member count
        std::array<bool, max_files> possible_members;
        uint64_t slot_bytes = 0;
    };
 
    static table_view host_table(lane_state& lane) {
        return view(lane.host->data());
    }
 
    static table_view device_tables(lane_state& lane) {
        return view(lane.tables->data());
    }
 
    void allocate(lane_state& lane) {
        if (lane.slots) return;
        auto const s = cuda::stream_ref{*lane.stream};
        auto const device = s.device();
        // Compressed inputs and inflated outputs: memory the decompression engine accepts on
        // devices that have one.
        auto const inflate_buffer = [&](size_t bytes) -> cuda::device_buffer<char> {
            if (hardware_limit_ == 0) {
                return cuda::make_device_buffer<char>(s, device, bytes, cuda::no_init);
            }
            return cuda::make_buffer<char>(
                s, device_gzip::hardware_inflate_pool(device), bytes, cuda::no_init
            );
        };
        if (compressed_capacity_ != 0) {
            lane.pinned.emplace(
                s,
                device_pinned_pool(device, pinned_placement::host),
                compressed_capacity_,
                cuda::no_init
            );
            lane.payload.emplace(inflate_buffer(std::max(compressed_capacity_, slot_capacity_)));
        }
        lane.host.emplace(
            s, device_pinned_pool(device, pinned_placement::host), host_table_bytes(), cuda::no_init
        );
        // nvCOMP's pointer, size and status arrays: the engine needs them in its memory too.
        lane.tables.emplace(inflate_buffer(host_table_bytes()));
        lane.headers.emplace(
            cuda::make_device_buffer<uint64_t>(s, device, header_capacity_, cuda::no_init)
        );
        lane.temp.emplace(
            cuda::make_device_buffer<char>(
                s, device, std::max(nvcomp_temp_bytes_, select_temp_bytes_) + 1, cuda::no_init
            )
        );
        // Sixteen bytes of tail padding let the header scan read whole words. Allocated last, so
        // a lane with slots has every buffer.
        lane.slots.emplace(inflate_buffer(slot_capacity_ + 16));
        if (streaming_window_ != 0) {
            cudaStreamAttrValue policy{};
            policy.accessPolicyWindow.base_ptr = lane.slots->data();
            policy.accessPolicyWindow.num_bytes = streaming_window_;
            policy.accessPolicyWindow.hitRatio = 0.0f;
            policy.accessPolicyWindow.hitProp = cudaAccessPropertyNormal;
            policy.accessPolicyWindow.missProp = cudaAccessPropertyStreaming;
            auto const status =
                cudaStreamSetAttribute(s.get(), cudaStreamAttributeAccessPolicyWindow, &policy);
            if (status != cudaSuccess) throw std::runtime_error(cudaGetErrorString(status));
        }
    }
 
    size_t workers_;
    // Lane reads wait on storage, so they keep a parallel file system's depth in flight.
    worker_pool readers_{
        ablated("read-depth") ? workers_ : std::max(workers_, parallel_read_depth)
    };
    size_t compressed_capacity_ = 0, slot_capacity_ = 0, streaming_window_ = 0;
    size_t hardware_limit_ = 0;  // device_gzip::hardware_inflate_limit of the stream's device
    PFN_cuMemBatchDecompressAsync_v12060 batch_decompress_ = nullptr;
    std::vector<CUmemDecompressParams> engine_params_;  // one batch's, reused
    size_t nvcomp_temp_bytes_ = 0, select_temp_bytes_ = 0;
    uint32_t header_capacity_ = 0;
    std::array<lane_state, lane_count> lanes_;
};
 
}  // namespace cuddl::detail