   }
 
   /***
-  **** --> reduce data between threads
-  ***/
+   **** --> reduce data between threads
+   ***/
   sh_pot[tid] = acc_i.w;
   sh_acc[tid] = make_float3(acc_i.x, acc_i.y, acc_i.z);
   sh_ds2[tid] = ds2_min;



  if (ty == 0) {
  #pragma unroll
    for (int i = 1; i < DIMy; i++){
        int idx = (i << DIM2x) + tx;
        direCount  += sh_ds2[idx];
        apprCount  += sh_ngb[idx];
    }
  }
  __syncthreads();

}


  extern "C" __global__ void
__launch_bounds__(NTHREAD)
  dev_approximate_gravity(const int n_active_groups,
      int    n_bodies,
      float eps2,
      uint2 node_begend,
      int    *active_groups,
      real4  *body_pos,
      real4  *multipole_data,
      float4 *acc_out,
      int    *ngb_out,
      int    *active_inout,
      int2   *interactions,
      float4  *boxSizeInfo,
      float4  *groupSizeInfo,
      float4  *boxCenterInfo,
      float4  *groupCenterInfo,
      real4   *body_vel,
      int     *MEM_BUF) {
    //                                                    int     grpOffset){


    const int blockDim2 = NTHREAD2;
    __shared__ int shmem[15*(1 << blockDim2)];
    //    __shared__ int shmem[24*(1 << blockDim2)]; is possible on FERMI
    //    int             lmem[LMEM_STACK_SIZE];



    /*********** check if this block is linked to a leaf **********/

    int bid = gridDim.x * blockIdx.y + blockIdx.x;

    while(true)
    {

      if(threadIdx.x == 0)
      {
        bid         = atomicAdd(&active_inout[n_bodies], 1);
        shmem[0]    = bid;
      }
      __syncthreads();

      bid   = shmem[0];

      if (bid >= n_active_groups) return;


      int tid = threadIdx.y * blockDim.x + threadIdx.x;

      int grpOffset = 0;

      //   volatile int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
      //   int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x + threadIdx.x*LMEM_STACK_SIZE];
      int *lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x];


      /*********** set necessary thread constants **********/

      real4 curGroupSize    = groupSizeInfo[active_groups[bid + grpOffset]];
      int   groupData       = __float_as_int(curGroupSize.w);
      uint body_i           =   groupData & CRITMASK;
      uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;

      real4 group_pos       = groupCenterInfo[active_groups[bid + grpOffset]];

      //   if(tid == 0)
      //   printf("[%f %f %f %f ] \n [%f %f %f %f ] %d %d \n",
      //           curGroupSize.x, curGroupSize.y, curGroupSize.z, curGroupSize.w,
      //           group_pos.x, group_pos.y, group_pos.z, group_pos.w, body_i, nb_i);


      int DIM2x = 0;
      while (((nb_i - 1) >> DIM2x) > 0) DIM2x++;

      DIM2x     = max(DIM2x,4);
      int DIM2y = blockDim2 - DIM2x;

      int tx = tid & ((1 << DIM2x) - 1);
      int ty = tid >> DIM2x;

      body_i += tx%nb_i;

      //float4 pos_i = tex1Dfetch(bodies_pos_ref, body_i);   // texture read: 4 floats


      float4 pos_i = body_pos[body_i];


      int ngb_i;

      float4 acc_i = {0.0f, 0.0f, 0.0f, 0.0f};

#ifdef INDSOFT
      eps2 = body_vel[body_i].w;
      float group_eps = eps2;

      volatile float *reduc = (float*) &shmem[0];
      reduc[threadIdx.x] = eps2;

      //Find the maximum softening value for the particles in this group
      __syncthreads();
      // do reduction in shared mem
      if(blockDim.x >= 512) if (tid < 256) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 256]);} __syncthreads();
      if(blockDim.x >= 256) if (tid < 128) {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 128]);} __syncthreads();
      if(blockDim.x >= 128) if (tid < 64)  {reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 64]);} __syncthreads();
      if(blockDim.x >= 64) if (tid < 32) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 32]);}
      if(blockDim.x >= 32) if (tid < 16) { reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 16]);}

      if(tid < 8)
      {
        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 8]);
        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 4]);
        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 2]);
        reduc[threadIdx.x] = group_eps = fmaxf(group_eps, reduc[threadIdx.x + 1]);
      }
      __syncthreads();

      group_eps = reduc[0];
#else
      float group_eps  = 0;
#endif

      int apprCount = 0;
      int direCount = 0;


      acc_i = approximate_gravity<blockDim2, 0>( DIM2x, DIM2y, tid, tx, ty,
          body_i, pos_i, group_pos,
          eps2, node_begend,
          multipole_data, body_pos,
          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
          group_eps, body_vel);
      if(apprCount < 0)
      {

        //Try to get access to the big stack, only one block per time is allowed
        if(threadIdx.x == 0)
        {
          int res = atomicExch(&active_inout[n_bodies+1], 1); //If the old value (res) is 0 we can go otherwise sleep
          int waitCounter  = 0;
          while(res != 0)
          {
            //Sleep
            for(int i=0; i < (1024); i++)
            {
              waitCounter += 1;
            }
            //Test again
            shmem[0] = waitCounter;
            res = atomicExch(&active_inout[n_bodies+1], 1); 
          }
        }

        __syncthreads();

        lmem = &MEM_BUF[gridDim.x*LMEM_STACK_SIZE*blockDim.x];    //Use the extra large buffer
        apprCount = direCount = 0;
        acc_i = approximate_gravity<blockDim2, 8>( DIM2x, DIM2y, tid, tx, ty,
            body_i, pos_i, group_pos,
            eps2, node_begend,
            multipole_data, body_pos,
            shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
            group_eps, body_vel);

        lmem = &MEM_BUF[blockIdx.x*LMEM_STACK_SIZE*blockDim.x]; //Back to normal location

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

      if (tid < nb_i) {
        acc_out     [body_i] = acc_i;
        ngb_out     [body_i] = ngb_i;
        active_inout[body_i] = 1;
        interactions[body_i].x = apprCount;
        interactions[body_i].y = direCount ;
      }


    }     //end while
  }


