     }
   }
 
+  __device__ __forceinline__ uint32_t bfe(uint32_t val, uint32_t o, uint32_t n)
+  {
+    uint32_t rv;
+    asm("bfe.u32 %0, %1, %2, %3;" : "=r"(rv) : "r"(val), "r"(o), "r"(n));
+    return rv;
+  }
+
+  __global__ __launch_bounds__(256) void reduceBase2b(uint32_t* __restrict__ hp2_d,
+                                                      uint32_t* __restrict__ sb2_d,
+                                                      const uint32_t n2,
+                                                      uint32_t* __restrict__ hp1_d,
+                                                      const uint32_t n1,
+                                                      const uint32_t* __restrict__ sb0_d,
+                                                      const uint32_t n0)
+  {
+    const uint32_t offset0_base = 32 * 160 * blockIdx.x + 4 * threadIdx.x;
+    const uint32_t offset1_base = 160 * blockIdx.x;
+
+    const uint32_t lane = threadIdx.x %  HP5_WARP_SIZE;
+    const uint32_t warp = threadIdx.x / HP5_WARP_SIZE;
+
+    __shared__ uint32_t sb1[160];
+
+    // Each warp reads 32*4=128 bytes
+
+    for (uint32_t i = 0; i < 5; i++) {
+      const uint32_t offset0 = offset0_base + 4 * 256 * i;
+
+      uint4 value;
+      if (offset0 + 3 < n0) {
+        asm volatile("ld.global.cs.v4.u32 {%0, %1, %2, %3}, [%4];" : "=r"(value.x), "=r"(value.y), "=r"(value.z), "=r"(value.w) : PTR_REG(sb0_d + offset0));
+      }
+      else {
+        if (offset0 < n0) {
+          asm volatile("ld.global.cs.u32 %0, [%1];" : "=r"(value.x) : PTR_REG(sb0_d + offset0));
+          if (offset0 + 1 < n0) {
+            asm volatile("ld.global.cs.u32 %0, [%1];" : "=r"(value.y) : PTR_REG(sb0_d + offset0 + 1));
+            if (offset0 + 2 < n0) {
+              asm volatile("ld.global.cs.u32 %0, [%1];" : "=r"(value.z) : PTR_REG(sb0_d + offset0 + 2));
+            }
+            else {
+              value.y = 0;
+            }
+          }
+          else {
+            value.y = 0;
+          }
+        }
+        else {
+          value.x = 0;
+        }
+        value.z = 0;
+      }
+
+      // Prefetch data for next row into L2
+      if (i + 1 < 5) {
+        asm volatile("prefetch.global.L2 [%0];"::PTR_REG(sb0_d + offset0 + 4 * 256));
+      }
+
+      const uint32_t warpMaskA = __ballot(value.x != 0);
+      const uint32_t warpMaskB = __ballot(value.y != 0);
+      const uint32_t warpMaskC = __ballot(value.z != 0);
+      const uint32_t warpMaskD = __ballot(value.w != 0);
+
+      if (lane < 4) {
+        uint32_t warpMaskASub = bfe(warpMaskA, 8 * lane, 8);
+        uint32_t warpMaskBSub = bfe(warpMaskB, 8 * lane, 8);
+        uint32_t warpMaskCSub = bfe(warpMaskC, 8 * lane, 8);
+        uint32_t warpMaskDSub = bfe(warpMaskD, 8 * lane, 8);
+
+        uint32_t warpMask = warpMaskASub<<24 | warpMaskBSub<<16 | warpMaskCSub<<8 | warpMaskDSub;
+
+        const uint32_t offset1 = offset1_base + 32 * i + 4 * warp + lane;
+        if (offset1 < n1) {
+          asm volatile("st.global.cs.u32 [%0], %1;"::PTR_REG(hp1_d + offset1), "r"(warpMask));
+        }
+        sb1[32 * i + 4 * warp + lane] = __popc(warpMask);
+      }
+    }
+
+    __syncthreads();
+    if (threadIdx.x < HP5_WARP_SIZE) { // First warp
+      const uint32_t offset2 = 32 * blockIdx.x + threadIdx.x;
+      if (offset2 < n2) {
+        uint4 hp = make_uint4(sb1[5 * threadIdx.x + 0],
+                              sb1[5 * threadIdx.x + 1],
+                              sb1[5 * threadIdx.x + 2],
+                              sb1[5 * threadIdx.x + 3]);
+#if 1
+        // Streaming store
+        asm volatile("st.global.cs.v4.u32 [%0], {%1,%2,%3,%4};"
+                     ::PTR_REG((uint4*)hp2_d + offset2), "r"(hp.x), "r"(hp.y), "r"(hp.z), "r"(hp.w));
+#else
+        ((uint4*)hp2_d)[offset2] = hp;
+#endif
+        sb2_d[offset2] = hp.x + hp.y + hp.z + hp.w + sb1[5 * threadIdx.x + 4];
+      }
+    }
+
+  }
+
   __global__ __launch_bounds__(256) void reduceBase2(uint32_t* __restrict__ hp2_d,
                                                      uint32_t* __restrict__ sb2_d,
                                                      const uint32_t n2,

#if 1
      // Prefetch data for next row into L2
      if (i + 1 < 20) {
        asm volatile("prefetch.global.L2 [%0];"::PTR_REG(sb0_d + offset0 + 256));
      }
#endif



    size_t i = 0;
    if (1 < L) {
      ::reduceBase2b<<<(levels[1] + 31) / 32, 8 * 32, 0, stream>>>(scratch_d + offsets[1],
                                                                  scratch_d + offsets[L + 1 + (sb ? 1 : 0)],
                                                                  levels[1],
                                                                  scratch_d + offsets[0],
