     }
 
     bid   = shmemv[0];
+
     if (bid >= n_active_groups) return;
 
-    int grpOffset = 0;
-
-    /*********** set necessary thread constants **********/
-#ifdef DO_BLOCK_TIMESTEP
-    real4 curGroupSize    = groupSizeInfo[active_groups[bid + grpOffset]];
+    int *lmem = &MEM_BUF[(CELL_LIST_MEM_PER_WARP<<nWarps2)*blockIdx.x + CELL_LIST_MEM_PER_WARP*warpId];
+    const bool success = treewalk<0,blockDim2,ACCUMULATE>(
+        bid, 
+        eps2,
+        node_begend,
+        active_groups,
+        group_body_pos,
+        groupSizeInfo,
+        groupCenterInfo,
+        shmem,
+        lmem,
+        acc_out,
+        interactions,
+        ngb_out,
+        active_inout);
+
+#if 0
+    if (bid % 10 == 0)
+      success = false;
+#endif
+
+#ifdef SHMODE
+    if (!success)
+      if (laneId == 0)
+        failedList[atomicAdd(&failed,1)] = bid;
+
+    if (failed + nWarps >= MAXFAILED)
+    {
+      __syncthreads();
+      if (warpId == 0)
+      {
+        int *lmem1 = &MEM_BUF[(CELL_LIST_MEM_PER_WARP<<nWarps2)*blockIdx.x];
+        const int n = failed;
+        failed = 0;
+        for (int it = 0; it < n; it++)
+        {
+          const bool success = treewalk<nWarp2,blockDim2,ACCUMULATE>(
+              failedList[it], 
+              eps2,
+              node_begend,
+              active_groups,
+              group_body_pos,
+              groupSizeInfo,
+              groupCenterInfo,
+              shmem,
+              lmem1,
+              acc_out,
+              interactions,
+              ngb_out,
+              active_inout);
+          assert(success);
+        }
+      }    
+      __syncthreads();
+    }
+
 #else
-    real4 curGroupSize    = groupSizeInfo[bid + grpOffset];
-#endif
-    const int   groupData       = __float_as_int(curGroupSize.w);
-    const uint body_addr        =   groupData & CRITMASK;
-    const uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;
-
-#ifdef DO_BLOCK_TIMESTEP
-    real4 group_pos       = groupCenterInfo[active_groups[bid + grpOffset]];
-#else
-    real4 group_pos       = groupCenterInfo[bid + grpOffset];
-#endif
-
-    uint body_i[2];
-    int ni = nb_i <= WARP_SIZE ? 1 : 2;
-    body_i[0] = body_addr + laneId%nb_i;
-    body_i[1] = body_addr + WARP_SIZE + laneId%(nb_i - WARP_SIZE);
-
-    float4 pos_i[2];
-    float4 acc_i[2];
-
-    pos_i[0] = group_body_pos[body_i[0]];
-    if(ni > 1) //Only read if we actually have ni == 2
-      pos_i[1] = group_body_pos[body_i[1]];
-
-    acc_i[0] = acc_i[1] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
-
-
-
-    int apprCount = 0;
-    int direCount = 0;
-
-    const float group_eps  = 0;
-    int ngb_i;
-#if 1
-    if (ni == 1)
-      approximate_gravity<0, blockDim2, 1>(
-          pos_i, group_pos,
-          eps2, node_begend,
-          multipole_data, body_pos,
-          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
-          group_eps, 
-          acc_i);
-    else
-      approximate_gravity<0, blockDim2, 2>(
-          pos_i, group_pos,
-          eps2, node_begend,
-          multipole_data, body_pos,
-          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
-          group_eps, 
-          acc_i);
-#endif
-
-    assert(apprCount >= 0);
-#if 0 /* this increase lmem spill count */
-    if(apprCount < 0)
+
+    //Try to get access to the big stack, only one block per time is allowed
+    if (!success)
     {
-
-      //Try to get access to the big stack, only one block per time is allowed
       if(laneId == 0)
       {
         int res = atomicExch(&active_inout[n_bodies+1], 1); //If the old value (res) is 0 we can go otherwise sleep

        {
          //Sleep
          for(int i=0; i < (1024); i++)
          {
            waitCounter += 1;
          }
          //Test again
          shmem[0] = waitCounter;
          res = atomicExch(&active_inout[n_bodies+1], 1); 
        }
      }

      lmem = &MEM_BUF[gridDim.x*(LMEM_STACK_SIZE*blockDim.x + LMEM_EXTRA_SIZE)];    //Use the extra large buffer
      apprCount = direCount = 0;
      acc_i[0] = acc_i[1] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
      if (ni == 1)    
        approximate_gravity<8, blockDim2, 1>(
            pos_i, group_pos,
            eps2, node_begend,
            multipole_data, body_pos,
            shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
            group_eps, 
            acc_i);
      else
        approximate_gravity<8, blockDim2, 2>(
            pos_i, group_pos,
            eps2, node_begend,
            multipole_data, body_pos,
            shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
            group_eps, 
            acc_i);

      lmem = &MEM_BUF[blockIdx.x*(LMEM_STACK_SIZE*blockDim.x + LMEM_EXTRA_SIZE)];

      if(laneId == 0)
      {
        atomicExch(&active_inout[n_bodies+1], 0); //Release the lock
      }
    }//end if apprCount < 0
