ex_*stencil_cpu: add legacy-transposed pass; BatchAccessor: document §8l

Added a tap-outer, voxel-inner variant of the Legacy path (same
leaf-only ReadAccessor, same probeLeaf + getValue mechanics, just the
nested loops swapped) as a new `legacy-transposed` benchmark pass in
both examples.  Checksums match the voxel-outer `legacy` pass byte-
for-byte on both synthetic and narrow-band workloads.

During the experiment we hit a GCC inlining pitfall: a runtime-args
inner lambda `[&](int di, int dj, int dk)` invoked 18 times via
parameter-pack fold did *not* get inlined (lambda body contains a
128-iteration loop; GCC's inline budget × 18 is exhausted).  Result:
18 explicit call instructions to a 542-byte processTap function with
6-register prologue/epilogue per call, plus tap offsets becoming
runtime register arguments (one spilled to stack) — accounting for
~13 % of the observed slowdown vs. Legacy.  Fix is a templated lambda
`[&]<int DI, int DJ, int DK>() [[gnu::always_inline]]` dispatched via
`.template operator()<...>()`.  Standalone processTap symbol vanishes;
transposed body grows 4.4 → 9.8 KB, matching Legacy's 10.5 KB.

Measured at ~32M active voxels on i9-285K (24 threads):
  - Narrowband taperLER: Legacy 2.2 ns/vox vs Transposed 2.1 ns/vox
    (marginal, within the ~10 % noise floor)
  - Synthetic 64M/50%:   Legacy 2.4 ns/vox vs Transposed 2.8 ns/vox
    (+19 %, outside noise)

Implementation verdict: LegacyStencilAccessor's voxel-outer moveTo
stays the default.  Tap-outer has no consistent perf advantage, and
voxel-outer wins on cleanliness (self-contained accessor, no scratch
arrays, no compiler-inlining fragility, 1:1 mapping to the stencil-
operator mental model).  `legacy-transposed` kept as a benchmark pass
for future reference.

BatchAccessor.md §8l captures the experiment, the inlining-pitfall
lesson, the measurement matrix, and the implementation-quality
rationale behind keeping voxel-outer as the production default.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Signed-off-by: Efty Sifakis <esifakis@nvidia.com>

Author: sifakis

nanovdb/nanovdb/examples/ex_narrowband_stencil_cpu/narrowband_stencil_cpu.cpp

@@ -356,7 +356,8 @@ static void runPerf(
     const std::string&                                                    passFilter = "all")
 {
     // wantPass(<name>) returns true if this pass should run under the current filter.
-    // Supported names: "decode", "stencil", "framing", "legacy".  "all" runs everything.
+    // Supported names: "decode", "stencil", "framing", "legacy",
+    //                  "legacy-transposed".  "all" runs everything.
     auto wantPass = [&](const char* name) {
         return passFilter == "all" || passFilter == name;
     };
@@ -539,6 +540,79 @@ static void runPerf(
                         [](uint64_t a, uint64_t b) { return a ^ b; });
     }  // end wantPass("legacy")
 
+    // ---- Legacy transposed: tap-outer, voxel-inner ----
+    // Same semantics as `legacy`, reordered.  For each of the 18 WENO5 taps,
+    // sweep all BlockWidth voxels — giving long runs of probeLeaf + getValue
+    // calls with the SAME compile-time tap offset but varying center voxels.
+    double legacyXposedUs = 0.0;
+    uint64_t legacyXposedChecksum = 0;
+    if (wantPass("legacy-transposed")) {
+    std::fill(sums.begin(), sums.end(), uint64_t(0));
+
+    using Weno5Taps = nanovdb::Weno5Stencil::Taps;
+    static constexpr int SIZE = int(std::tuple_size_v<Weno5Taps>);
+
+    legacyXposedUs = timeForEach([&] {
+        nanovdb::util::forEach(size_t(0), size_t(nBlocks), size_t(1),
+            [&](const nanovdb::util::Range1D& range) {
+                alignas(64) uint32_t leafIndex[BlockWidth];
+                alignas(64) uint16_t voxelOffset[BlockWidth];
+                alignas(64) nanovdb::Coord centers[BlockWidth];
+                alignas(64) uint64_t s[BlockWidth];
+                nanovdb::ReadAccessor<BuildT, 0, -1, -1> acc(grid->tree().root());
+                uint64_t* bs0 = sums.data();
+
+                for (size_t bID = range.begin(); bID != range.end(); ++bID) {
+                    CPUVBM::decodeInverseMaps(
+                        grid, firstLeafID[bID],
+                        &jumpMap[bID * CPUVBM::JumpMapLength],
+                        firstOffset + bID * BlockWidth,
+                        leafIndex, voxelOffset);
+
+                    for (int i = 0; i < BlockWidth; ++i) {
+                        s[i] = 0;
+                        if (leafIndex[i] == CPUVBM::UnusedLeafIndex) continue;
+                        const uint16_t vo = voxelOffset[i];
+                        const uint32_t li = leafIndex[i];
+                        const nanovdb::Coord cOrigin = firstLeaf[li].origin();
+                        centers[i] = cOrigin + nanovdb::Coord(
+                            (vo >> 6) & 7, (vo >> 3) & 7, vo & 7);
+                    }
+
+                    auto processTap = [&]<int DI, int DJ, int DK>()
+                        [[gnu::always_inline]]
+                    {
+                        for (int i = 0; i < BlockWidth; ++i) {
+                            if (leafIndex[i] == CPUVBM::UnusedLeafIndex) continue;
+                            const nanovdb::Coord c = centers[i]
+                                + nanovdb::Coord(DI, DJ, DK);
+                            const LeafT* leaf = acc.probeLeaf(c);
+                            if (!leaf) continue;
+                            const uint32_t offset = (uint32_t(c[0] & 7) << 6)
+                                                  | (uint32_t(c[1] & 7) << 3)
+                                                  |  uint32_t(c[2] & 7);
+                            s[i] += leaf->data()->getValue(offset);
+                        }
+                    };
+
+                    [&]<size_t... Is>(std::index_sequence<Is...>) {
+                        (processTap.template operator()<
+                            std::tuple_element_t<Is, Weno5Taps>::di,
+                            std::tuple_element_t<Is, Weno5Taps>::dj,
+                            std::tuple_element_t<Is, Weno5Taps>::dk>(), ...);
+                    }(std::make_index_sequence<SIZE>{});
+
+                    uint64_t* bs = bs0 + bID * BlockWidth;
+                    for (int i = 0; i < BlockWidth; ++i) bs[i] = s[i];
+                }
+            });
+    });
+
+    legacyXposedChecksum =
+        std::accumulate(sums.begin(), sums.end(), uint64_t(0),
+                        [](uint64_t a, uint64_t b) { return a ^ b; });
+    }  // end wantPass("legacy-transposed")
+
     std::printf("\nEnd-to-end stencil gather (%u blocks, %lu active voxels):\n",
         nBlocks, nVoxels);
     std::printf("  decodeInverseMaps only: %7.1f ms  (%5.1f ns/voxel)\n",
@@ -552,9 +626,14 @@ static void runPerf(
     std::printf("  LegacyStencilAccessor : %7.1f ms  (%5.1f ns/voxel)  [%+5.1f ms over decode]  checksum=0x%016lx\n",
         legacyUs  / 1e3, legacyUs  * 1e3 / double(nVoxels),
         (legacyUs - decodeUs) / 1e3, legacyChecksum);
+    std::printf("  Legacy transposed     : %7.1f ms  (%5.1f ns/voxel)  [%+5.1f ms over decode]  checksum=0x%016lx\n",
+        legacyXposedUs / 1e3, legacyXposedUs * 1e3 / double(nVoxels),
+        (legacyXposedUs - decodeUs) / 1e3, legacyXposedChecksum);
 
     if (stencilChecksum != legacyChecksum)
-        std::cerr << "  WARNING: checksums differ — accessor results disagree!\n";
+        std::cerr << "  WARNING: stencil/legacy checksums differ — accessor results disagree!\n";
+    if (legacyChecksum != legacyXposedChecksum)
+        std::cerr << "  WARNING: legacy/legacy-transposed checksums differ — ordering bug!\n";
 }
 
 // ============================================================
@@ -574,7 +653,7 @@ static void printUsage(const char* argv0)
         << "  --grid=<name>        Select grid by name (default: first FloatGrid)\n"
         << "  --pass=<name>        Run one perf pass:\n"
         << "                         all (default), verify, decode, stencil,\n"
-        << "                         framing, legacy\n"
+        << "                         framing, legacy, legacy-transposed\n"
         << "  --threads=<n>        Limit TBB parallelism (0 = TBB default)\n"
         << "  --skip-validation    Skip the sidecar ordering sanity check\n";
 }

nanovdb/nanovdb/examples/ex_stencil_gather_cpu/stencil_gather_cpu.cpp

@@ -262,7 +262,8 @@ static void runPerf(
     const std::string&                                                    passFilter = "all")
 {
     // wantPass(<name>) returns true if this pass should run under the current filter.
-    // Supported names: "decode", "stencil", "framing", "legacy".  "all" runs everything.
+    // Supported names: "decode", "stencil", "framing", "legacy",
+    //                  "legacy-transposed".  "all" runs everything.
     auto wantPass = [&](const char* name) {
         return passFilter == "all" || passFilter == name;
     };
@@ -445,6 +446,79 @@ static void runPerf(
                         [](uint64_t a, uint64_t b) { return a ^ b; });
     }  // end wantPass("legacy")
 
+    // ---- Legacy transposed: tap-outer, voxel-inner ----
+    // Same semantics as `legacy`, reordered.  For each of the 18 WENO5 taps,
+    // sweep all BlockWidth voxels — giving long runs of probeLeaf + getValue
+    // calls with the SAME compile-time tap offset but varying center voxels.
+    double legacyXposedUs = 0.0;
+    uint64_t legacyXposedChecksum = 0;
+    if (wantPass("legacy-transposed")) {
+    std::fill(sums.begin(), sums.end(), uint64_t(0));
+
+    using Weno5Taps = nanovdb::Weno5Stencil::Taps;
+    static constexpr int SIZE = int(std::tuple_size_v<Weno5Taps>);
+
+    legacyXposedUs = timeForEach([&] {
+        nanovdb::util::forEach(size_t(0), size_t(nBlocks), size_t(1),
+            [&](const nanovdb::util::Range1D& range) {
+                alignas(64) uint32_t leafIndex[BlockWidth];
+                alignas(64) uint16_t voxelOffset[BlockWidth];
+                alignas(64) nanovdb::Coord centers[BlockWidth];
+                alignas(64) uint64_t s[BlockWidth];
+                nanovdb::ReadAccessor<BuildT, 0, -1, -1> acc(grid->tree().root());
+                uint64_t* bs0 = sums.data();
+
+                for (size_t bID = range.begin(); bID != range.end(); ++bID) {
+                    CPUVBM::decodeInverseMaps(
+                        grid, firstLeafID[bID],
+                        &jumpMap[bID * CPUVBM::JumpMapLength],
+                        firstOffset + bID * BlockWidth,
+                        leafIndex, voxelOffset);
+
+                    for (int i = 0; i < BlockWidth; ++i) {
+                        s[i] = 0;
+                        if (leafIndex[i] == CPUVBM::UnusedLeafIndex) continue;
+                        const uint16_t vo = voxelOffset[i];
+                        const uint32_t li = leafIndex[i];
+                        const nanovdb::Coord cOrigin = firstLeaf[li].origin();
+                        centers[i] = cOrigin + nanovdb::Coord(
+                            (vo >> 6) & 7, (vo >> 3) & 7, vo & 7);
+                    }
+
+                    auto processTap = [&]<int DI, int DJ, int DK>()
+                        [[gnu::always_inline]]
+                    {
+                        for (int i = 0; i < BlockWidth; ++i) {
+                            if (leafIndex[i] == CPUVBM::UnusedLeafIndex) continue;
+                            const nanovdb::Coord c = centers[i]
+                                + nanovdb::Coord(DI, DJ, DK);
+                            const LeafT* leaf = acc.probeLeaf(c);
+                            if (!leaf) continue;
+                            const uint32_t offset = (uint32_t(c[0] & 7) << 6)
+                                                  | (uint32_t(c[1] & 7) << 3)
+                                                  |  uint32_t(c[2] & 7);
+                            s[i] += leaf->data()->getValue(offset);
+                        }
+                    };
+
+                    [&]<size_t... Is>(std::index_sequence<Is...>) {
+                        (processTap.template operator()<
+                            std::tuple_element_t<Is, Weno5Taps>::di,
+                            std::tuple_element_t<Is, Weno5Taps>::dj,
+                            std::tuple_element_t<Is, Weno5Taps>::dk>(), ...);
+                    }(std::make_index_sequence<SIZE>{});
+
+                    uint64_t* bs = bs0 + bID * BlockWidth;
+                    for (int i = 0; i < BlockWidth; ++i) bs[i] = s[i];
+                }
+            });
+    });
+
+    legacyXposedChecksum =
+        std::accumulate(sums.begin(), sums.end(), uint64_t(0),
+                        [](uint64_t a, uint64_t b) { return a ^ b; });
+    }  // end wantPass("legacy-transposed")
+
     std::printf("\nEnd-to-end stencil gather (%u blocks, %lu active voxels):\n",
         nBlocks, nVoxels);
     std::printf("  decodeInverseMaps only: %7.1f ms  (%5.1f ns/voxel)\n",
@@ -458,9 +532,14 @@ static void runPerf(
     std::printf("  LegacyStencilAccessor : %7.1f ms  (%5.1f ns/voxel)  [%+5.1f ms over decode]  checksum=0x%016lx\n",
         legacyUs  / 1e3, legacyUs  * 1e3 / double(nVoxels),
         (legacyUs - decodeUs) / 1e3, legacyChecksum);
+    std::printf("  Legacy transposed     : %7.1f ms  (%5.1f ns/voxel)  [%+5.1f ms over decode]  checksum=0x%016lx\n",
+        legacyXposedUs / 1e3, legacyXposedUs * 1e3 / double(nVoxels),
+        (legacyXposedUs - decodeUs) / 1e3, legacyXposedChecksum);
 
     if (stencilChecksum != legacyChecksum)
-        std::cerr << "  WARNING: checksums differ — accessor results disagree!\n";
+        std::cerr << "  WARNING: stencil/legacy checksums differ — accessor results disagree!\n";
+    if (legacyChecksum != legacyXposedChecksum)
+        std::cerr << "  WARNING: legacy/legacy-transposed checksums differ — ordering bug!\n";
 }
 
 // ============================================================

nanovdb/nanovdb/util/BatchAccessor.md

@@ -1467,6 +1467,71 @@ follow-ons implied by this work but not pursued here:
   legacy compatibility shim (to be retired after a deprecation window)
   rather than a permanent toggle.
 
+### 8l. Follow-up: tap-outer loop ordering in the Legacy path
+
+Tested whether flipping `legacy`'s loop nest to tap-outer,
+voxel-inner helps on spatially-coherent workloads where many voxels in
+a batch are likely to share the same `valueMask` word.  Added a
+`legacy-transposed` benchmark pass in both `ex_stencil_gather_cpu` and
+`ex_narrowband_stencil_cpu`; checksums match `legacy` byte-for-byte on
+both workloads.
+
+#### 8l.1 Inlining pitfall
+
+First attempt used a runtime-args inner lambda
+`[&](int di, int dj, int dk) { for (int i = 0; i < 128; ++i) ... }`
+invoked 18 times via a parameter-pack fold.  GCC refused to inline the
+18 instantiations — the lambda body contains a 128-iteration loop with
+`probeLeaf` + `getValue` inside, which blew past the per-caller inline
+budget × 18.  Result: 18 explicit `call` instructions to a 542-byte
+`processTap` function with a 6-register prologue/epilogue per call,
+and tap offsets `(di, dj, dk)` as runtime register arguments (one
+spilled to stack) — so the compiler also couldn't specialise the loop
+body per tap.  That alone accounted for ~10 ms (~13 %) of the observed
+slowdown vs. Legacy.
+
+Fix: templated lambda
+`[&]<int DI, int DJ, int DK>() [[gnu::always_inline]] { ... }`
+dispatched via `.template operator()<di, dj, dk>()` inside the fold.
+The standalone `processTap` symbol disappears; transposed body grows
+from 4.4 KB → 9.8 KB (matching Legacy's 10.5 KB), and only cold-path
+tree-walk helpers remain as call targets.
+
+#### 8l.2 Results
+
+Measured at ~32M active voxel scale on i9-285K (24 threads, no HT):
+
+| Workload | Legacy (voxel-outer) | Transposed (tap-outer) | Δ |
+|----------|---------------------:|-----------------------:|--:|
+| Narrowband taperLER.vdb | 2.2 ns/vox | 2.1 ns/vox | −3 to −6 % (within noise) |
+| Synthetic 64M/50%       | 2.4 ns/vox | 2.8 ns/vox | +19 % |
+
+The narrowband tap-outer edge is marginal and within the ~10 %
+run-to-run noise floor observed on this host.  Synthetic's tap-outer
+slowdown is clearly outside noise.  Not a consistent win.
+
+#### 8l.3 Implementation verdict: voxel-outer stays the default
+
+`LegacyStencilAccessor`'s voxel-outer `moveTo(center)` kept as the
+production default:
+
+- **Clean abstraction**: `moveTo(center)` + indexed tap access maps
+  1:1 to the stencil operator's mental model.  A tap-outer batched
+  form would need external accumulator state and a centers-array
+  input, with no natural class boundary.
+- **No scratch arrays**: voxel-outer keeps the per-voxel accumulator
+  in a register and the 18-tap buffer inside the accessor; tap-outer
+  needs stack-local `centers[128]` and `s[128]`.
+- **Compiler robustness**: voxel-outer's 18-call single source
+  location is reliably collapsed by GCC.  Tap-outer relies on an
+  explicit `[[gnu::always_inline]]` workaround that, if lost during
+  future refactors, would silently regress performance by ~13 %.
+
+`legacy-transposed` retained as a benchmark pass for reference and as
+a datapoint reinforcing why the hybrid `StencilAccessor` is structured
+tap-outer (SIMD direction-computation inherently amortises across
+lanes at the same tap).
+
 ---
 
 ## 9. Relationship to Phase 1 Prototype
@@ -1559,6 +1624,16 @@ follow-ons implied by this work but not pursued here:
   Scope: benchmark-only; the library default is unchanged (right default
   for `probeValue`/`probeLeaf`/mixed workloads).
 
+- **Tap-outer loop ordering evaluation in the Legacy path (§8l)**:
+  added `legacy-transposed` benchmark pass (checksums match byte-for-byte)
+  and tested on both workloads at matched ~32M-voxel scale.  Narrowband:
+  marginal tap-outer edge, within noise.  Synthetic: tap-outer ~19 %
+  slower.  Uncovered a GCC inlining pitfall for runtime-args inner lambdas
+  (fixed via templated lambda + `[[gnu::always_inline]]`).
+  **Verdict**: voxel-outer `LegacyStencilAccessor` remains the default —
+  cleaner abstraction, no scratch arrays, no compiler-inlining fragility,
+  and no consistent perf advantage to tap-outer.
+
 ### Remaining
 
 - **`[[gnu::always_inline]]` on `Simd.h` helpers** (§8f) vs