 
   prefix[tid] += value;
   __syncthreads();
-  
+
   x = prefix[DIM - 1];
   __syncthreads();
   return x;
 #endif
 }
 
-template<int SHIFT>
+#if 1
+__device__ __forceinline__ uint shfl_scan_add_step(uint partial, uint up_offset)
+{
+  uint result;
+  asm(
+      "{.reg .u32 r0;"
+      ".reg .pred p;"
+      "shfl.up.b32 r0|p, %1, %2, 0;"
+      "@p add.u32 r0, r0, %3;"
+      "mov.u32 %0, r0;}"
+      : "=r"(result) : "r"(partial), "r"(up_offset), "r"(partial));
+  return result;
+}
+
+  template <const int levels>
+__device__ __forceinline__ uint inclusive_scan_warp(int mysum)
+{
+  for(int i = 0; i < levels; ++i)
+    mysum = shfl_scan_add_step(mysum, 1 << i);
+  return mysum;
+}
+
+template<int DIM2>
+__device__ __forceinline__ int calc_prefix(int* prefix, int tid, int value) 
+{
+  if (DIM2 != 6)  /* should never be called */
+    return calc_prefix1<DIM2>(prefix, tid, value);
+  else
+  {
+    const int WARP_SIZE2 = 5;
+    const int WARP_SIZE  = 1 << WARP_SIZE2;
+    prefix[tid] = inclusive_scan_warp<WARP_SIZE2>(value);
+    __syncthreads();
+
+    prefix[tid] += prefix[WARP_SIZE - 1] & BTEST(tid >= WARP_SIZE);
+    __syncthreads();
+
+    const int x = prefix[(1 << DIM2)- 1];
+    __syncthreads();       /* must be here, otherwise the code crashes */
+
+    return x;
+  }
+}
+#endif
+
+
+  template<int SHIFT>
 __forceinline__ __device__ int ACCS(const int i)
 {
   return (i & ((LMEM_STACK_SIZE << SHIFT) - 1))*blockDim.x + threadIdx.x;

} 

template<int DIM2>
__device__ int calc_prefix(int* prefix, int tid, int value) {
  int  x;
  
  const int DIM = 1 << DIM2;

  prefix[tid] = value;

  //  begin,    end,   size
  // -----------------------
#ifdef _ORIG_SHMEM_

  int *approxS = (int*)&shmem  [     0];            //  0*DIM,  2*DIM,  2*DIM
  int *directS = (int*)&approxS[ 2*DIM];            //  2*DIM,  3*DIM,  1*DIM
  int *nodesS = (int*)&directS [   DIM];            //  3*DIM, 12*DIM,  9*DIM
  int *prefix = (int*)&nodesS  [9 *DIM];            // 12*DIM, 14*DIM,  2*DIM
  int *sh_body = &approxS[DIM];

  int *prefix0 = &prefix[  0];
  int *prefix1 = &prefix[DIM];

  const int NJMAX = DIM*2;
  int    *body_list = (int*   )&nodesS   [  DIM]; //  4*DIM,  6*DIM,  2*DIM
  float  *sh_mass   = (float* )&body_list[NJMAX]; //  6*DIM,  7*DIM,  1*DIM
  float3 *sh_pos    = (float3*)&sh_mass  [  DIM]; //  7*DIM, 10*DIM   3*DIM

  int *approxM = approxS;
  int *directM = directS;
  int * nodesM =  nodesS;
