From a5a722ac2cb78674f7765e3739563dd957fc8cec Mon Sep 17 00:00:00 2001
From: James Willis <james.s.willis@durham.ac.uk>
Date: Tue, 29 Nov 2016 14:35:05 +0000
Subject: [PATCH] Added vectorised version of non-sym density function that
interleaves 2 sets of vector operations.
---
src/hydro/Gadget2/hydro_iact.h | 124 +++++++++++++++++++++++++++++++++
1 file changed, 124 insertions(+)
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 08fb2b37db..12cc778d74 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -382,6 +382,130 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
#endif
}
+/**
+ * @brief Density interaction computed using 2 interleaved vectors (non-symmetric vectorized version).
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_2_vec_density(float *R2, float *Dx, float *Dy, float *Dz, vector hi_inv,
+ vector vix, vector viy, vector viz, float *Vjx, float *Vjy, float *Vjz, float *Mj, vector *rhoSum, vector *rho_dhSum, vector *wcountSum, vector *wcount_dhSum, vector *div_vSum, vector *curlvxSum,vector *curlvySum, vector *curlvzSum, vector mask, vector mask2, int knlMask, int knlMask2) {
+
+ vector r, ri, r2, xi, wi, wi_dx;
+ vector mj;
+ vector dx, dy, dz, dvx, dvy, dvz;
+ vector vjx, vjy, vjz;
+ vector dvdr;
+ vector curlvrx, curlvry, curlvrz;
+ vector r_2, ri2, r2_2, xi2, wi2, wi_dx2;
+ vector mj2;
+ vector dx2, dy2, dz2, dvx2, dvy2, dvz2;
+ vector vjx2, vjy2, vjz2;
+ vector dvdr2;
+ vector curlvrx2, curlvry2, curlvrz2;
+
+ /* Get the masses. */
+ mj.v = vec_load(Mj);
+ mj2.v = vec_load(&Mj[VEC_SIZE]);
+ vjx.v = vec_load(Vjx);
+ vjx2.v = vec_load(&Vjx[VEC_SIZE]);
+ vjy.v = vec_load(Vjy);
+ vjy2.v = vec_load(&Vjy[VEC_SIZE]);
+ vjz.v = vec_load(Vjz);
+ vjz2.v = vec_load(&Vjz[VEC_SIZE]);
+ dx.v = vec_load(Dx);
+ dx2.v = vec_load(&Dx[VEC_SIZE]);
+ dy.v = vec_load(Dy);
+ dy2.v = vec_load(&Dy[VEC_SIZE]);
+ dz.v = vec_load(Dz);
+ dz2.v = vec_load(&Dz[VEC_SIZE]);
+
+ /* Get the radius and inverse radius. */
+ r2.v = vec_load(R2);
+ r2_2.v = vec_load(&R2[VEC_SIZE]);
+ VEC_RECIPROCAL_SQRT(r2.v, ri.v);
+ VEC_RECIPROCAL_SQRT(r2_2.v, ri2.v);
+ r.v = vec_mul(r2.v, ri.v);
+ r_2.v = vec_mul(r2_2.v, ri2.v);
+
+ xi.v = vec_mul(r.v, hi_inv.v);
+ xi2.v = vec_mul(r_2.v, hi_inv.v);
+
+ //kernel_deval_2_vec(&xi, &wi, &wi_dx,&xi2, &wi2, &wi_dx2);
+ kernel_deval_vec(&xi, &wi, &wi_dx);
+ kernel_deval_vec(&xi2, &wi2, &wi_dx2);
+
+ /* Compute dv. */
+ dvx.v = vec_sub(vix.v, vjx.v);
+ dvx2.v = vec_sub(vix.v, vjx2.v);
+ dvy.v = vec_sub(viy.v, vjy.v);
+ dvy2.v = vec_sub(viy.v, vjy2.v);
+ dvz.v = vec_sub(viz.v, vjz.v);
+ dvz2.v = vec_sub(viz.v, vjz2.v);
+
+ /* Compute dv dot r */
+ dvdr.v = vec_fma(dvx.v, dx.v, vec_fma(dvy.v, dy.v, vec_mul(dvz.v, dz.v)));
+ dvdr2.v = vec_fma(dvx2.v, dx2.v, vec_fma(dvy2.v, dy2.v, vec_mul(dvz2.v, dz2.v)));
+ dvdr.v = vec_mul(dvdr.v, ri.v);
+ dvdr2.v = vec_mul(dvdr2.v, ri2.v);
+
+ /* Compute dv cross r */
+ curlvrx.v = vec_fma(dvy.v, dz.v, vec_mul(vec_set1(-1.0f),vec_mul(dvz.v, dy.v)));
+ curlvrx2.v = vec_fma(dvy2.v, dz2.v, vec_mul(vec_set1(-1.0f),vec_mul(dvz2.v, dy2.v)));
+ curlvry.v = vec_fma(dvz.v, dx.v, vec_mul(vec_set1(-1.0f), vec_mul(dvx.v, dz.v)));
+ curlvry2.v = vec_fma(dvz2.v, dx2.v, vec_mul(vec_set1(-1.0f), vec_mul(dvx2.v, dz2.v)));
+ curlvrz.v = vec_fma(dvx.v, dy.v, vec_mul(vec_set1(-1.0f), vec_mul(dvy.v, dx.v)));
+ curlvrz2.v = vec_fma(dvx2.v, dy2.v, vec_mul(vec_set1(-1.0f), vec_mul(dvy2.v, dx2.v)));
+ curlvrx.v = vec_mul(curlvrx.v,ri.v);
+ curlvrx2.v = vec_mul(curlvrx2.v,ri2.v);
+ curlvry.v = vec_mul(curlvry.v,ri.v);
+ curlvry2.v = vec_mul(curlvry2.v,ri2.v);
+ curlvrz.v = vec_mul(curlvrz.v,ri.v);
+ curlvrz2.v = vec_mul(curlvrz2.v,ri2.v);
+
+ /* Mask updates to intermediate vector sums for particle pi. */
+#ifdef HAVE_AVX512_F
+ rhoSum->v = _mm512_mask_add_ps(rhoSum->v, knlMask, vec_mul(mj.v, wi.v), rhoSum->v);
+ rhoSum->v = _mm512_mask_add_ps(rhoSum->v, knlMask2, vec_mul(mj2.v, wi2.v), rhoSum->v);
+
+ rho_dhSum->v = _mm512_mask_sub_ps(rho_dhSum->v, knlMask, rho_dhSum->v, vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(xi.v, wi_dx.v))));
+ rho_dhSum->v = _mm512_mask_sub_ps(rho_dhSum->v, knlMask2, rho_dhSum->v, vec_mul(mj2.v, vec_fma(vec_set1(hydro_dimension), wi2.v, vec_mul(xi2.v, wi_dx2.v))));
+
+ wcountSum->v = _mm512_mask_add_ps(wcountSum->v, knlMask, wi.v, wcountSum->v);
+ wcountSum->v = _mm512_mask_add_ps(wcountSum->v, knlMask2, wi2.v, wcountSum->v);
+
+ wcount_dhSum->v = _mm512_mask_sub_ps(wcount_dhSum->v, knlMask, wcount_dhSum->v, vec_mul(xi.v, wi_dx.v));
+ wcount_dhSum->v = _mm512_mask_sub_ps(wcount_dhSum->v, knlMask2, wcount_dhSum->v, vec_mul(xi2.v, wi_dx2.v));
+
+ div_vSum->v = _mm512_mask_sub_ps(div_vSum->v, knlMask, div_vSum->v, vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)));
+ div_vSum->v = _mm512_mask_sub_ps(div_vSum->v, knlMask2, div_vSum->v, vec_mul(mj2.v, vec_mul(dvdr2.v, wi_dx2.v)));
+
+ curlvxSum->v = _mm512_mask_add_ps(curlvxSum->v, knlMask, vec_mul(mj.v, vec_mul(curlvrx.v, wi_dx.v)), curlvxSum->v);
+ curlvxSum->v = _mm512_mask_add_ps(curlvxSum->v, knlMask2, vec_mul(mj2.v, vec_mul(curlvrx2.v, wi_dx2.v)), curlvxSum->v);
+
+ curlvySum->v = _mm512_mask_add_ps(curlvySum->v, knlMask, vec_mul(mj.v, vec_mul(curlvry.v, wi_dx.v)), curlvySum->v);
+ curlvySum->v = _mm512_mask_add_ps(curlvySum->v, knlMask2, vec_mul(mj2.v, vec_mul(curlvry2.v, wi_dx2.v)), curlvySum->v);
+
+ curlvzSum->v = _mm512_mask_add_ps(curlvzSum->v, knlMask, vec_mul(mj.v, vec_mul(curlvrz.v, wi_dx.v)), curlvzSum->v);
+ curlvzSum->v = _mm512_mask_add_ps(curlvzSum->v, knlMask2, vec_mul(mj2.v, vec_mul(curlvrz2.v, wi_dx2.v)), curlvzSum->v);
+#else
+ rhoSum->v += vec_and(vec_mul(mj.v, wi.v),mask.v);
+ rhoSum->v += vec_and(vec_mul(mj2.v, wi2.v),mask2.v);
+ rho_dhSum->v -= vec_and(vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(xi.v, wi_dx.v))),mask.v);
+ rho_dhSum->v -= vec_and(vec_mul(mj2.v, vec_fma(vec_set1(hydro_dimension), wi2.v, vec_mul(xi2.v, wi_dx2.v))),mask2.v);
+ wcountSum->v += vec_and(wi.v,mask.v);
+ wcountSum->v += vec_and(wi2.v,mask2.v);
+ wcount_dhSum->v -= vec_and(vec_mul(xi.v, wi_dx.v),mask.v);
+ wcount_dhSum->v -= vec_and(vec_mul(xi2.v, wi_dx2.v),mask2.v);
+ div_vSum->v -= vec_and(vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)),mask.v);
+ div_vSum->v -= vec_and(vec_mul(mj2.v, vec_mul(dvdr2.v, wi_dx2.v)),mask2.v);
+ curlvxSum->v += vec_and(vec_mul(mj.v, vec_mul(curlvrx.v, wi_dx.v)),mask.v);
+ curlvxSum->v += vec_and(vec_mul(mj2.v, vec_mul(curlvrx2.v, wi_dx2.v)),mask2.v);
+ curlvySum->v += vec_and(vec_mul(mj.v, vec_mul(curlvry.v, wi_dx.v)),mask.v);
+ curlvySum->v += vec_and(vec_mul(mj2.v, vec_mul(curlvry2.v, wi_dx2.v)),mask2.v);
+ curlvzSum->v += vec_and(vec_mul(mj.v, vec_mul(curlvrz.v, wi_dx.v)),mask.v);
+ curlvzSum->v += vec_and(vec_mul(mj2.v, vec_mul(curlvrz2.v, wi_dx2.v)),mask2.v);
+#endif
+}
+
/**
* @brief Force loop
*/
--
GitLab