 
     }
     else {
-
-      reduce<<<static_cast<uint32_t>(levels[0]), blockSize, 0, stream>>>(scratch_d + offsets[0],
-                                                                         counts_d,
-                                                                         n);
-      for (size_t i = 1; i < levels.size(); i++) {
-        reduce<<<static_cast<uint32_t>(levels[i]), blockSize, 0, stream>>>(scratch_d + offsets[i],
-                                                                           scratch_d + offsets[i - 1],
-                                                                           levels[i - 1]);
-      }
-
+      uint32_t L = static_cast<uint32_t>(levels.size());
+        
+      // From input, populate level 0
+      reduce << <levels[0], blockSize, 0, stream >> > (scratch_d + offsets[0],
+                                                       counts_d,
+                                                       n);
+      // From level i-1, populate level i, up to including L-1.
+      for (uint32_t i = 1; i < L; i++) {
+        reduce << <levels[i], blockSize, 0, stream >> > (scratch_d + offsets[i],
+                                                         scratch_d + offsets[i - 1],
+                                                         levels[i - 1]);
+      }
+
+      // Run scan on last level L-1, and write off total sum to last element of output (offsets_d+N).
       auto l = levels.size() - 1;
       if (sum_d == nullptr) {
         scan<false, true, false><<<1, blockSize, 0, stream>>>(scratch_d + offsets[l],

                                                             scratch_d + offsets[l],
                                                             levels[l]);
      }


    }

                                     size_t N,
                                     cudaStream_t stream)
{
  if (N <= std::numeric_limits<uint32_t>::max()) {
    uint32_t n = static_cast<uint32_t>(N);

    std::vector<uint32_t> levels;
    calcLevelSizes(levels, n);

    std::vector<uint32_t> offsets;
    offsets.push_back(0);
    size_t alignment = 128 / sizeof(uint32_t);
    for (size_t i = 0; i < levels.size(); i++) {
      size_t levelSize = (levels[i] + alignment - 1) & ~alignment;
      offsets.push_back(offsets[i] + levelSize);
    }

    uint32_t blockSize = SCAN_WARPS * WARPSIZE;

    if (levels.empty()) {

