 
 #endif
 
-  template<class OP, class T, const int BLOCKDIM2>
-__device__ __inline__ T inclusive_scan_block(T v_in, const unsigned int idx) 
+  template<const int BLOCKDIM2>
+__device__ __inline__ int inclusive_scan_block(volatile int* shdata, int v_in, const unsigned int idx) 
 {
-  __shared__ T shdata[WARP_SIZE];
   const unsigned int laneId = idx & (WARP_SIZE - 1);
   const unsigned int warpId = idx >> WARP_SIZE2;
 
-  T val = inclusive_scan_warp<OP, T, WARP_SIZE2>(v_in, idx);
+  int val = inclusive_scan_warp<WARP_SIZE2>(v_in, idx);
   if (31 == laneId) shdata[warpId] = val;
   __syncthreads();
 
-  if (0 == warpId) shdata[idx] = inclusive_scan_warp<OP, T, BLOCKDIM2 - WARP_SIZE2>(shdata[idx], idx);
-  __syncthreads();
-
-  if (warpId > 0) val = OP::apply(shdata[warpId - 1], val);
+  if (0 == warpId) shdata[idx] = inclusive_scan_warp<BLOCKDIM2 - WARP_SIZE2>(shdata[idx], idx);
+  __syncthreads();
+
+  if (warpId > 0) val += shdata[warpId - 1];
 
   return val; 
 }
 
+  template<const int BLOCKDIM2>
+__device__ __inline__ int2 inclusive_scan_blockS(volatile int* shdata, int v_in, const unsigned int idx) 
+{
+  const unsigned int laneId = idx & (WARP_SIZE - 1);
+  const unsigned int warpId = idx >> WARP_SIZE2;
+
+  int val = inclusive_scan_warp<WARP_SIZE2>(v_in, idx);
+  if (31 == laneId) shdata[warpId] = val;
+  __syncthreads();
+
+  if (0 == warpId) shdata[idx] = inclusive_scan_warp<BLOCKDIM2 - WARP_SIZE2>(shdata[idx], idx);
+  __syncthreads();
+
+  if (warpId > 0) val += shdata[warpId - 1];
+  __syncthreads();
+  if ((1 << BLOCKDIM2) - 1 == idx) shdata[0] = val;
+  __syncthreads();
+
+  return make_int2(val, shdata[0]);
+}
+
+template<const int BLOCKDIM2>
+__device__ int inclusive_scan_array(volatile int *ptr, volatile int* shdata, const int N, const unsigned int idx) 
+{
+  int y = 0;
+
+  for (int p = 0; p < N; p += blockDim.x) 
+  {
+    __syncthreads();
+    ptr[p + idx] = inclusive_scan_block<BLOCKDIM2>(shdata, ptr[p+idx], idx);
+    ptr[p + idx] += y;
+
+    __syncthreads();
+
+    y = ptr[blockDim.x - 1];
+  }
+
+  return y;
+}
 
 
 #endif



#define TEXTURES
// #define OLDPREFIX


template<int DIM2, int SHIFT>

        sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory

        // step 1
#if 0
        int v0 = inclusive_scan_block<DIM2>(prefix_shmem, (int)flag, tid);       // inclusive scan on flags to construct array
        prefix0[tid] = v0;
        __syncthreads();

        if (flag) prefix1[prefix0[tid] - 1] = tid;                             //with tidś whose leaves have to be opened
        __syncthreads();                                                      //thread barrier, make sure all warps completed the job
#else
        offset = inclusive_scan_block<DIM2>(prefix_shmem, (int)flag, tid);       // inclusive scan on flags to construct array
        if (flag) prefix1[offset - 1] = tid;                             //with tidś whose leaves have to be opened
#endif

        // step 2
#if 0
        int v0 = inclusive_scan_block<DIM2>(prefix_shmem, nbody, tid);        // inclusive scan to compute memory offset for each body
        prefix0[tid] = v0;
        __syncthreads();
        int n_bodies = prefix0[blockDim.x - 1];                            //Total number of bides extract from the leaves
        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs

        direct [tid]  = prefix0[tid];                                       //Store a copy of inclusive scan in direct
        prefix0[tid] -= nbody;                                              //convert inclusive int oexclusive scan
        prefix0[tid] += 1;                                                  //add unity, since later prefix0[tid] == 0 used to check barrier
#else
        offset2 = inclusive_scan_blockS<DIM2>(prefix_shmem, nbody, tid);        // inclusive scan to compute memory offset for each body
        int offset1  = offset2.x;
        int n_bodies = offset2.y;

        direct[tid] = offset1;
        offset1    -= nbody;
        offset1    += 1;
        prefix0[tid] = offset1;
        __syncthreads();
#endif

        int nl_pre = 0;                                                     //Number of leaves that have already been processed

