                                       threshold);
   }
 
-  
+  template<class FieldType>
+  __global__ __launch_bounds__(32 * 6) void reduceBaseIndexedMultiWarp(uint8_t* __restrict__           index_cases_d,
+                                                                       uchar4* __restrict__            out_vertex_level0_d,
+                                                                       uint32_t* __restrict__          out_vertex_sideband_d,
+                                                                       uchar4* __restrict__            out_index_level0_d,
+                                                                       uint32_t* __restrict__          out_index_sideband_d,
+                                                                       const uint8_t* __restrict__     index_count,
+                                                                       const FieldType* __restrict__   field_d,
+                                                                       const size_t                    field_row_stride,
+                                                                       const size_t                    field_slice_stride,
+                                                                       const float                     threshold,
+                                                                       const uint3                     chunks,
+                                                                       const uint3                     grid_max_index)
+  {
+    const uint32_t warp = threadIdx.x / 32;
+    const uint32_t thread = threadIdx.x % 32;
+    const uint32_t chunk_ix = blockIdx.x;
+
+    // Figure out which chunk we are in.
+    // FIXME: Try to use grid extents to mach chunk grid to avoid all these modulo/divisions.
+    uint3 chunk = make_uint3(chunk_ix % chunks.x,
+                             (chunk_ix / chunks.x) % chunks.y,
+                             (chunk_ix / chunks.x) / chunks.y);
+    uint3 cell = make_uint3(31 * chunk.x + thread,
+                            5 * chunk.y,
+                            5 * chunk.z);
+
+    __shared__ float shmem[32 * 6 * 2];
+
+    reduceBaseKernelMultiWarp<FieldType, true>(index_cases_d + static_cast<size_t>(5 * 32 * 5 * blockIdx.x + 5 * thread), // Index doesn't need more than 32 bits
+                                               out_index_level0_d + static_cast<size_t>(32 * 5 * blockIdx.x + thread),
+                                               out_vertex_level0_d + static_cast<size_t>(32 * 5 * blockIdx.x + thread),
+                                               out_index_sideband_d + static_cast<size_t>(32 * 5 * blockIdx.x + thread),
+                                               out_vertex_sideband_d + static_cast<size_t>(32 * 5 * blockIdx.x + thread),
+                                               shmem,
+                                               index_count,
+                                               field_d + cell.z * field_slice_stride
+                                               + cell.y * field_row_stride
+                                               + min(grid_max_index.x, cell.x),
+                                               field_row_stride,
+                                               field_slice_stride,
+                                               cell,
+                                               grid_max_index,
+                                               warp,
+                                               thread,
+                                               threshold);
+  }
 
 
 

                                                                           make_uint3(ctx->grid_size.x - 1,
                                                                                      ctx->grid_size.y - 1,
                                                                                      ctx->grid_size.z - 1));

      CHECKED_CUDA(cudaEventRecord(ctx->baseEvent, stream));
      CHECKED_CUDA(cudaStreamWaitEvent(ctx->indexStream, ctx->baseEvent, 0));

      //        case it might be better to just merge the two launches and use e.g. block.y to tell which
      //        pyramid to reduce.
      bool sb = true;
      for (unsigned i = i0; i < ctx->levels - 3; i++) {
        const auto blocks = (ctx->level_sizes[i] + 31) / 32;
        reduce1 << <blocks, 5 * 32, 0, stream >> > (ctx->index_pyramid + ctx->level_offsets[i],    // Each block will write 32 uvec4's into this
                                                    ctx->index_sidebands[sb ? 1 : 0],              // Each block will write 32 uint32's into this
