2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2012, The GROMACS Development Team
6 * Copyright (c) 2012,2013, by the GROMACS development team, led by
7 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
8 * others, as listed in the AUTHORS file in the top-level source
9 * directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 #ifndef _nbnxn_kernel_simd_utils_x86_256d_h_
38 #define _nbnxn_kernel_simd_utils_x86_256d_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_exclfilter gmx_mm_pr
49 static const int filter_stride = 2;
51 /* Transpose 2 double precision registers */
52 static gmx_inline void
53 gmx_mm_transpose2_op_pd(__m128d in0, __m128d in1,
54 __m128d *out0, __m128d *out1)
56 *out0 = _mm_unpacklo_pd(in0, in1);
57 *out1 = _mm_unpackhi_pd(in0, in1);
60 /* Sum the elements within each input register and store the sums in out */
61 static gmx_inline __m256d
62 gmx_mm_transpose_sum4_pr(__m256d in0, __m256d in1,
63 __m256d in2, __m256d in3)
65 in0 = _mm256_hadd_pd(in0, in1);
66 in2 = _mm256_hadd_pd(in2, in3);
68 return _mm256_add_pd(_mm256_permute2f128_pd(in0, in2, 0x20), _mm256_permute2f128_pd(in0, in2, 0x31));
71 static gmx_inline __m256
72 gmx_mm256_invsqrt_ps_single(__m256 x)
74 const __m256 half = _mm256_set_ps(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5);
75 const __m256 three = _mm256_set_ps(3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0);
77 __m256 lu = _mm256_rsqrt_ps(x);
79 return _mm256_mul_ps(half, _mm256_mul_ps(_mm256_sub_ps(three, _mm256_mul_ps(_mm256_mul_ps(lu, lu), x)), lu));
82 /* Put two 128-bit 4-float registers into one 256-bit 8-float register */
83 static gmx_inline __m256
84 gmx_2_m128_to_m256(__m128 in0, __m128 in1)
86 return _mm256_insertf128_ps(_mm256_castps128_ps256(in0), in1, 1);
89 /* Put two 128-bit 2-double registers into one 256-bit 4-double register */
90 static gmx_inline __m256d
91 gmx_2_m128d_to_m256d(__m128d in0, __m128d in1)
93 return _mm256_insertf128_pd(_mm256_castpd128_pd256(in0), in1, 1);
96 /* Do 2 double precision invsqrt operations.
97 * Doing the SIMD rsqrt and the first Newton Raphson iteration
98 * in single precision gives full double precision accuracy.
100 static gmx_inline void
101 gmx_mm_invsqrt2_pd(__m256d in0, __m256d in1,
102 __m256d *out0, __m256d *out1)
104 const __m256d half = _mm256_set1_pd(0.5);
105 const __m256d three = _mm256_set1_pd(3.0);
109 s = gmx_2_m128_to_m256(_mm256_cvtpd_ps(in0), _mm256_cvtpd_ps(in1));
110 ir = gmx_mm256_invsqrt_ps_single(s);
111 lu0 = _mm256_cvtps_pd(_mm256_castps256_ps128(ir));
112 lu1 = _mm256_cvtps_pd(_mm256_extractf128_ps(ir, 1));
113 *out0 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu0, lu0), in0)), lu0));
114 *out1 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu1, lu1), in1)), lu1));
117 static gmx_inline void
118 load_lj_pair_params(const real *nbfp, const int *type, int aj,
119 __m256d *c6_S, __m256d *c12_S)
121 __m128d clj_S[UNROLLJ], c6t_S[2], c12t_S[2];
124 for (p = 0; p < UNROLLJ; p++)
126 clj_S[p] = _mm_load_pd(nbfp+type[aj+p]*nbfp_stride);
128 gmx_mm_transpose2_op_pd(clj_S[0], clj_S[1], &c6t_S[0], &c12t_S[0]);
129 gmx_mm_transpose2_op_pd(clj_S[2], clj_S[3], &c6t_S[1], &c12t_S[1]);
130 *c6_S = gmx_2_m128d_to_m256d(c6t_S[0], c6t_S[1]);
131 *c12_S = gmx_2_m128d_to_m256d(c12t_S[0], c12t_S[1]);
134 /* The load_table functions below are performance critical. They
135 * always take the ti parameter, which should contain a buffer that
136 * is aligned with prepare_table_load_buffer(), but it is only used
137 * with full-width AVX_256. */
139 static gmx_inline void
140 load_table_f(const real *tab_coul_F, __m128i ti_S, int *ti,
141 __m256d *ctab0_S, __m256d *ctab1_S)
143 __m128d ctab_S[4], tr_S[4];
146 _mm_store_si128((__m128i *)ti, ti_S);
147 for (j = 0; j < 4; j++)
149 ctab_S[j] = _mm_loadu_pd(tab_coul_F+ti[j]);
151 /* Shuffle the force table entries to a convenient order */
152 gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], &tr_S[0], &tr_S[1]);
153 gmx_mm_transpose2_op_pd(ctab_S[2], ctab_S[3], &tr_S[2], &tr_S[3]);
154 *ctab0_S = gmx_2_m128d_to_m256d(tr_S[0], tr_S[2]);
155 *ctab1_S = gmx_2_m128d_to_m256d(tr_S[1], tr_S[3]);
156 /* The second force table entry should contain the difference */
157 *ctab1_S = _mm256_sub_pd(*ctab1_S, *ctab0_S);
160 static gmx_inline void
161 load_table_f_v(const real *tab_coul_F, const real *tab_coul_V,
162 __m128i ti_S, int *ti,
163 __m256d *ctab0_S, __m256d *ctab1_S, __m256d *ctabv_S)
165 __m128d ctab_S[8], tr_S[4];
168 _mm_store_si128((__m128i *)ti, ti_S);
169 for (j = 0; j < 4; j++)
171 ctab_S[j] = _mm_loadu_pd(tab_coul_F+ti[j]);
173 /* Shuffle the force table entries to a convenient order */
174 gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], &tr_S[0], &tr_S[1]);
175 gmx_mm_transpose2_op_pd(ctab_S[2], ctab_S[3], &tr_S[2], &tr_S[3]);
176 *ctab0_S = gmx_2_m128d_to_m256d(tr_S[0], tr_S[2]);
177 *ctab1_S = gmx_2_m128d_to_m256d(tr_S[1], tr_S[3]);
178 /* The second force table entry should contain the difference */
179 *ctab1_S = _mm256_sub_pd(*ctab1_S, *ctab0_S);
181 for (j = 0; j < 4; j++)
183 ctab_S[4+j] = _mm_loadu_pd(tab_coul_V+ti[j]);
185 /* Shuffle the energy table entries to a single register */
186 *ctabv_S = gmx_2_m128d_to_m256d(_mm_shuffle_pd(ctab_S[4], ctab_S[5], _MM_SHUFFLE2(0, 0)), _mm_shuffle_pd(ctab_S[6], ctab_S[7], _MM_SHUFFLE2(0, 0)));
189 static gmx_inline gmx_exclfilter
190 gmx_load1_exclfilter(int e)
192 return _mm256_castsi256_pd(_mm256_set1_epi32(e));
195 static gmx_inline gmx_exclfilter
196 gmx_load_exclusion_filter(const unsigned *i)
198 return gmx_load_pr((real *) (i));
201 static gmx_inline gmx_mm_pb
202 gmx_checkbitmask_pb(gmx_exclfilter m0, gmx_exclfilter m1)
204 /* With <= 16 bits used the cast and conversion should not be
205 * required, since only mantissa bits are set and that would give
206 * a non-zero float, but with the Intel compiler this does not
207 * work correctly. Because AVX does not have int->double
208 * conversion, we convert via float. */
209 return _mm256_cmp_pd(_mm256_castps_pd(_mm256_cvtepi32_ps(_mm256_castpd_si256(_mm256_and_pd(m0, m1)))), _mm256_setzero_pd(), 0x0c);
212 #endif /* _nbnxn_kernel_simd_utils_x86_s256d_h_ */