From 9e6d47cb883a741883c818f705973473ed9ce0dd Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Tue, 14 Jul 2026 05:04:27 +0000 Subject: [PATCH 1/8] NanoVDB CUDA: parallel bit-identical GridChecksum (slicing-by-4 + GF(2) combine) Port of Agent A's exp/14. Replaces the two serial anti-patterns in blockedCRC32 with parallel equivalents: - Per-block CRC via slicing-by-4: the 256-entry base LUT is staged into shared memory and three derived slice tables are built in place, so the dependent CRC chain advances four bytes per step with four shared-memory lookups instead of one byte per L2 lookup. (Not __constant__: divergent indices serialize there.) - The final fold (a single thread CRCing the whole block-CRC array - megabytes for multi-GB grids) becomes parallel per-chunk CRCs plus a GF(2) crc32_combine (zlib construction: one shift operator per fixed chunk length, applied per fold). Small inputs keep the single-thread path. Checksum VALUES are bit-identical (slicing-by-4 and the GF(2) combine are both value-preserving); all checksum_full goldens pass on dragon/emu/crawler/ wdas_cloud. The pre-v32.6.0 crc32TailOld path is untouched. A/B vs the Tranche-1 base (RTX PRO 6000, 3 rounds): dragon +84.7%, emu +92.0%, crawler +93.2%, wdas_cloud +94.1% (~6.5x / 12.5x / 15x / 17x); transient memory unchanged. Co-Authored-By: Claude Fable 5 Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/GridChecksum.cuh | 116 ++++++++++++++++++-- 1 file changed, 108 insertions(+), 8 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh b/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh index 2d9bc02af5..5c84326b2e 100644 --- a/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh @@ -109,6 +109,91 @@ inline unique_ptr createCrc32Lut(size_t extra = 0, cudaStream_t stream return lut; } +/// @brief Cuda kernel computing per-block CRC32 checksums via slicing-by-4. +/// The 256-entry base LUT is staged into shared memory and three derived +/// slice tables are built in place, so the four divergent lookups per 4-byte +/// step hit shared memory and the dependent update chain advances four bytes +/// per step instead of one. Bit-identical to the byte-serial crc32(). +/// The final block absorbs any remainder of @c totalSize. +__global__ inline void crc32SlicedKernel(const void *d_data, uint32_t* d_blockCRC, uint64_t blockCount, uint32_t log2BlockSize, uint64_t totalSize, const uint32_t *d_lut) +{ + __shared__ uint32_t sLut[4][256]; + for (uint32_t i = threadIdx.x; i < 256; i += blockDim.x) sLut[0][i] = d_lut[i]; + __syncthreads(); + for (int k = 1; k < 4; ++k) { + for (uint32_t i = threadIdx.x; i < 256; i += blockDim.x) { + const uint32_t c = sLut[k-1][i]; + sLut[k][i] = (c >> 8) ^ sLut[0][c & 0xffu]; + } + __syncthreads(); + } + const uint64_t tid = blockIdx.x * (uint64_t)blockDim.x + threadIdx.x; + if (tid >= blockCount) return; + const uint8_t *p = (const uint8_t*)d_data + (tid << log2BlockSize); + uint64_t n = uint64_t(1) << log2BlockSize; + if (tid + 1 == blockCount) n += totalSize - (blockCount << log2BlockSize); + uint32_t crc = ~0u; + // blocks start at power-of-two offsets of a 32B-aligned buffer -> 4B aligned + const uint32_t *w = (const uint32_t*)p; + for (uint64_t i = 0, nw = n >> 2; i < nw; ++i) { + const uint32_t x = crc ^ w[i]; + crc = sLut[3][x & 0xffu] ^ sLut[2][(x >> 8) & 0xffu] ^ sLut[1][(x >> 16) & 0xffu] ^ sLut[0][x >> 24]; + } + for (uint64_t i = n & ~3ull; i < n; ++i) { + crc ^= p[i]; + for (int j = 0; j < 8; ++j) crc = (crc >> 1) ^ (0xEDB88320u & (-(crc & 1u))); + } + d_blockCRC[tid] = ~crc; +} + +/// @brief y = M x over GF(2), M given as 32 column words +__device__ inline uint32_t crc32Gf2MatTimes(const uint32_t *mat, uint32_t vec) +{ + uint32_t sum = 0; + while (vec) { + if (vec & 1u) sum ^= *mat; + vec >>= 1; + ++mat; + } + return sum; +} + +/// @brief Single-thread kernel folding per-chunk CRCs into the CRC of the +/// whole stream, using crc(A||B) = shift(crc(A), len(B)) ^ crc(B) where +/// shift is the GF(2) operator advancing a CRC past len(B) zero bytes +/// (zlib crc32_combine). All chunks share one length so a single operator +/// (built once by binary exponentiation) serves every fold; the final, +/// possibly shorter, chunk gets its own. O(chunkCount x 32) - negligible. +/// Bit-identical to a serial crc32 over the concatenated stream. +__global__ inline void crc32CombineKernel(const uint32_t *d_chunkCRC, uint64_t chunkCount, uint64_t chunkBytes, uint64_t lastChunkBytes, uint32_t *d_crc) +{ + uint32_t oddBuf[32], evenBuf[32], acc[32], accLast[32]; + auto buildOp = [&](uint32_t *dst, uint64_t bits) { + uint32_t *cur = oddBuf, *nxt = evenBuf; + cur[0] = 0xEDB88320u;// operator for a single zero bit + for (int n = 1; n < 32; ++n) cur[n] = 1u << (n - 1); + for (int n = 0; n < 32; ++n) dst[n] = 1u << n;// identity + while (bits) { + if (bits & 1ull) { + for (int n = 0; n < 32; ++n) dst[n] = crc32Gf2MatTimes(cur, dst[n]); + } + bits >>= 1; + if (bits) { + for (int n = 0; n < 32; ++n) nxt[n] = crc32Gf2MatTimes(cur, cur[n]); + uint32_t *t = cur; cur = nxt; nxt = t; + } + } + }; + buildOp(acc, chunkBytes * 8ull); + if (lastChunkBytes != chunkBytes) buildOp(accLast, lastChunkBytes * 8ull); + uint32_t crc = d_chunkCRC[0]; + for (uint64_t i = 1; i < chunkCount; ++i) { + const uint32_t *op = (i + 1 == chunkCount && lastChunkBytes != chunkBytes) ? accLast : acc; + crc = crc32Gf2MatTimes(op, crc) ^ d_chunkCRC[i]; + } + *d_crc = crc; +} + /// @brief Compute CRC32 checksum of 4K block /// @param d_data device pointer to start of data /// @param size number of bytes @@ -121,14 +206,29 @@ inline void blockedCRC32(const void *d_data, size_t size, const uint32_t *d_lut, const uint64_t checksumCount = size >> NANOVDB_CRC32_LOG2_BLOCK_SIZE;// 4 KB (4096 byte) unique_ptr buffer(checksumCount, stream);// for checksums of 4 KB blocks uint32_t *d_checksums = buffer.get(); - lambdaKernel<<>>(checksumCount, [=] __device__(size_t tid) { - uint32_t blockSize = 1 << NANOVDB_CRC32_LOG2_BLOCK_SIZE; - if (tid+1 == checksumCount) blockSize += size - (checksumCount<>>(1, [=] __device__(size_t) {// Compute CRC32 of all the 4K blocks - *d_crc = crc32((const uint8_t*)d_checksums, checksumCount*sizeof(uint32_t), d_lut); - }); cudaCheckError(); + crc32SlicedKernel<<>>( + d_data, d_checksums, checksumCount, NANOVDB_CRC32_LOG2_BLOCK_SIZE, size, d_lut); + cudaCheckError(); + // CRC of the block-checksum array itself. The former single-thread pass + // over checksumCount*4 bytes (megabytes for multi-GB grids) is replaced + // by parallel per-chunk CRCs plus a GF(2) combine - bit-identical result. + const uint64_t checksumBytes = checksumCount*sizeof(uint32_t); + constexpr uint64_t log2ChunkSize = 12, chunkSize = uint64_t(1) << log2ChunkSize; + if (checksumBytes <= 2*chunkSize) {// small: single-thread CRC is fine + lambdaKernel<<<1, 1, 0, stream>>>(1, [=] __device__(size_t) { + *d_crc = crc32((const uint8_t*)d_checksums, checksumBytes, d_lut); + }); cudaCheckError(); + } else { + const uint64_t chunkCount = checksumBytes >> log2ChunkSize;// final chunk absorbs the remainder + const uint64_t lastChunkBytes = chunkSize + (checksumBytes - (chunkCount << log2ChunkSize)); + unique_ptr chunkBuffer(chunkCount, stream); + uint32_t *d_chunkCRC = chunkBuffer.get(); + crc32SlicedKernel<<>>( + d_checksums, d_chunkCRC, chunkCount, log2ChunkSize, checksumBytes, d_lut); + cudaCheckError(); + crc32CombineKernel<<<1, 1, 0, stream>>>(d_chunkCRC, chunkCount, chunkSize, lastChunkBytes, d_crc); + cudaCheckError(); + } }// void cudaBlockedCRC32(const void *d_data, size_t size, const uint32_t *d_lut, uint32_t *d_crc, cudaStream_t stream) /// @brief Compute CRC32 checksum of 4K block From 44d5e657a78f0d2b72bba99e07f9174b87a9c36f Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Tue, 14 Jul 2026 05:08:39 +0000 Subject: [PATCH 2/8] NanoVDB CUDA: cooperative GridStats reductions (warp-per-leaf, block-per-node) Port of Agent A's exp/13. Replaces the thread-per-node serial reductions: - processLeaf: one warp per leaf - lanes stride the 512 voxel slots (mask-gated) and merge partial statistics through shared memory; lane 0 keeps the bbox update and the final store. - processInternal: one 128-thread block per internal node striding the 4,096/32,768-entry child table with a shared-memory tree reduction of Stats + bbox. The former kernel gave each thread up to 32,768 serial child visits and ran an entire lower level of a 181M-voxel grid on ~9 warps of a 188-SM GPU. min/max are order-independent, so stats_minmax output is bit-identical - all stats_minmax goldens pass on dragon/emu/crawler/wdas_cloud. avg/stddev (stats_all) use a Welford-correct parallel merge but are NOT bit-identical to the serial result (floating-point reduction order differs); stats_all stays a bench-only measurement, never a byte-exact claim. A/B vs the Tranche-1 base (RTX PRO 6000, 3 rounds): stats_minmax dragon +87.3% emu +82.5% crawler +80.3% wdas +83.5% (5-8x) stats_all dragon +71.9% emu +51.7% crawler +29.8% wdas +35.1% transient memory unchanged. Co-Authored-By: Claude Fable 5 Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/GridStats.cuh | 129 ++++++++++++++++------- 1 file changed, 91 insertions(+), 38 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/GridStats.cuh b/nanovdb/nanovdb/tools/cuda/GridStats.cuh index f8dd01c025..3504fe41d0 100644 --- a/nanovdb/nanovdb/tools/cuda/GridStats.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridStats.cuh @@ -58,65 +58,117 @@ public: namespace {// define cuda kernels in an unnamed namespace +// One warp per leaf: lanes stride the 512 voxel slots (mask-gated) and merge +// their partial statistics through shared memory; lane 0 handles the bbox +// update and the final store. Launch with 128 threads (4 warps) per block. template __global__ void processLeaf(NodeManager *d_nodeMgr, StatsT *d_stats) { - const uint32_t tid = blockIdx.x * blockDim.x + threadIdx.x; + constexpr uint32_t WarpsPerBlock = 4; + __shared__ StatsT sStats[WarpsPerBlock * 32]; + + const uint32_t warpID = threadIdx.x >> 5, lane = threadIdx.x & 31u; + const uint32_t tid = blockIdx.x * WarpsPerBlock + warpID;// leaf index if (tid >= d_nodeMgr->leafCount()) return; auto &d_leaf = d_nodeMgr->leaf(tid); - if (d_leaf.updateBBox()) {// updates active bounding box (also updates data->mFlags) and return true if non-empty + bool nonEmpty = false; + if (lane == 0) nonEmpty = d_leaf.updateBBox();// updates active bounding box (also updates data->mFlags) + nonEmpty = __shfl_sync(0xffffffffu, nonEmpty, 0); + + if (nonEmpty) { if constexpr(StatsT::hasStats()) { StatsT stats; - for (auto it = d_leaf.cbeginValueOn(); it; ++it) stats.add(*it); - if constexpr(StatsT::hasAverage()) { - d_stats[tid] = stats; - *reinterpret_cast(&d_leaf.mMinimum) = tid; - } else { - stats.setStats(d_leaf); + const auto &mask = d_leaf.valueMask(); + for (uint32_t i = lane; i < 512; i += 32) + if (mask.isOn(i)) stats.add(d_leaf.getValue(i)); + StatsT *sWarp = sStats + (warpID << 5); + sWarp[lane] = stats; + __syncwarp(); + for (uint32_t d = 16; d; d >>= 1) { + if (lane < d) sWarp[lane].add(sWarp[lane + d]); + __syncwarp(); + } + if (lane == 0) { + if constexpr(StatsT::hasAverage()) { + d_stats[tid] = sWarp[0]; + *reinterpret_cast(&d_leaf.mMinimum) = tid; + } else { + sWarp[0].setStats(d_leaf); + } } } } - d_leaf.mFlags &= ~uint8_t(1u);// enable rendering + if (lane == 0) d_leaf.mFlags &= ~uint8_t(1u);// enable rendering }// processLeaf +// One block per internal node: threads stride the node's child table (4096 or +// 32768 entries) and merge partial bboxes/statistics through shared memory - +// the former one-thread-per-node kernel serialized up to 32768 child visits +// per thread and left the device nearly empty at typical node counts. template __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) { using ChildT = typename NanoNode::type; - uint32_t nodeID = blockIdx.x * blockDim.x + threadIdx.x;// thread id (reused below to avoid compiler warning) + using NodeT = typename NanoNode::type; + constexpr uint32_t Threads = 128; + __shared__ StatsT sStats[Threads]; + __shared__ CoordBBox sBBox[Threads]; + __shared__ int sSlot;// any one child's d_stats slot, claimed for this node + + const uint32_t nodeID = blockIdx.x; + const uint32_t tID = threadIdx.x; if (nodeID >= d_nodeMgr->nodeCount(LEVEL)) return; auto &d_node = d_nodeMgr->template node(nodeID); - auto &bbox = d_node.mBBox; - bbox = CoordBBox();// empty bbox - StatsT stats; + if (tID == 0) sSlot = -1; + __syncthreads(); - for (auto it = d_node.beginChild(); it; ++it) { - auto &child = *it; - bbox.expand( child.bbox() ); - if constexpr(StatsT::hasAverage()) { - nodeID = *reinterpret_cast(&child.mMinimum); - StatsT &s = d_stats[nodeID]; - s.setStats(child); - stats.add(s); - } else if constexpr(StatsT::hasMinMax()) { - stats.add(child.minimum()); - stats.add(child.maximum()); + CoordBBox bbox;// empty + StatsT stats; + int mySlot = -1; + for (uint32_t i = tID; i < NodeT::SIZE; i += Threads) { + if (d_node.childMask().isOn(i)) { + auto &child = *d_node.getChild(i); + bbox.expand( child.bbox() ); + if constexpr(StatsT::hasAverage()) { + const int slot = *reinterpret_cast(&child.mMinimum); + StatsT &s = d_stats[slot]; + s.setStats(child); + stats.add(s); + mySlot = slot; + } else if constexpr(StatsT::hasMinMax()) { + stats.add(child.minimum()); + stats.add(child.maximum()); + } + } else if (d_node.valueMask().isOn(i)) { + const Coord ijk = d_node.offsetToGlobalCoord(i); + bbox[0].minComponent(ijk); + bbox[1].maxComponent(ijk + Coord(ChildT::DIM - 1)); + if constexpr(StatsT::hasStats()) stats.add(d_node.data()->getValue(i), ChildT::NUM_VALUES); } } - for (auto it = d_node.cbeginValueOn(); it; ++it) { - const Coord ijk = it.getCoord(); - bbox[0].minComponent(ijk); - bbox[1].maxComponent(ijk + Coord(ChildT::DIM - 1)); - if constexpr(StatsT::hasStats()) stats.add(*it, ChildT::NUM_VALUES); + if constexpr(StatsT::hasAverage()) + if (mySlot >= 0) atomicMax(&sSlot, mySlot); + sStats[tID] = stats; + sBBox[tID] = bbox; + __syncthreads(); + for (uint32_t d = Threads >> 1; d; d >>= 1) { + if (tID < d) { + sStats[tID].add(sStats[tID + d]); + sBBox[tID].expand(sBBox[tID + d]); + } + __syncthreads(); } - if constexpr(StatsT::hasAverage()) { - d_stats[nodeID] = stats; - *reinterpret_cast(&d_node.mMinimum) = nodeID; - } else if constexpr(StatsT::hasMinMax()) { - stats.setStats(d_node); + if (tID == 0) { + d_node.mBBox = sBBox[0]; + if constexpr(StatsT::hasAverage()) { + d_stats[sSlot] = sStats[0]; + *reinterpret_cast(&d_node.mMinimum) = sSlot; + } else if constexpr(StatsT::hasMinMax()) { + sStats[0].setStats(d_node); + } + d_node.mFlags &= ~uint64_t(1u);// enable rendering } - d_node.mFlags &= ~uint64_t(1u);// enable rendering }// processInternal template @@ -201,11 +253,12 @@ void GridStats::update(NanoGrid *d_grid, cudaStream_t st if constexpr(StatsT::hasAverage()) cudaCheck(util::cuda::mallocAsync((void**)&d_stats, nodeCount[0]*sizeof(StatsT), stream)); - processLeaf<<>>(d_nodeMgr, d_stats); + // warp per leaf (4 warps per 128-thread block); block per internal node + processLeaf<<>>(d_nodeMgr, d_stats); - processInternal<<>>(d_nodeMgr, d_stats); + if (nodeCount[1]) processInternal<<>>(d_nodeMgr, d_stats); - processInternal<<>>(d_nodeMgr, d_stats); + if (nodeCount[2]) processInternal<<>>(d_nodeMgr, d_stats); processRootAndGrid<<<1, 1, 0, stream>>>(d_nodeMgr, d_stats); From 5fcaa5a347649158800dd111c1b856a860a390f8 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Tue, 14 Jul 2026 05:13:03 +0000 Subject: [PATCH 3/8] NanoVDB CUDA: coalesced IndexToGrid node/leaf remap Port of Agent A's exp/20. Consecutive threads now process consecutive table entries and leaf values - the former mapping gave each thread a consecutive RUN (a 32-entry stride across the warp for internal nodes, 8 for leaves) on every iteration. The 4 KB upper-node value/child masks are copied cooperatively across the block instead of member-wise by thread 0 while the rest idled. Combines with the Tranche-1 scoped zero-init already on this file: the leaf kernel's stats/padding-gap zeroing is preserved inside the thread-0 block, and the coalesced value loop still writes every mValues[i], so output stays bit-identical - all indextogrid goldens pass on dragon/emu/crawler/wdas_cloud. A/B vs the Tranche-1 base (RTX PRO 6000, 3 rounds): dragon +23%, emu +21%, crawler +18%, wdas_cloud +14%; transient memory unchanged. The remap gain is on top of Tranche 1's scoped init, so absolute indextogrid is now faster than both master and A's exp/20. Co-Authored-By: Claude Fable 5 Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh | 34 +++++++++++++--------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh b/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh index 6d8bde11e6..63d5151339 100644 --- a/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh +++ b/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh @@ -218,11 +218,11 @@ __global__ void processNodesKernel(typename IndexToGrid::NodeAccessor auto &srcNode = nodeAcc->template srcNode(blockIdx.x); auto &dstNode = nodeAcc->template dstNode(blockIdx.x); - if (threadIdx.x == 0 && threadIdx.y == 0) { + const int tid = threadIdx.y*blockDim.x + threadIdx.x; + const int nThreads = blockDim.x*blockDim.y; + if (tid == 0) { dstNode.mBBox = srcNode.mBBox; dstNode.mFlags = srcNode.mFlags; - dstNode.mValueMask = srcNode.mValueMask; - dstNode.mChildMask = srcNode.mChildMask; auto &srcGrid = nodeAcc->srcGrid(); if (srcGrid.hasMinMax()) { dstNode.mMinimum = srcValues[srcNode.mMinimum]; @@ -233,9 +233,16 @@ __global__ void processNodesKernel(typename IndexToGrid::NodeAccessor if (srcGrid.hasStdDeviation()) dstNode.mStdDevi = srcValues[srcNode.mStdDevi]; } } - const int off = blockDim.x*blockDim.y*threadIdx.x + blockDim.x*threadIdx.y; - for (int threadIdx_z=0; threadIdx_z::NodeAccessor static_assert(!BuildTraits::is_special, "Invalid destination type!"); auto &srcLeaf = nodeAcc->template srcNode<0>(blockIdx.x); auto &dstLeaf = nodeAcc->template dstNode(blockIdx.x); - if (threadIdx.x == 0 && threadIdx.y == 0) { + const int tid = threadIdx.y*blockDim.x + threadIdx.x; + const int nThreads = blockDim.x*blockDim.y; + if (tid == 0) { dstLeaf.mBBoxMin = srcLeaf.mBBoxMin; for (int i=0; i<3; ++i) dstLeaf.mBBoxDif[i] = srcLeaf.mBBoxDif[i]; dstLeaf.mFlags = srcLeaf.mFlags; @@ -284,12 +293,11 @@ __global__ void processLeafsKernel(typename IndexToGrid::NodeAccessor if (srcGrid.hasStdDeviation()) dstLeaf.mStdDevi = srcValues[srcLeaf.getDev()]; } } - const int off = blockDim.x*blockDim.y*threadIdx.x + blockDim.x*threadIdx.y; - auto *dst = dstLeaf.mValues + off; - for (int threadIdx_z=0; threadIdx_z Date: Tue, 14 Jul 2026 22:47:47 +0000 Subject: [PATCH 4/8] NanoVDB CUDA: GridStats review fixes - guard leaf launch + unique per-node stats slot Addresses PR review on the cooperative GridStats: - processLeaf is now launched only when nodeCount[0] > 0, matching the guards on the internal-node launches (a leaf-less grid would otherwise be a <<<0,...>>> launch). - The average path now gives each node its own d_stats slot instead of reusing one of its children's. A fully-tiled internal node (active value tiles, no children) never claimed a child slot, so the old code wrote d_stats[-1] (out of bounds) and propagated an invalid slot to its parent. d_stats is now sized for all nodes (leaves, then lower, then upper) and each node writes slot leafCount + levelOffset + nodeID; the sSlot/atomicMax machinery is removed. stats_minmax is unaffected (it never touches d_stats) and stays byte-exact; stats_all output is unchanged for grids without childless internal nodes and correct for those that have them. compute-sanitizer memcheck clean. Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/GridStats.cuh | 27 ++++++++++++------------ 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/GridStats.cuh b/nanovdb/nanovdb/tools/cuda/GridStats.cuh index 3504fe41d0..4ba6658f8c 100644 --- a/nanovdb/nanovdb/tools/cuda/GridStats.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridStats.cuh @@ -114,28 +114,22 @@ __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) constexpr uint32_t Threads = 128; __shared__ StatsT sStats[Threads]; __shared__ CoordBBox sBBox[Threads]; - __shared__ int sSlot;// any one child's d_stats slot, claimed for this node const uint32_t nodeID = blockIdx.x; const uint32_t tID = threadIdx.x; if (nodeID >= d_nodeMgr->nodeCount(LEVEL)) return; auto &d_node = d_nodeMgr->template node(nodeID); - if (tID == 0) sSlot = -1; - __syncthreads(); CoordBBox bbox;// empty StatsT stats; - int mySlot = -1; for (uint32_t i = tID; i < NodeT::SIZE; i += Threads) { if (d_node.childMask().isOn(i)) { auto &child = *d_node.getChild(i); bbox.expand( child.bbox() ); if constexpr(StatsT::hasAverage()) { - const int slot = *reinterpret_cast(&child.mMinimum); - StatsT &s = d_stats[slot]; + StatsT &s = d_stats[*reinterpret_cast(&child.mMinimum)]; s.setStats(child); stats.add(s); - mySlot = slot; } else if constexpr(StatsT::hasMinMax()) { stats.add(child.minimum()); stats.add(child.maximum()); @@ -147,8 +141,6 @@ __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) if constexpr(StatsT::hasStats()) stats.add(d_node.data()->getValue(i), ChildT::NUM_VALUES); } } - if constexpr(StatsT::hasAverage()) - if (mySlot >= 0) atomicMax(&sSlot, mySlot); sStats[tID] = stats; sBBox[tID] = bbox; __syncthreads(); @@ -162,8 +154,15 @@ __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) if (tID == 0) { d_node.mBBox = sBBox[0]; if constexpr(StatsT::hasAverage()) { - d_stats[sSlot] = sStats[0]; - *reinterpret_cast(&d_node.mMinimum) = sSlot; + // Each node writes its OWN unique d_stats slot (leaves occupy + // [0, leafCount), then lower nodes, then upper nodes), so a node + // with no children still has a valid slot. The previous scheme + // reused a child's slot, writing d_stats[-1] for a childless + // (fully-tiled) internal node - an out-of-bounds write. + const uint32_t slot = d_nodeMgr->leafCount() + + (LEVEL == 2 ? d_nodeMgr->nodeCount(1) : 0u) + nodeID; + d_stats[slot] = sStats[0]; + *reinterpret_cast(&d_node.mMinimum) = slot; } else if constexpr(StatsT::hasMinMax()) { sStats[0].setStats(d_node); } @@ -251,10 +250,12 @@ void GridStats::update(NanoGrid *d_grid, cudaStream_t st StatsT *d_stats = nullptr; - if constexpr(StatsT::hasAverage()) cudaCheck(util::cuda::mallocAsync((void**)&d_stats, nodeCount[0]*sizeof(StatsT), stream)); + // One d_stats slot per node (leaves, then lower, then upper) so every node + // has its own slot - see processInternal. + if constexpr(StatsT::hasAverage()) cudaCheck(util::cuda::mallocAsync((void**)&d_stats, (nodeCount[0]+nodeCount[1]+nodeCount[2])*sizeof(StatsT), stream)); // warp per leaf (4 warps per 128-thread block); block per internal node - processLeaf<<>>(d_nodeMgr, d_stats); + if (nodeCount[0]) processLeaf<<>>(d_nodeMgr, d_stats); if (nodeCount[1]) processInternal<<>>(d_nodeMgr, d_stats); From 2e30d612830c6f69244dd3acaca62d492ab2fbf9 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Sun, 19 Jul 2026 01:53:03 +0000 Subject: [PATCH 5/8] NanoVDB CUDA: precompute GridChecksum combine operators on the host crc32CombineKernel rebuilt its two GF(2) shift operators by binary exponentiation on a single device thread every call - a fixed ~480 us that dominated the checksum for small and medium grids (profiling showed it as large as the entire tail-streaming pass). The operators depend only on the chunk length, so they are now built on the host in microseconds and uploaded, and the device kernel is reduced to the O(chunkCount) fold. Bit-identical: all checksum goldens pass on dragon/emu/crawler/wdas_cloud. A/B vs the parallel-checksum base (RTX PRO 6000, 3 rounds): dragon +40.7% (1.69x) emu +15.3% crawler +9.8% wdas_cloud +1.9% The gain shrinks with grid size because it removes a fixed cost; transient memory unchanged. Signed-off-by: Jonathan Swartz Co-Authored-By: Claude Opus 4.8 (1M context) --- nanovdb/nanovdb/tools/cuda/GridChecksum.cuh | 68 ++++++++++++--------- pendingchanges/nanovdbchecksumcombine.txt | 5 ++ 2 files changed, 43 insertions(+), 30 deletions(-) create mode 100644 pendingchanges/nanovdbchecksumcombine.txt diff --git a/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh b/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh index 5c84326b2e..f5585e8b96 100644 --- a/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridChecksum.cuh @@ -147,7 +147,7 @@ __global__ inline void crc32SlicedKernel(const void *d_data, uint32_t* d_blockCR } /// @brief y = M x over GF(2), M given as 32 column words -__device__ inline uint32_t crc32Gf2MatTimes(const uint32_t *mat, uint32_t vec) +__host__ __device__ inline uint32_t crc32Gf2MatTimes(const uint32_t *mat, uint32_t vec) { uint32_t sum = 0; while (vec) { @@ -158,37 +158,35 @@ __device__ inline uint32_t crc32Gf2MatTimes(const uint32_t *mat, uint32_t vec) return sum; } -/// @brief Single-thread kernel folding per-chunk CRCs into the CRC of the -/// whole stream, using crc(A||B) = shift(crc(A), len(B)) ^ crc(B) where -/// shift is the GF(2) operator advancing a CRC past len(B) zero bytes -/// (zlib crc32_combine). All chunks share one length so a single operator -/// (built once by binary exponentiation) serves every fold; the final, -/// possibly shorter, chunk gets its own. O(chunkCount x 32) - negligible. -/// Bit-identical to a serial crc32 over the concatenated stream. -__global__ inline void crc32CombineKernel(const uint32_t *d_chunkCRC, uint64_t chunkCount, uint64_t chunkBytes, uint64_t lastChunkBytes, uint32_t *d_crc) +/// @brief Build into @c dst[32] the GF(2) operator that advances a CRC past +/// @c bits zero bits, by binary exponentiation of the single-bit operator +/// (zlib crc32_combine construction). O(log bits x 32) - a few microseconds on +/// the host, which is where it is called so the device combine stays a cheap fold. +__host__ __device__ inline void crc32BuildShiftOp(uint32_t *dst, uint64_t bits) +{ + uint32_t oddBuf[32], evenBuf[32]; + uint32_t *cur = oddBuf, *nxt = evenBuf; + cur[0] = 0xEDB88320u;// operator for a single zero bit + for (int n = 1; n < 32; ++n) cur[n] = 1u << (n - 1); + for (int n = 0; n < 32; ++n) dst[n] = 1u << n;// identity + while (bits) { + if (bits & 1ull) for (int n = 0; n < 32; ++n) dst[n] = crc32Gf2MatTimes(cur, dst[n]); + bits >>= 1; + if (bits) { for (int n = 0; n < 32; ++n) nxt[n] = crc32Gf2MatTimes(cur, cur[n]); uint32_t *t = cur; cur = nxt; nxt = t; } + } +} + +/// @brief Single-thread kernel folding per-chunk CRCs into the CRC of the whole +/// stream, using crc(A||B) = shift(crc(A), len(B)) ^ crc(B). The two GF(2) shift +/// operators - @c d_acc for a full chunk and @c d_accLast for the final, +/// possibly shorter, chunk - are precomputed on the host, so this is just an +/// O(chunkCount) fold. Bit-identical to a serial crc32 over the concatenated +/// stream. (@c d_accLast equals @c d_acc when the last chunk is full.) +__global__ inline void crc32CombineKernel(const uint32_t *d_chunkCRC, uint64_t chunkCount, const uint32_t *d_acc, const uint32_t *d_accLast, uint32_t *d_crc) { - uint32_t oddBuf[32], evenBuf[32], acc[32], accLast[32]; - auto buildOp = [&](uint32_t *dst, uint64_t bits) { - uint32_t *cur = oddBuf, *nxt = evenBuf; - cur[0] = 0xEDB88320u;// operator for a single zero bit - for (int n = 1; n < 32; ++n) cur[n] = 1u << (n - 1); - for (int n = 0; n < 32; ++n) dst[n] = 1u << n;// identity - while (bits) { - if (bits & 1ull) { - for (int n = 0; n < 32; ++n) dst[n] = crc32Gf2MatTimes(cur, dst[n]); - } - bits >>= 1; - if (bits) { - for (int n = 0; n < 32; ++n) nxt[n] = crc32Gf2MatTimes(cur, cur[n]); - uint32_t *t = cur; cur = nxt; nxt = t; - } - } - }; - buildOp(acc, chunkBytes * 8ull); - if (lastChunkBytes != chunkBytes) buildOp(accLast, lastChunkBytes * 8ull); uint32_t crc = d_chunkCRC[0]; for (uint64_t i = 1; i < chunkCount; ++i) { - const uint32_t *op = (i + 1 == chunkCount && lastChunkBytes != chunkBytes) ? accLast : acc; + const uint32_t *op = (i + 1 == chunkCount) ? d_accLast : d_acc; crc = crc32Gf2MatTimes(op, crc) ^ d_chunkCRC[i]; } *d_crc = crc; @@ -226,7 +224,17 @@ inline void blockedCRC32(const void *d_data, size_t size, const uint32_t *d_lut, crc32SlicedKernel<<>>( d_checksums, d_chunkCRC, chunkCount, log2ChunkSize, checksumBytes, d_lut); cudaCheckError(); - crc32CombineKernel<<<1, 1, 0, stream>>>(d_chunkCRC, chunkCount, chunkSize, lastChunkBytes, d_crc); + // Precompute the two GF(2) shift operators on the host (data-independent, + // ~microseconds) and upload them, instead of rebuilding them on a single + // device thread inside the combine - that build dominated the checksum + // for small and medium grids. + uint32_t hOps[64]; + crc32BuildShiftOp(hOps, chunkSize * 8ull);// operator for a full chunk + crc32BuildShiftOp(hOps + 32, lastChunkBytes * 8ull);// operator for the last chunk + unique_ptr opBuffer(64, stream); + uint32_t *d_ops = opBuffer.get(); + cudaCheck(cudaMemcpyAsync(d_ops, hOps, 64*sizeof(uint32_t), cudaMemcpyHostToDevice, stream)); + crc32CombineKernel<<<1, 1, 0, stream>>>(d_chunkCRC, chunkCount, d_ops, d_ops + 32, d_crc); cudaCheckError(); } }// void cudaBlockedCRC32(const void *d_data, size_t size, const uint32_t *d_lut, uint32_t *d_crc, cudaStream_t stream) diff --git a/pendingchanges/nanovdbchecksumcombine.txt b/pendingchanges/nanovdbchecksumcombine.txt new file mode 100644 index 0000000000..1d864cae66 --- /dev/null +++ b/pendingchanges/nanovdbchecksumcombine.txt @@ -0,0 +1,5 @@ +NanoVDB: + Improvements: + - Sped up the CUDA GridChecksum combine: the per-chunk CRCs are now folded with + GF(2) shift operators precomputed on the host, instead of rebuilding them on a + single device thread on every call (~1.7x faster on small and medium grids). From 06ad21bf7069186ea19f90f2793f71c2c9e3c835 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Sun, 19 Jul 2026 03:53:42 +0000 Subject: [PATCH 6/8] Improve comments for processInternal/Leaf Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/GridStats.cuh | 46 ++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 3 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/GridStats.cuh b/nanovdb/nanovdb/tools/cuda/GridStats.cuh index 4ba6658f8c..be2ea9d326 100644 --- a/nanovdb/nanovdb/tools/cuda/GridStats.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridStats.cuh @@ -73,27 +73,47 @@ __global__ void processLeaf(NodeManager *d_nodeMgr, StatsT *d_stats) auto &d_leaf = d_nodeMgr->leaf(tid); bool nonEmpty = false; - if (lane == 0) nonEmpty = d_leaf.updateBBox();// updates active bounding box (also updates data->mFlags) - nonEmpty = __shfl_sync(0xffffffffu, nonEmpty, 0); + // updates leaf's active bounding box (also updates data->mFlags), only run on lane 0 + if (lane == 0) nonEmpty = d_leaf.updateBBox(); + // broadcast nonEmpty value to all lanes; register-to-register warp shuffle with all lanes participating + nonEmpty = __shfl_sync(0xffffffffu, nonEmpty, 0); // 0xffffffffu is "all 32 lanes" mask if (nonEmpty) { if constexpr(StatsT::hasStats()) { + // 1) Per-lane partial: each of the 32 lanes folds its strided share of the + // 512 voxel slots (lane L visits L, L+32, L+64, ...) into a local StatsT, + // skipping inactive voxels via the value mask. This is the coalesced pass - + // consecutive lanes touch consecutive slots on each stride. StatsT stats; const auto &mask = d_leaf.valueMask(); for (uint32_t i = lane; i < 512; i += 32) if (mask.isOn(i)) stats.add(d_leaf.getValue(i)); - StatsT *sWarp = sStats + (warpID << 5); + // 2) Warp reduction: stage the 32 partials in this warp's slice of shared + // memory, then merge them pairwise in log2(32)=5 steps (16->8->4->2->1). + // StatsT::add is an associative merge (exact for min/max; a Welford + // combine for mean/variance), so the tree order is safe. __syncwarp + // after each step orders the shared-memory writes within the warp. + StatsT *sWarp = sStats + (warpID << 5);// this warp's 32-slot scratch region sWarp[lane] = stats; __syncwarp(); + // pairwise reduction for (uint32_t d = 16; d; d >>= 1) { if (lane < d) sWarp[lane].add(sWarp[lane + d]); __syncwarp(); } + // 3) Publish: lane 0 now holds the leaf's fully merged stats in sWarp[0]. if (lane == 0) { if constexpr(StatsT::hasAverage()) { + // mean/variance: the parent lower node must Welford-merge this leaf's + // full accumulator (count+mean+M2), not just its extrema, so publish + // the whole StatsT into the per-node scratch array at this leaf's slot + // and stash that slot index in mMinimum (reused as a uint32 handle) so + // the parent kernel can locate it. d_stats[tid] = sWarp[0]; *reinterpret_cast(&d_leaf.mMinimum) = tid; } else { + // min/max only: extrema compose directly from child to parent, so write + // them straight into the leaf - no scratch slot needed. sWarp[0].setStats(d_leaf); } } @@ -120,27 +140,46 @@ __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) if (nodeID >= d_nodeMgr->nodeCount(LEVEL)) return; auto &d_node = d_nodeMgr->template node(nodeID); + // 1) Per-thread partial: each thread folds its strided share of the child table + // (NodeT::SIZE = 4096 for lower, 32768 for upper) into a local bbox + stats. + // Children were finalized by an earlier launch (leaves, then lower, then upper), + // so a child's bbox()/stats are valid to read here. Each entry is one of three: CoordBBox bbox;// empty StatsT stats; for (uint32_t i = tID; i < NodeT::SIZE; i += Threads) { if (d_node.childMask().isOn(i)) { + // (a) a child node: union its finalized bbox and merge its statistics. auto &child = *d_node.getChild(i); bbox.expand( child.bbox() ); if constexpr(StatsT::hasAverage()) { + // The child published its full accumulator into d_stats and stashed the + // slot index in its mMinimum. We must read that handle FIRST (to index + // d_stats), then setStats commits the child's real min/max/avg/std into + // the child node - which also overwrites the mMinimum-as-slot handle with + // the true minimum. Finally merge the child's accumulator into ours. StatsT &s = d_stats[*reinterpret_cast(&child.mMinimum)]; s.setStats(child); stats.add(s); } else if constexpr(StatsT::hasMinMax()) { + // min/max compose directly - the child already holds its final extrema. stats.add(child.minimum()); stats.add(child.maximum()); } } else if (d_node.valueMask().isOn(i)) { + // (b) an active tile: one constant value covering the whole child region + // (no child node). Grow the bbox to span the tile's extent, and add the + // value with multiplicity NUM_VALUES (the voxel count it stands in for) + // so the mean/variance are weighted correctly. const Coord ijk = d_node.offsetToGlobalCoord(i); bbox[0].minComponent(ijk); bbox[1].maxComponent(ijk + Coord(ChildT::DIM - 1)); if constexpr(StatsT::hasStats()) stats.add(d_node.data()->getValue(i), ChildT::NUM_VALUES); } + // (c) otherwise inactive - contributes nothing. } + // 2) Block reduction: merge the 128 per-thread partials (both stats and bbox) via a + // shared-memory tree in log2(128)=7 steps. This spans 4 warps, so it needs + // __syncthreads (block-wide barrier), not the __syncwarp used in processLeaf. sStats[tID] = stats; sBBox[tID] = bbox; __syncthreads(); @@ -151,6 +190,7 @@ __global__ void processInternal(NodeManager *d_nodeMgr, StatsT *d_stats) } __syncthreads(); } + // 3) Publish: thread 0 now holds the node's merged bbox + stats. if (tID == 0) { d_node.mBBox = sBBox[0]; if constexpr(StatsT::hasAverage()) { From ad251341d03d22e58cadc6acf4806286a676036c Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Mon, 20 Jul 2026 00:28:00 +0000 Subject: [PATCH 7/8] Comment cleanup and swap some constants for semantically defined values Signed-off-by: Jonathan Swartz --- nanovdb/nanovdb/tools/cuda/GridStats.cuh | 2 +- nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh | 15 +++++---------- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/nanovdb/nanovdb/tools/cuda/GridStats.cuh b/nanovdb/nanovdb/tools/cuda/GridStats.cuh index be2ea9d326..f2c85f89d5 100644 --- a/nanovdb/nanovdb/tools/cuda/GridStats.cuh +++ b/nanovdb/nanovdb/tools/cuda/GridStats.cuh @@ -86,7 +86,7 @@ __global__ void processLeaf(NodeManager *d_nodeMgr, StatsT *d_stats) // consecutive lanes touch consecutive slots on each stride. StatsT stats; const auto &mask = d_leaf.valueMask(); - for (uint32_t i = lane; i < 512; i += 32) + for (uint32_t i = lane; i < NanoLeaf::SIZE; i += 32) if (mask.isOn(i)) stats.add(d_leaf.getValue(i)); // 2) Warp reduction: stage the 32 partials in this warp's slice of shared // memory, then merge them pairwise in log2(32)=5 steps (16->8->4->2->1). diff --git a/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh b/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh index 63d5151339..31cc0aa507 100644 --- a/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh +++ b/nanovdb/nanovdb/tools/cuda/IndexToGrid.cuh @@ -233,15 +233,12 @@ __global__ void processNodesKernel(typename IndexToGrid::NodeAccessor if (srcGrid.hasStdDeviation()) dstNode.mStdDevi = srcValues[srcNode.mStdDevi]; } } - // Cooperative, coalesced mask copies (upper-node masks are 4 KB each; a - // single thread used to copy them member-wise while the block idled) - for (int w = tid; w < SrcNodeT::SIZE/64; w += nThreads) { + // Cooperative, coalesced mask copies + for (int w = tid; w < srcNode.mValueMask.wordCount(); w += nThreads) { dstNode.mValueMask.words()[w] = srcNode.mValueMask.words()[w]; dstNode.mChildMask.words()[w] = srcNode.mChildMask.words()[w]; } - // Consecutive threads process consecutive table entries (coalesced); the - // former mapping gave each thread a consecutive RUN, i.e. a 32-entry - // stride across the warp on every iteration. + // Consecutive threads process consecutive table entries (coalesced) for (int i = tid; i < SrcNodeT::SIZE; i += nThreads) { if (srcNode.mChildMask.isOn(i)) { if constexpr(sizeof(SrcNodeT)==sizeof(DstNodeT) && sizeof(SrcChildT)==sizeof(DstChildT)) { @@ -293,10 +290,8 @@ __global__ void processLeafsKernel(typename IndexToGrid::NodeAccessor if (srcGrid.hasStdDeviation()) dstLeaf.mStdDevi = srcValues[srcLeaf.getDev()]; } } - // Consecutive threads write consecutive values (coalesced); the former - // mapping gave each thread a consecutive run of 8, i.e. an 8-entry - // stride across the warp on every iteration. - for (int i = tid; i < 512; i += nThreads) + // Consecutive threads write consecutive values (coalesced) + for (int i = tid; i < NanoLeaf::SIZE; i += nThreads) dstLeaf.mValues[i] = srcValues[srcLeaf.getValue(i)]; }// processLeafsKernel From 65c9d4884c195d2a75c65b858484fdaf5c75364d Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Mon, 20 Jul 2026 00:38:44 +0000 Subject: [PATCH 8/8] Update pendingchanges Signed-off-by: Jonathan Swartz --- pendingchanges/nanovdbchecksumcombine.txt | 5 ----- pendingchanges/nanovdboptimizations.txt | 9 +++++++++ 2 files changed, 9 insertions(+), 5 deletions(-) delete mode 100644 pendingchanges/nanovdbchecksumcombine.txt create mode 100644 pendingchanges/nanovdboptimizations.txt diff --git a/pendingchanges/nanovdbchecksumcombine.txt b/pendingchanges/nanovdbchecksumcombine.txt deleted file mode 100644 index 1d864cae66..0000000000 --- a/pendingchanges/nanovdbchecksumcombine.txt +++ /dev/null @@ -1,5 +0,0 @@ -NanoVDB: - Improvements: - - Sped up the CUDA GridChecksum combine: the per-chunk CRCs are now folded with - GF(2) shift operators precomputed on the host, instead of rebuilding them on a - single device thread on every call (~1.7x faster on small and medium grids). diff --git a/pendingchanges/nanovdboptimizations.txt b/pendingchanges/nanovdboptimizations.txt new file mode 100644 index 0000000000..d4a810db60 --- /dev/null +++ b/pendingchanges/nanovdboptimizations.txt @@ -0,0 +1,9 @@ +NanoVDB: + Improvements: + - Parallelized the CUDA GridChecksum CRC32 with a slicing-by-4 table lookup and + an associative GF(2) combine (fold operators precomputed on the host); results + are bit-identical and measured 17x faster on large grids. + - Reworked the CUDA GridStats into cooperative reductions (warp-per-leaf and + block-per-internal-node), replacing the one-thread-per-node kernel. + - Made the CUDA IndexToGrid node and leaf remapping cooperative and coalesced, so + masks and tables are copied by the whole block rather than a single thread.