     if (10 < L) offsets[10] = offsetBlob.dataA;
     if (11 < L) offsets[11] = offsetBlob.dataB;
 
-    // Traverse apex.
-    uint32_t offset = 0;
-    uint32_t key = index;
-    offset = processHistoElement(key, 5 * offset, hp_d[1]);
-    offset = processHistoElement(key, 5 * offset, hp_d[2 + offset]);
-    offset = processHistoElement(key, 5 * offset, hp_d[7 + offset]);
-    for (uint32_t i = L; 0 < i; i--) {
-      offset = processDataElement(key, 5 * offset, hp_d[offsets[i - 1] + offset]);
-    }
-    out_d[index] = offset;
+    uint32_t N = hp_d[0].x;
+    uint4 T = hp_d[1];
+
+    for (uint32_t k = HP5_EXTRACT_BLOCKSIZE * blockIdx.x; k < N; k += HP5_EXTRACT_BLOCKS * HP5_EXTRACT_BLOCKSIZE) {
+      const uint32_t index = k + threadIdx.x;
+      if (index < N) {
+        uint32_t offset = 0;
+        uint32_t key = index;
+        offset = processHistoElement(key, 5 * offset, T);
+        offset = processHistoElement(key, 5 * offset, hp_d[2 + offset]);
+        offset = processHistoElement(key, 5 * offset, hp_d[7 + offset]);
+        for (uint32_t i = L; 0 < i; i--) {
+          offset = processDataElement(key, 5 * offset, hp_d[offsets[i - 1] + offset]);
+        }
+        out_d[index] = offset;
+      }
+    }
   }
 
   void scratchLayout(std::vector<uint32_t>& levels, std::vector<uint32_t>& offsets, uint32_t N)

  }


  OffsetBlob offsetBlob;
  std::memcpy(&offsetBlob.data0, offsets.data(), sizeof(uint32_t)*std::min(size_t(12), L));

  switch (L)
  {
  case 0: ::extract<0><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 1: ::extract<1><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 2: ::extract<2><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 3: ::extract<3><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 4: ::extract<4><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 5: ::extract<5><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 6: ::extract<6><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 7: ::extract<7><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 8: ::extract<8><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 9: ::extract<9><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 10: ::extract<10><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  case 11: ::extract<11><<<(N + 127) / 128, 128, 0, stream>>>(out_d, reinterpret_cast<uint4*>(scratch_d), offsetBlob); break;
  default:
    abort();
    break;

               const uint4* __restrict__ hp_d,
               OffsetBlob offsetBlob)
  {

    __shared__ uint32_t offsets[12];
    if (0 < L) offsets[0] = offsetBlob.data0;
