   x = prefix[DIM - 1];
   __syncthreads();
   return x;
-#else
-
-  int offset = 0;
-  int tid2 = tid << 1;
-
-#pragma unroll
-  for (int d = DIM >> 1; d > 0; d >>= 1) {
-    __syncthreads();
-
-    int iflag = BTEST(tid < d);
-    int ai = (((tid2 + 1) << offset) - 1) & iflag;
-    int bi = (((tid2 + 2) << offset) - 1) & iflag;
-
-    prefix[bi] += prefix[ai] & iflag;
-    offset++;
-  }
-
-  // clear the last element
-  if (tid == 0) prefix[DIM - 1] = 0;
-
-  // traverse down the tree building the scan in place
-#pragma unroll
-  for (int d = 1; d < DIM; d <<= 1) {
-    offset--;
-    __syncthreads();
-
-    int iflag = BTEST(tid < d);
-    int ai = (((tid2 + 1) << offset) - 1) & iflag;
-    int bi = (((tid2 + 2) << offset) - 1) & iflag;
-
-    int t       = prefix[ai];
-    if (tid < d) {
-      prefix[ai]  = (prefix[bi] & iflag) + (t & BTEST(tid >= d));
-      prefix[bi] += t & iflag;
-    }
-  }
-  __syncthreads();
-
-  prefix[tid] += value;
-  __syncthreads();
-
-  x = prefix[DIM - 1];
-  __syncthreads();
-  return x;
-#endif
 }
 
-#if 1
+/****** KEPLER __shfl prefix sum ******/
+
 __device__ __forceinline__ uint shfl_scan_add_step(uint partial, uint up_offset)
 {
   uint result;

**** --> prefix calculation via Horn(2005) data-parallel algoritm
***/
#define BTEST(x) (-(int)(x))
template<int DIM2>
__device__ int calc_prefix(int N, int* prefix_in, int tid) {
  int x, y = 0;

  const int DIM = 1 << DIM2;
  
  for (int p = 0; p < N; p += DIM) {
    int *prefix = &prefix_in[p];

    x = prefix[tid -  1]; __syncthreads(); prefix[tid] += x & BTEST(tid >=  1); __syncthreads();
    x = prefix[tid -  2]; __syncthreads(); prefix[tid] += x & BTEST(tid >=  2); __syncthreads();
    x = prefix[tid -  4]; __syncthreads(); prefix[tid] += x & BTEST(tid >=  4); __syncthreads();
    x = prefix[tid -  8]; __syncthreads(); prefix[tid] += x & BTEST(tid >=  8); __syncthreads();
    x = prefix[tid - 16]; __syncthreads(); prefix[tid] += x & BTEST(tid >= 16); __syncthreads();
    if (DIM2 >= 6) {x = prefix[tid - 32]; __syncthreads(); prefix[tid] += x & BTEST(tid >= 32); __syncthreads();}
    if (DIM2 >= 7) {x = prefix[tid - 64]; __syncthreads(); prefix[tid] += x & BTEST(tid >= 64); __syncthreads();}
    if (DIM2 >= 8) {x = prefix[tid -128]; __syncthreads(); prefix[tid] += x & BTEST(tid >=128); __syncthreads();}
    

    prefix[tid] += y;
    __syncthreads();

    y = prefix[DIM-1];
    __syncthreads();
  }

  return y;
} 

template<int DIM2>
__device__ int calc_prefix1(int* prefix, int tid, int value)
{
  int  x;

        if (flag && nodesM[offset + 6] == 0) nodesM[offset + 6] = child + 6;
        if (flag && nodesM[offset + 7] == 0) nodesM[offset + 7] = child + 7;
        __syncthreads();

        n_offset += n_total;    //Increase the offset in the array by the number of newly added nodes

