2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
35 #ifndef _nbnxn_kernel_simd_utils_x86_128d_h_
36 #define _nbnxn_kernel_simd_utils_x86_128d_h_
38 #include "gromacs/legacyheaders/types/simple.h"
40 /* This files contains all functions/macros for the SIMD kernels
41 * which have explicit dependencies on the j-cluster size and/or SIMD-width.
42 * The functionality which depends on the j-cluster size is:
45 * energy group pair energy storage
48 #define gmx_mm_extract_epi32(x, imm) _mm_cvtsi128_si32(_mm_srli_si128((x), 4 * (imm)))
50 typedef gmx_simd_int32_t gmx_exclfilter;
51 /* This is set to a constant for now, since the code does not adapt automatically just
52 * because we set the SIMD widths to other values.
54 static const int filter_stride = 2;
56 /* Transpose 2 double precision registers */
57 static gmx_inline void
58 gmx_mm_transpose2_op_pd(__m128d in0, __m128d in1,
59 __m128d *out0, __m128d *out1)
61 *out0 = _mm_unpacklo_pd(in0, in1);
62 *out1 = _mm_unpackhi_pd(in0, in1);
65 /* Sum the elements within each input register and store the sums in out */
66 static gmx_inline __m128d
67 gmx_mm_transpose_sum2_pr(__m128d in0, __m128d in1)
71 gmx_mm_transpose2_op_pd(in0, in1, &tr0, &tr1);
73 return _mm_add_pd(tr0, tr1);
77 gmx_mm128_invsqrt_ps_single(__m128 x)
79 const __m128 half = _mm_set_ps(0.5, 0.5, 0.5, 0.5);
80 const __m128 three = _mm_set_ps(3.0, 3.0, 3.0, 3.0);
82 __m128 lu = _mm_rsqrt_ps(x);
84 return _mm_mul_ps(half, _mm_mul_ps(_mm_sub_ps(three, _mm_mul_ps(_mm_mul_ps(lu, lu), x)), lu));
87 /* Do 2 double precision invsqrt operations.
88 * Doing the SIMD rsqrt and the first Newton Raphson iteration
89 * in single precision gives full double precision accuracy.
91 static gmx_inline void
92 gmx_mm_invsqrt2_pd(__m128d in0, __m128d in1,
93 __m128d *out0, __m128d *out1)
95 const __m128d half = _mm_set1_pd(0.5);
96 const __m128d three = _mm_set1_pd(3.0);
100 s = _mm_movelh_ps(_mm_cvtpd_ps(in0), _mm_cvtpd_ps(in1));
101 ir = gmx_mm128_invsqrt_ps_single(s);
102 lu0 = _mm_cvtps_pd(ir);
103 lu1 = _mm_cvtps_pd(_mm_movehl_ps(ir, ir));
104 *out0 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu0, lu0), in0)), lu0));
105 *out1 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu1, lu1), in1)), lu1));
108 static gmx_inline void
109 load_lj_pair_params(const real *nbfp, const int *type, int aj,
110 __m128d *c6_S, __m128d *c12_S)
112 __m128d clj_S[UNROLLJ];
115 for (p = 0; p < UNROLLJ; p++)
117 clj_S[p] = _mm_load_pd(nbfp+type[aj+p]*nbfp_stride);
119 gmx_mm_transpose2_op_pd(clj_S[0], clj_S[1], c6_S, c12_S);
122 /* The load_table functions below are performance critical.
123 * The routines issue UNROLLI*UNROLLJ _mm_load_ps calls.
124 * As these all have latencies, scheduling is crucial.
125 * The Intel compilers and CPUs seem to do a good job at this.
126 * But AMD CPUs perform significantly worse with gcc than with icc.
127 * Performance is improved a bit by using the extract function UNROLLJ times,
128 * instead of doing an _mm_store_si128 for every i-particle.
129 * This is only faster when we use FDV0 formatted tables, where we also need
130 * to multiple the index by 4, which can be done by a SIMD bit shift.
131 * With single precision AVX, 8 extracts are much slower than 1 store.
132 * Because of this, the load_table_f function always takes the ti
133 * parameter, which should contain a buffer that is aligned with
134 * prepare_table_load_buffer(), but it is only used with full-width
137 static gmx_inline void
138 load_table_f(const real *tab_coul_F, gmx_simd_int32_t ti_S, int gmx_unused *ti,
139 __m128d *ctab0_S, __m128d *ctab1_S)
144 /* Without SSE4.1 the extract macro needs an immediate: unroll */
145 idx[0] = gmx_mm_extract_epi32(ti_S, 0);
146 ctab_S[0] = _mm_loadu_pd(tab_coul_F+idx[0]);
147 idx[1] = gmx_mm_extract_epi32(ti_S, 1);
148 ctab_S[1] = _mm_loadu_pd(tab_coul_F+idx[1]);
150 /* Shuffle the force table entries to a convenient order */
151 gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], ctab0_S, ctab1_S);
152 /* The second force table entry should contain the difference */
153 *ctab1_S = _mm_sub_pd(*ctab1_S, *ctab0_S);
156 static gmx_inline void
157 load_table_f_v(const real *tab_coul_F, const real *tab_coul_V,
158 gmx_simd_int32_t ti_S, int gmx_unused *ti,
159 __m128d *ctab0_S, __m128d *ctab1_S, __m128d *ctabv_S)
164 /* Without SSE4.1 the extract macro needs an immediate: unroll */
165 idx[0] = gmx_mm_extract_epi32(ti_S, 0);
166 ctab_S[0] = _mm_loadu_pd(tab_coul_F+idx[0]);
167 idx[1] = gmx_mm_extract_epi32(ti_S, 1);
168 ctab_S[1] = _mm_loadu_pd(tab_coul_F+idx[1]);
170 /* Shuffle the force table entries to a convenient order */
171 gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], ctab0_S, ctab1_S);
172 /* The second force table entry should contain the difference */
173 *ctab1_S = _mm_sub_pd(*ctab1_S, *ctab0_S);
175 ctab_S[2] = _mm_loadu_pd(tab_coul_V+idx[0]);
176 ctab_S[3] = _mm_loadu_pd(tab_coul_V+idx[1]);
178 /* Shuffle the energy table entries to a single register */
179 *ctabv_S = _mm_shuffle_pd(ctab_S[2], ctab_S[3], _MM_SHUFFLE2(0, 0));
182 static gmx_inline gmx_exclfilter
183 gmx_load1_exclfilter(int e)
185 return _mm_set1_epi32(e);
188 static gmx_inline gmx_exclfilter
189 gmx_load_exclusion_filter(const unsigned *i)
191 /* For now this has to be an explicit-float load since we use stride==2 */
192 return gmx_simd_load_fi(i);
195 static gmx_inline gmx_simd_bool_t
196 gmx_checkbitmask_pb(gmx_exclfilter m0, gmx_exclfilter m1)
198 return _mm_castsi128_pd(_mm_cmpeq_epi32(_mm_andnot_si128(m0, m1), _mm_setzero_si128()));
201 #endif /* _nbnxn_kernel_simd_utils_x86_s128d_h_ */