             sh_mass[laneId] = posj.w;
             sh_pos [laneId] = make_float3(posj.x, posj.y, posj.z);
 
-#pragma unroll 16
-            for (int j = 0; j < WARP_SIZE; j++)
+            for (int i = 0; i < WARP_SIZE; i++)
               for (int k = 0; k < NI; k++)
-                acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[offs + j], sh_pos[offs + j], eps2);
+                acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], eps2);
             direCount += WARP_SIZE*NI;
           }
 

      sh_pos [laneId] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
    }

#pragma unroll 16
    for (int j = 0; j < WARP_SIZE; j++) 
      for (int k = 0; k < NI; k++)
        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[offs + j], sh_pos[offs + j], eps2);
    direCount += WARP_SIZE*NI;
  }
}

      sh_pos [laneId] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);

    }
#ifndef _QUADRUPOLE_
    for (int i = 0; i < WARP_SIZE; i++)
      for (int k = 0; k < NI; k++)
        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i],eps2);
#else
    for (int i = 0; i < WARP_SIZE; i++)
    {
      float4 Q0, Q1;
      Q0 = Q1 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
      if (i < n_approx)
      {
        const int address = approxM[i] * 3;
        Q0 = tex1Dfetch(texMultipole, address + 1);
        Q1 = tex1Dfetch(texMultipole, address + 2);
      }
      for (int k = 0; k < NI; k++)
        acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], Q0, Q1, eps2);
    }
#endif
    apprCount += WARP_SIZE*NI;

  } //if n_approx > 0