#endif

    if (laneId < nb_i) 
    {
      const int addr = body_i[0];
      acc_out     [addr] = acc_i[0];
      //       ngb_out     [addr] = ngb_i;
      ngb_out     [addr] = addr; //JB Fixed this for demo 
      active_inout[addr] = 1;
      interactions[addr].x = apprCount / ni;
      interactions[addr].y = direCount / ni ;
      if (ni == 2)
      {
        const int addr = body_i[1];
        acc_out     [addr] = acc_i[1];
        //         ngb_out     [addr] = ngb_i;
        ngb_out     [addr] = addr; //JB Fixed this for demo 
        active_inout[addr] = 1;     
        interactions[addr].x = apprCount / ni; 
        interactions[addr].y = direCount / ni; 
      }
    }
  }     //end while
}
#else
  extern "C"
__launch_bounds__(NTHREAD,1024/NTHREAD)
  __global__ void

  }  /* level completed */
#endif

  if (approxCounter > 0)
  {
    approxAcc<NI,false>(acc_i, pos_i, laneIdx < approxCounter ? approxCellIdx : -1, eps2);
    if (INTCOUNT)
      interactionCounters.x += approxCounter * NI;
    approxCounter = 0;
  }

  if (directCounter > 0)
  {
    directAcc<NI,false>(acc_i, pos_i, laneIdx < directCounter ? directPtclIdx : -1, eps2);
    if (INTCOUNT)
      interactionCounters.y += directCounter * NI;
    directCounter = 0;
  }

  return interactionCounters;
}

template<int SHIFT2, int BLOCKDIM2, bool ACCUMULATE>
static __device__ 
bool treewalk(
    const int bid,
    const float eps2,
    const uint2 node_begend,
    const int    *active_groups,
    const real4  *group_body_pos,        
    const float4  *groupSizeInfo,
    const float4  *groupCenterInfo,
    int *shmem,
    int *lmem,
    float4 *acc_out,
    int2   *interactions,
    int    *ngb_out,
    int    *active_inout)
{

  /*********** set necessary thread constants **********/
#ifdef DO_BLOCK_TIMESTEP
  real4 curGroupSize    = groupSizeInfo[active_groups[bid]];
#else
  real4 curGroupSize    = groupSizeInfo[bid];
#endif
  const int   groupData       = __float_as_int(curGroupSize.w);
  const uint body_addr        =   groupData & CRITMASK;
  const uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;

#ifdef DO_BLOCK_TIMESTEP
  real4 group_pos       = groupCenterInfo[active_groups[bid]];
#else
  real4 group_pos       = groupCenterInfo[bid];
#endif


  uint body_i[2];
  const int ni = nb_i <= WARP_SIZE ? 1 : 2;
  body_i[0] = body_addr + laneId%nb_i;
  body_i[1] = body_addr + WARP_SIZE + laneId%(nb_i - WARP_SIZE);

  float4 pos_i[2];
  float4 acc_i[2];

  pos_i[0] = group_body_pos[body_i[0]];
  if(ni > 1) //Only read if we actually have ni == 2
    pos_i[1] = group_body_pos[body_i[1]];

  acc_i[0] = acc_i[1] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);


#if 0
  const bool INTCOUNT = false;
#else
  const bool INTCOUNT = true;
#endif
  uint2 counters = {0};
  {
    if (ni == 1)
      counters = approximate_gravity<SHIFT2, BLOCKDIM2, 1,INTCOUNT>(
          acc_i,
          pos_i,
          group_pos,
          eps2, 
          node_begend,
          shmem, 
          lmem, 
          curGroupSize);
    else
      counters = approximate_gravity<SHIFT2, BLOCKDIM2, 2,INTCOUNT>(
          acc_i,
          pos_i,
          group_pos,
          eps2, 
          node_begend,
          shmem, 
          lmem, 
          curGroupSize);
  }
  if(counters.x == 0xFFFFFFFF && counters.y == 0xFFFFFFFF)
    return false;

  if (laneId < nb_i) 
  {
    const int addr = body_i[0];
    if (ACCUMULATE)
    {
      acc_out     [addr].x += acc_i[0].x;
      acc_out     [addr].y += acc_i[0].y;
      acc_out     [addr].z += acc_i[0].z;
      acc_out     [addr].w += acc_i[0].w;
    }
    else
      acc_out     [addr] = acc_i[0];
    //       ngb_out     [addr] = ngb_i;
    ngb_out     [addr] = addr; //JB Fixed this for demo 
    active_inout[addr] = 1;
    interactions[addr].x = counters.x / ni;
    interactions[addr].y = counters.y / ni ;
    if (ni == 2)
    {
      const int addr = body_i[1];
      if (ACCUMULATE)
      {
        acc_out     [addr].x += acc_i[1].x;
        acc_out     [addr].y += acc_i[1].y;
        acc_out     [addr].z += acc_i[1].z;
        acc_out     [addr].w += acc_i[1].w;
      }
      else
      {
        acc_out     [addr] = acc_i[1];
      }

      //         ngb_out     [addr] = ngb_i;
      ngb_out     [addr] = addr; //JB Fixed this for demo 
      active_inout[addr] = 1;     
      interactions[addr].x = counters.x / ni; 
      interactions[addr].y = counters.y / ni; 
    }
  }

  return true;
}

template<bool ACCUMULATE, int BLOCKDIM2>
static __device__
void approximate_gravity_main(
      const int n_active_groups,
      int    n_bodies,
      float eps2,
