             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
+#if 1    /*  a bit faster */
+          const float4  Q0  = tex1Dfetch(texMultipole, address + 1);
+          const float4  Q1  = tex1Dfetch(texMultipole, address + 2);
+          for (int i = 0; i < WARP_SIZE; i++)
+          {
+            const float4 jQ0 = make_float4(__shfl(Q0.x, i), __shfl(Q0.y, i), __shfl(Q0.z, i), 0.0f);
+            const float4 jQ1 = make_float4(__shfl(Q1.x, i), __shfl(Q1.y, i), __shfl(Q1.z, i), 0.0f);
+            for (int k = 0; k < NI; k++)
+              acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], jQ0, jQ1, eps2);
+          }
+#else
           for (int i = 0; i < WARP_SIZE; i++)
           {
             const int address = approxM[n_approx + i] * 3;


  float Qii = oct_q11 + oct_q22 + oct_q33;
  float QijRiRj =
    (oct_q11*dr.x*dr.x + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z) +
    2.0f*(oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z);

  acc.w        -= D0 + 0.5f*D1*Qii + 0.5f*D2*QijRiRj;

            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
#endif /* _QUADRUPOLE_ */
          apprCount += WARP_SIZE*NI;
        }
