Uncrustify all files
[alexxy/gromacs.git] / src / gromacs / simd / macros.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013, 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.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35
36 /* The macros in this file are intended to be used for writing
37  * architecture-independent SIMD intrinsics code.
38  * To support a new architecture, adding macros here should be (nearly)
39  * all that is needed.
40  */
41
42 #ifdef GMX_SIMD_MACROS_H
43 #error "gromacs/simd/macros.h included twice"
44 #else
45 #define GMX_SIMD_MACROS_H
46
47 /* NOTE: SSE2 acceleration does not include floor or blendv */
48
49
50 /* Uncomment the next line, without other SIMD active, for testing plain-C */
51 /* #define GMX_SIMD_REFERENCE_PLAIN_C */
52 #ifdef GMX_SIMD_REFERENCE_PLAIN_C
53 /* Plain C SIMD reference implementation, also serves as documentation */
54 #define GMX_HAVE_SIMD_MACROS
55
56 /* In general the reference SIMD supports any SIMD width, including 1.
57  * See types/nb_verlet.h for details
58  */
59 #define GMX_SIMD_REF_WIDTH  4
60
61 /* Include plain-C reference implementation, also serves as documentation */
62 #include "gromacs/simd/macros_ref.h"
63
64 #define GMX_SIMD_WIDTH_HERE  GMX_SIMD_REF_WIDTH
65
66 /* float/double SIMD register type */
67 #define gmx_mm_pr  gmx_simd_ref_pr
68
69 /* boolean SIMD register type */
70 #define gmx_mm_pb  gmx_simd_ref_pb
71
72 /* integer SIMD register type, only for table indexing and exclusion masks */
73 #define gmx_epi32  gmx_simd_ref_epi32
74 #define GMX_SIMD_EPI32_WIDTH  GMX_SIMD_REF_EPI32_WIDTH
75
76 /* Load GMX_SIMD_WIDTH_HERE reals for memory starting at r */
77 #define gmx_load_pr       gmx_simd_ref_load_pr
78 /* Set all SIMD register elements to *r */
79 #define gmx_load1_pr      gmx_simd_ref_load1_pr
80 #define gmx_set1_pr       gmx_simd_ref_set1_pr
81 #define gmx_setzero_pr    gmx_simd_ref_setzero_pr
82 #define gmx_store_pr      gmx_simd_ref_store_pr
83
84 #define gmx_add_pr        gmx_simd_ref_add_pr
85 #define gmx_sub_pr        gmx_simd_ref_sub_pr
86 #define gmx_mul_pr        gmx_simd_ref_mul_pr
87 /* For the FMA macros below, aim for c=d in code, so FMA3 uses 1 instruction */
88 #define gmx_madd_pr       gmx_simd_ref_madd_pr
89 #define gmx_nmsub_pr      gmx_simd_ref_nmsub_pr
90
91 #define gmx_max_pr        gmx_simd_ref_max_pr
92 #define gmx_blendzero_pr  gmx_simd_ref_blendzero_pr
93
94 #define gmx_round_pr      gmx_simd_ref_round_pr
95
96 /* Not required, only used to speed up the nbnxn tabulated PME kernels */
97 #define GMX_SIMD_HAVE_FLOOR
98 #ifdef GMX_SIMD_HAVE_FLOOR
99 #define gmx_floor_pr      gmx_simd_ref_floor_pr
100 #endif
101
102 /* Not required, only used when blendv is faster than comparison */
103 #define GMX_SIMD_HAVE_BLENDV
104 #ifdef GMX_SIMD_HAVE_BLENDV
105 #define gmx_blendv_pr     gmx_simd_ref_blendv_pr
106 #endif
107
108 /* Copy the sign of a to b, assumes b >= 0 for efficiency */
109 #define gmx_cpsgn_nonneg_pr  gmx_simd_ref_cpsgn_nonneg_pr
110
111 /* Very specific operation required in the non-bonded kernels */
112 #define gmx_masknot_add_pr   gmx_simd_ref_masknot_add_pr
113
114 /* Comparison */
115 #define gmx_cmplt_pr      gmx_simd_ref_cmplt_pr
116
117 /* Logical operations on SIMD booleans */
118 #define gmx_and_pb        gmx_simd_ref_and_pb
119 #define gmx_or_pb         gmx_simd_ref_or_pb
120
121 /* Returns a single int (0/1) which tells if any of the 4 booleans is True */
122 #define gmx_anytrue_pb    gmx_simd_ref_anytrue_pb
123
124 /* Conversions only used for PME table lookup */
125 #define gmx_cvttpr_epi32  gmx_simd_ref_cvttpr_epi32
126 #define gmx_cvtepi32_pr   gmx_simd_ref_cvtepi32_pr
127
128 /* These two function only need to be approximate, Newton-Raphson iteration
129  * is used for full accuracy in gmx_invsqrt_pr and gmx_inv_pr.
130  */
131 #define gmx_rsqrt_pr      gmx_simd_ref_rsqrt_pr
132 #define gmx_rcp_pr        gmx_simd_ref_rcp_pr
133
134 /* sqrt+inv+sin+cos+acos+atan2 are used for bonded potentials, exp for PME */
135 #define GMX_SIMD_HAVE_EXP
136 #ifdef GMX_SIMD_HAVE_EXP
137 #define gmx_exp_pr        gmx_simd_ref_exp_pr
138 #endif
139 #define GMX_SIMD_HAVE_TRIGONOMETRIC
140 #ifdef GMX_SIMD_HAVE_TRIGONOMETRIC
141 #define gmx_sqrt_pr       gmx_simd_ref_sqrt_pr
142 #define gmx_sincos_pr     gmx_simd_ref_sincos_pr
143 #define gmx_acos_pr       gmx_simd_ref_acos_pr
144 #define gmx_atan2_pr      gmx_simd_ref_atan2_pr
145 #endif
146
147 #endif /* GMX_SIMD_REFERENCE_PLAIN_C */
148
149
150 /* The same SIMD macros can be translated to SIMD intrinsics (and compiled
151  * to instructions for) different SIMD width and float precision.
152  *
153  * On x86: The gmx_ prefix is replaced by _mm_ or _mm256_ (SSE or AVX).
154  * The _pr suffix is replaced by _ps or _pd (for single or double precision).
155  * Compiler settings will decide if 128-bit intrinsics will
156  * be translated into SSE or AVX instructions.
157  */
158
159
160 #ifdef GMX_USE_HALF_WIDTH_SIMD_HERE
161 #if defined GMX_X86_AVX_256 || defined __MIC__
162 /* We have half SIMD width support, continue */
163 #else
164 #error "half SIMD width intrinsics are not supported"
165 #endif
166 #endif
167
168 #if defined GMX_TARGET_X86 && !defined __MIC__
169
170 #ifdef GMX_X86_SSE2
171 /* This is for general x86 SIMD instruction sets that also support SSE2 */
172 #define GMX_HAVE_SIMD_MACROS
173
174 /* Include the highest supported x86 SIMD intrisics + math functions */
175 #ifdef GMX_X86_AVX_256
176 #include "general_x86_avx_256.h"
177 #ifdef GMX_DOUBLE
178 #include "math_x86_avx_256_double.h"
179 #else  /* GMX_DOUBLE */
180 #include "math_x86_avx_256_single.h"
181 #endif /* GMX_DOUBLE */
182 #else  /* GMX_X86_AVX_256 */
183 #ifdef GMX_X86_AVX_128_FMA
184 #include "general_x86_avx_128_fma.h"
185 #ifdef GMX_DOUBLE
186 #include "math_x86_avx_128_fma_double.h"
187 #else  /* GMX_DOUBLE */
188 #include "math_x86_avx_128_fma_single.h"
189 #endif /* GMX_DOUBLE */
190 #else  /* GMX_X86_AVX_128_FMA */
191 #ifdef GMX_X86_SSE4_1
192 #include "general_x86_sse4_1.h"
193 #ifdef GMX_DOUBLE
194 #include "math_x86_sse4_1_double.h"
195 #else  /* GMX_DOUBLE */
196 #include "math_x86_sse4_1_single.h"
197 #endif /* GMX_DOUBLE */
198 #else  /* GMX_X86_SSE4_1 */
199 #ifdef GMX_X86_SSE2
200 #include "general_x86_sse2.h"
201 #ifdef GMX_DOUBLE
202 #include "math_x86_sse2_double.h"
203 #else  /* GMX_DOUBLE */
204 #include "math_x86_sse2_single.h"
205 #endif /* GMX_DOUBLE */
206 #else  /* GMX_X86_SSE2 */
207 #error No x86 acceleration defined
208 #endif /* GMX_X86_SSE2 */
209 #endif /* GMX_X86_SSE4_1 */
210 #endif /* GMX_X86_AVX_128_FMA */
211 #endif /* GMX_X86_AVX_256 */
212
213 /* exp and trigonometric functions are included above */
214 #define GMX_SIMD_HAVE_EXP
215 #define GMX_SIMD_HAVE_TRIGONOMETRIC
216
217 #if !defined GMX_X86_AVX_256 || defined GMX_USE_HALF_WIDTH_SIMD_HERE
218
219 #ifndef GMX_DOUBLE
220
221 #define GMX_SIMD_WIDTH_HERE  4
222
223 #define gmx_mm_pr  __m128
224
225 #define gmx_mm_pb  __m128
226
227 #define gmx_epi32  __m128i
228 #define GMX_SIMD_EPI32_WIDTH  4
229
230 #define gmx_load_pr       _mm_load_ps
231 #define gmx_load1_pr      _mm_load1_ps
232 #define gmx_set1_pr       _mm_set1_ps
233 #define gmx_setzero_pr    _mm_setzero_ps
234 #define gmx_store_pr      _mm_store_ps
235
236 #define gmx_add_pr        _mm_add_ps
237 #define gmx_sub_pr        _mm_sub_ps
238 #define gmx_mul_pr        _mm_mul_ps
239 #ifdef GMX_X86_AVX_128_FMA
240 #define GMX_SIMD_HAVE_FMA
241 #define gmx_madd_pr(a, b, c)   _mm_macc_ps(a, b, c)
242 #define gmx_nmsub_pr(a, b, c)  _mm_nmacc_ps(a, b, c)
243 #else
244 #define gmx_madd_pr(a, b, c)   _mm_add_ps(c, _mm_mul_ps(a, b))
245 #define gmx_nmsub_pr(a, b, c)  _mm_sub_ps(c, _mm_mul_ps(a, b))
246 #endif
247 #define gmx_max_pr        _mm_max_ps
248 #define gmx_blendzero_pr  _mm_and_ps
249
250 #define gmx_cmplt_pr      _mm_cmplt_ps
251 #define gmx_and_pb        _mm_and_ps
252 #define gmx_or_pb         _mm_or_ps
253
254 #ifdef GMX_X86_SSE4_1
255 #define gmx_round_pr(x)   _mm_round_ps(x, 0x0)
256 #define GMX_SIMD_HAVE_FLOOR
257 #define gmx_floor_pr      _mm_floor_ps
258 #else
259 #define gmx_round_pr(x)   _mm_cvtepi32_ps(_mm_cvtps_epi32(x))
260 #endif
261
262 #ifdef GMX_X86_SSE4_1
263 #define GMX_SIMD_HAVE_BLENDV
264 #define gmx_blendv_pr     _mm_blendv_ps
265 #endif
266
267 static gmx_inline gmx_mm_pr gmx_cpsgn_nonneg_pr(gmx_mm_pr a, gmx_mm_pr b)
268 {
269     /* The value -0.0 has only the sign-bit set */
270     gmx_mm_pr sign_mask = _mm_set1_ps(-0.0);
271     return _mm_or_ps(_mm_and_ps(a, sign_mask), b);
272 };
273
274 static gmx_inline gmx_mm_pr gmx_masknot_add_pr(gmx_mm_pb a, gmx_mm_pr b, gmx_mm_pr c)
275 {
276     return _mm_add_ps(b, _mm_andnot_ps(a, c));
277 };
278
279 #define gmx_anytrue_pb    _mm_movemask_ps
280
281 #define gmx_cvttpr_epi32  _mm_cvttps_epi32
282 #define gmx_cvtepi32_pr   _mm_cvtepi32_ps
283
284 #define gmx_rsqrt_pr      _mm_rsqrt_ps
285 #define gmx_rcp_pr        _mm_rcp_ps
286
287 #define gmx_exp_pr        gmx_mm_exp_ps
288 #define gmx_sqrt_pr       gmx_mm_sqrt_ps
289 #define gmx_sincos_pr     gmx_mm_sincos_ps
290 #define gmx_acos_pr       gmx_mm_acos_ps
291 #define gmx_atan2_pr      gmx_mm_atan2_ps
292 #define gmx_erfc_pr       gmx_mm_erfc_ps
293
294 #else /* ifndef GMX_DOUBLE */
295
296 #define GMX_SIMD_WIDTH_HERE  2
297
298 #define gmx_mm_pr  __m128d
299
300 #define gmx_mm_pb  __m128d
301
302 #define gmx_epi32  __m128i
303 #define GMX_SIMD_EPI32_WIDTH  4
304
305 #define gmx_load_pr       _mm_load_pd
306 #define gmx_load1_pr      _mm_load1_pd
307 #define gmx_set1_pr       _mm_set1_pd
308 #define gmx_setzero_pr    _mm_setzero_pd
309 #define gmx_store_pr      _mm_store_pd
310
311 #define gmx_add_pr        _mm_add_pd
312 #define gmx_sub_pr        _mm_sub_pd
313 #define gmx_mul_pr        _mm_mul_pd
314 #ifdef GMX_X86_AVX_128_FMA
315 #define GMX_SIMD_HAVE_FMA
316 #define gmx_madd_pr(a, b, c)   _mm_macc_pd(a, b, c)
317 #define gmx_nmsub_pr(a, b, c)  _mm_nmacc_pd(a, b, c)
318 #else
319 #define gmx_madd_pr(a, b, c)   _mm_add_pd(c, _mm_mul_pd(a, b))
320 #define gmx_nmsub_pr(a, b, c)  _mm_sub_pd(c, _mm_mul_pd(a, b))
321 #endif
322 #define gmx_max_pr        _mm_max_pd
323 #define gmx_blendzero_pr  _mm_and_pd
324
325 #ifdef GMX_X86_SSE4_1
326 #define gmx_round_pr(x)   _mm_round_pd(x, 0x0)
327 #define GMX_SIMD_HAVE_FLOOR
328 #define gmx_floor_pr      _mm_floor_pd
329 #else
330 #define gmx_round_pr(x)   _mm_cvtepi32_pd(_mm_cvtpd_epi32(x))
331 /* gmx_floor_pr is not used in code for pre-SSE4_1 hardware */
332 #endif
333
334 #ifdef GMX_X86_SSE4_1
335 #define GMX_SIMD_HAVE_BLENDV
336 #define gmx_blendv_pr     _mm_blendv_pd
337 #endif
338
339 static gmx_inline gmx_mm_pr gmx_cpsgn_nonneg_pr(gmx_mm_pr a, gmx_mm_pr b)
340 {
341     gmx_mm_pr sign_mask = _mm_set1_pd(-0.0);
342     return _mm_or_pd(_mm_and_pd(a, sign_mask), b);
343 };
344
345 static gmx_inline gmx_mm_pr gmx_masknot_add_pr(gmx_mm_pb a, gmx_mm_pr b, gmx_mm_pr c)
346 {
347     return _mm_add_pd(b, _mm_andnot_pd(a, c));
348 };
349
350 #define gmx_cmplt_pr      _mm_cmplt_pd
351
352 #define gmx_and_pb        _mm_and_pd
353 #define gmx_or_pb         _mm_or_pd
354
355 #define gmx_anytrue_pb    _mm_movemask_pd
356
357 #define gmx_cvttpr_epi32  _mm_cvttpd_epi32
358 #define gmx_cvtepi32_pr   _mm_cvtepi32_pd
359
360 #define gmx_rsqrt_pr(r)   _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(r)))
361 #define gmx_rcp_pr(r)     _mm_cvtps_pd(_mm_rcp_ps(_mm_cvtpd_ps(r)))
362
363 #define gmx_exp_pr        gmx_mm_exp_pd
364 #define gmx_sqrt_pr       gmx_mm_sqrt_pd
365 #define gmx_sincos_pr     gmx_mm_sincos_pd
366 #define gmx_acos_pr       gmx_mm_acos_pd
367 #define gmx_atan2_pr      gmx_mm_atan2_pd
368 #define gmx_erfc_pr       gmx_mm_erfc_pd
369
370 #endif /* ifndef GMX_DOUBLE */
371
372 #else
373 /* We have GMX_X86_AVX_256 and not GMX_USE_HALF_WIDTH_SIMD_HERE,
374  * so we use 256-bit SIMD.
375  */
376
377 #ifndef GMX_DOUBLE
378
379 #define GMX_SIMD_WIDTH_HERE  8
380
381 #define gmx_mm_pr  __m256
382
383 #define gmx_mm_pb  __m256
384
385 #define gmx_epi32  __m256i
386 #define GMX_SIMD_EPI32_WIDTH  8
387
388 #define gmx_load_pr       _mm256_load_ps
389 #define gmx_load1_pr(x)   _mm256_set1_ps((x)[0])
390 #define gmx_set1_pr       _mm256_set1_ps
391 #define gmx_setzero_pr    _mm256_setzero_ps
392 #define gmx_store_pr      _mm256_store_ps
393
394 #define gmx_add_pr        _mm256_add_ps
395 #define gmx_sub_pr        _mm256_sub_ps
396 #define gmx_mul_pr        _mm256_mul_ps
397 #define gmx_madd_pr(a, b, c)   _mm256_add_ps(c, _mm256_mul_ps(a, b))
398 #define gmx_nmsub_pr(a, b, c)  _mm256_sub_ps(c, _mm256_mul_ps(a, b))
399 #define gmx_max_pr        _mm256_max_ps
400 #define gmx_blendzero_pr  _mm256_and_ps
401
402 #define gmx_round_pr(x)   _mm256_round_ps(x, 0x0)
403 #define GMX_SIMD_HAVE_FLOOR
404 #define gmx_floor_pr      _mm256_floor_ps
405
406 #define GMX_SIMD_HAVE_BLENDV
407 #define gmx_blendv_pr     _mm256_blendv_ps
408
409 static gmx_inline gmx_mm_pr gmx_cpsgn_nonneg_pr(gmx_mm_pr a, gmx_mm_pr b)
410 {
411     gmx_mm_pr sign_mask = _mm256_set1_ps(-0.0);
412     return _mm256_or_ps(_mm256_and_ps(a, sign_mask), b);
413 };
414
415 static gmx_inline gmx_mm_pr gmx_masknot_add_pr(gmx_mm_pb a, gmx_mm_pr b, gmx_mm_pr c)
416 {
417     return _mm256_add_ps(b, _mm256_andnot_ps(a, c));
418 };
419
420 /* Less-than (we use ordered, non-signaling, but that's not required) */
421 #define gmx_cmplt_pr(x, y) _mm256_cmp_ps(x, y, 0x11)
422 #define gmx_and_pb        _mm256_and_ps
423 #define gmx_or_pb         _mm256_or_ps
424
425 #define gmx_anytrue_pb    _mm256_movemask_ps
426
427 #define gmx_cvttpr_epi32  _mm256_cvttps_epi32
428
429 #define gmx_rsqrt_pr      _mm256_rsqrt_ps
430 #define gmx_rcp_pr        _mm256_rcp_ps
431
432 #define gmx_exp_pr        gmx_mm256_exp_ps
433 #define gmx_sqrt_pr       gmx_mm256_sqrt_ps
434 #define gmx_sincos_pr     gmx_mm256_sincos_ps
435 #define gmx_acos_pr       gmx_mm256_acos_ps
436 #define gmx_atan2_pr      gmx_mm256_atan2_ps
437 #define gmx_erfc_pr       gmx_mm256_erfc_ps
438
439 #else /* ifndef GMX_DOUBLE */
440
441 #define GMX_SIMD_WIDTH_HERE  4
442
443 #define gmx_mm_pr  __m256d
444
445 #define gmx_mm_pb  __m256d
446
447 /* We use 128-bit integer registers because of missing 256-bit operations */
448 #define gmx_epi32  __m128i
449 #define GMX_SIMD_EPI32_WIDTH  4
450
451 #define gmx_load_pr       _mm256_load_pd
452 #define gmx_load1_pr(x)   _mm256_set1_pd((x)[0])
453 #define gmx_set1_pr       _mm256_set1_pd
454 #define gmx_setzero_pr    _mm256_setzero_pd
455 #define gmx_store_pr      _mm256_store_pd
456
457 #define gmx_add_pr        _mm256_add_pd
458 #define gmx_sub_pr        _mm256_sub_pd
459 #define gmx_mul_pr        _mm256_mul_pd
460 #define gmx_madd_pr(a, b, c)   _mm256_add_pd(c, _mm256_mul_pd(a, b))
461 #define gmx_nmsub_pr(a, b, c)  _mm256_sub_pd(c, _mm256_mul_pd(a, b))
462 #define gmx_max_pr        _mm256_max_pd
463 #define gmx_blendzero_pr  _mm256_and_pd
464
465 #define gmx_round_pr(x)   _mm256_round_pd(x, 0x0)
466 #define GMX_SIMD_HAVE_FLOOR
467 #define gmx_floor_pr      _mm256_floor_pd
468
469 #define GMX_SIMD_HAVE_BLENDV
470 #define gmx_blendv_pr     _mm256_blendv_pd
471
472 static gmx_inline gmx_mm_pr gmx_cpsgn_nonneg_pr(gmx_mm_pr a, gmx_mm_pr b)
473 {
474     gmx_mm_pr sign_mask = _mm256_set1_pd(-0.0);
475     return _mm256_or_pd(_mm256_and_pd(a, sign_mask), b);
476 };
477
478 static gmx_inline gmx_mm_pr gmx_masknot_add_pr(gmx_mm_pb a, gmx_mm_pr b, gmx_mm_pr c)
479 {
480     return _mm256_add_pd(b, _mm256_andnot_pd(a, c));
481 };
482
483 /* Less-than (we use ordered, non-signaling, but that's not required) */
484 #define gmx_cmplt_pr(x, y) _mm256_cmp_pd(x, y, 0x11)
485
486 #define gmx_and_pb        _mm256_and_pd
487 #define gmx_or_pb         _mm256_or_pd
488
489 #define gmx_anytrue_pb    _mm256_movemask_pd
490
491 #define gmx_cvttpr_epi32  _mm256_cvttpd_epi32
492
493 #define gmx_rsqrt_pr(r)   _mm256_cvtps_pd(_mm_rsqrt_ps(_mm256_cvtpd_ps(r)))
494 #define gmx_rcp_pr(r)     _mm256_cvtps_pd(_mm_rcp_ps(_mm256_cvtpd_ps(r)))
495
496 #define gmx_exp_pr        gmx_mm256_exp_pd
497 #define gmx_sqrt_pr       gmx_mm256_sqrt_pd
498 #define gmx_sincos_pr     gmx_mm256_sincos_pd
499 #define gmx_acos_pr       gmx_mm256_acos_pd
500 #define gmx_atan2_pr      gmx_mm256_atan2_pd
501 #define gmx_erfc_pr       gmx_mm256_erfc_pd
502
503 #endif /* ifndef GMX_DOUBLE */
504
505 #endif /* 128- or 256-bit x86 SIMD */
506
507 #endif /* GMX_X86_SSE2 */
508
509 #endif /* GMX_TARGET_X86 */
510
511 #ifdef GMX_CPU_ACCELERATION_IBM_QPX
512
513 /* This hack works on the compilers that can reach this code. A real
514    solution with broader scope will be proposed in master branch. */
515 #define gmx_always_inline __attribute__((always_inline))
516
517 /* This is for the A2 core on BlueGene/Q that supports IBM's QPX
518    vector built-in functions */
519 #include <mass_simd.h>
520 #define GMX_HAVE_SIMD_MACROS
521 #ifdef __clang__
522 #include <qpxmath.h>
523 #endif
524
525 /* No need to version the code by the precision, because the QPX AXU
526    extends to and truncates from double precision for free. */
527
528 #define GMX_SIMD_WIDTH_HERE  4
529 typedef vector4double gmx_mm_pr;
530 typedef vector4double gmx_mm_pb;
531 typedef vector4double gmx_epi32;
532 #define GMX_SIMD_EPI32_WIDTH  4
533
534 static gmx_inline gmx_mm_pr gmx_always_inline gmx_load_pr(const real *a)
535 {
536 #ifdef NDEBUG
537     return vec_ld(0, (real *) a);
538 #else
539     return vec_lda(0, (real *) a);
540 #endif
541 }
542
543 static gmx_inline gmx_mm_pr gmx_always_inline gmx_load1_pr(const real *a)
544 {
545     return vec_splats(*a);
546 }
547
548 static gmx_inline gmx_mm_pr gmx_always_inline gmx_set1_pr(real a)
549 {
550     return vec_splats(a);
551 }
552
553 static gmx_inline gmx_mm_pr gmx_always_inline gmx_setzero_pr()
554 {
555     return vec_splats(0.0);
556 }
557
558 static gmx_inline void gmx_always_inline gmx_store_pr(real *a, gmx_mm_pr b)
559 {
560 #ifdef NDEBUG
561     vec_st(b, 0, a);
562 #else
563     vec_sta(b, 0, a);
564 #endif
565 }
566
567 static gmx_inline gmx_mm_pr gmx_always_inline gmx_add_pr(gmx_mm_pr a, gmx_mm_pr b)
568 {
569     return vec_add(a, b);
570 }
571
572 static gmx_inline gmx_mm_pr gmx_always_inline gmx_sub_pr(gmx_mm_pr a, gmx_mm_pr b)
573 {
574     return vec_sub(a, b);
575 }
576
577 static gmx_inline gmx_mm_pr gmx_always_inline gmx_mul_pr(gmx_mm_pr a, gmx_mm_pr b)
578 {
579     return vec_mul(a, b);
580 }
581
582 static gmx_inline gmx_mm_pr gmx_always_inline gmx_madd_pr(gmx_mm_pr a, gmx_mm_pr b, gmx_mm_pr c)
583 {
584     return vec_madd(a, b, c);
585 }
586
587 static gmx_inline gmx_mm_pr gmx_always_inline gmx_nmsub_pr(gmx_mm_pr a, gmx_mm_pr b, gmx_mm_pr c)
588 {
589     return vec_nmsub(a, b, c);
590 }
591
592 static gmx_inline gmx_mm_pr gmx_always_inline gmx_max_pr(gmx_mm_pr a, gmx_mm_pr b)
593 {
594     return vec_sel(b, a, vec_sub(a, b));
595 }
596
597 static gmx_inline gmx_mm_pr gmx_always_inline gmx_blendzero_pr(gmx_mm_pr a, gmx_mm_pr b)
598 {
599     return vec_sel(gmx_setzero_pr(), a, b);
600 }
601
602 static gmx_inline gmx_mm_pb gmx_always_inline gmx_cmplt_pr(gmx_mm_pr a, gmx_mm_pr b)
603 {
604     return vec_cmplt(a, b);
605 }
606
607 static gmx_inline gmx_mm_pb gmx_always_inline gmx_and_pb(gmx_mm_pb a, gmx_mm_pb b)
608 {
609     return vec_and(a, b);
610 }
611
612 static gmx_inline gmx_mm_pb gmx_always_inline gmx_or_pb(gmx_mm_pb a, gmx_mm_pb b)
613 {
614     return vec_or(a, b);
615 }
616
617 static gmx_inline gmx_mm_pr gmx_always_inline gmx_round_pr(gmx_mm_pr a)
618 {
619     return vec_round(a);
620 }
621
622 #define GMX_SIMD_HAVE_FLOOR
623 static gmx_inline gmx_mm_pr gmx_always_inline gmx_floor_pr(gmx_mm_pr a)
624 {
625     return vec_floor(a);
626 }
627
628 #define GMX_SIMD_HAVE_BLENDV
629 static gmx_inline gmx_mm_pr gmx_always_inline gmx_blendv_pr(gmx_mm_pr a, gmx_mm_pr b, gmx_mm_pr c)
630 {
631     return vec_sel(b, a, gmx_cmplt_pr(gmx_setzero_pr(), c));
632 }
633
634 static gmx_inline gmx_mm_pr gmx_always_inline gmx_cpsgn_nonneg_pr(gmx_mm_pr a, gmx_mm_pr b)
635 {
636     return vec_cpsgn(a, b);
637 };
638
639 static gmx_inline gmx_mm_pr gmx_always_inline gmx_masknot_add_pr(gmx_mm_pb a, gmx_mm_pr b, gmx_mm_pr c)
640 {
641     return vec_add(b, vec_sel(c, gmx_setzero_pr(), a));
642 };
643
644 static gmx_inline gmx_bool gmx_always_inline
645 GMX_SIMD_IS_TRUE(real x)
646 {
647     return x >= 0.0;
648 }
649
650 static gmx_inline gmx_epi32 gmx_always_inline gmx_cvttpr_epi32(gmx_mm_pr a)
651 {
652     return vec_ctiwuz(a);
653 }
654 /* Don't want this, we have floor */
655 /* #define gmx_cvtepi32_pr   vec_cvtepi32 */
656
657 /* A2 core on BG/Q delivers relative error of 2^-14, whereas Power ISA
658    Architecture only promises 2^-8. So probably no need for
659    Newton-Raphson iterates at single or double. */
660 static gmx_inline gmx_mm_pr gmx_always_inline gmx_rsqrt_pr(gmx_mm_pr a)
661 {
662     return vec_rsqrte(a);
663 }
664
665 /* A2 core on BG/Q delivers relative error of 2^-14, whereas Power ISA
666    Architecture only promises 2^-5. So probably no need for
667    Newton-Raphson iterates at single or double. */
668 static gmx_inline gmx_mm_pr gmx_always_inline gmx_rcp_pr(gmx_mm_pr a)
669 {
670     return vec_re(a);
671 }
672
673 /* Note that here, and below, we use the built-in SLEEF port when
674    compiling on BlueGene/Q with clang */
675
676 #define GMX_SIMD_HAVE_EXP
677 static gmx_inline gmx_mm_pr gmx_always_inline gmx_exp_pr(gmx_mm_pr a)
678 {
679 #ifdef __clang__
680 #ifndef GMX_DOUBLE
681     return xexpf(a);
682 #else
683     return xexp(a);
684 #endif
685 #else
686 #ifndef GMX_DOUBLE
687     return expf4(a);
688 #else
689     return expd4(a);
690 #endif
691 #endif
692 }
693
694 static gmx_inline gmx_mm_pr gmx_always_inline gmx_sqrt_pr(gmx_mm_pr a)
695 {
696 #ifdef NDEBUG
697     return vec_swsqrt_nochk(a);
698 #else
699     return vec_swsqrt(a);
700 #endif
701 }
702
703 #define GMX_SIMD_HAVE_TRIGONOMETRIC
704 static gmx_inline int gmx_always_inline gmx_sincos_pr(gmx_mm_pr a, gmx_mm_pr *b, gmx_mm_pr *c)
705 {
706 #ifdef __clang__
707 #ifndef GMX_DOUBLE
708     xsincosf(a, b, c);
709 #else
710     xsincos(a, b, c);
711 #endif
712 #else
713 #ifndef GMX_DOUBLE
714     sincosf4(a, b, c);
715 #else
716     sincosd4(a, b, c);
717 #endif
718 #endif
719     return 1;
720 }
721
722 static gmx_inline gmx_mm_pr gmx_always_inline gmx_acos_pr(gmx_mm_pr a)
723 {
724 #ifdef __clang__
725 #ifndef GMX_DOUBLE
726     return xacosf(a);
727 #else
728     return xacos(a);
729 #endif
730 #else
731 #ifndef GMX_DOUBLE
732     return acosf4(a);
733 #else
734     return acosd4(a);
735 #endif
736 #endif
737 }
738
739 /* NB The order of parameters here is correct; the
740    documentation of atan2[df]4 in SIMD MASS is wrong. */
741 static gmx_inline gmx_mm_pr gmx_always_inline gmx_atan2_pr(gmx_mm_pr a, gmx_mm_pr b)
742 {
743 #ifdef __clang__
744 #ifndef GMX_DOUBLE
745     return xatan2f(a, b);
746 #else
747     return xatan2(a, b);
748 #endif
749 #else
750 #ifndef GMX_DOUBLE
751     return atan2f4(a, b);
752 #else
753     return atan2d4(a, b);
754 #endif
755 #endif
756 }
757
758 static gmx_inline gmx_mm_pr gmx_always_inline gmx_erfc_pr(gmx_mm_pr a)
759 {
760     /* The BG/Q qpxmath.h vector math library intended for use with
761        bgclang does not have erfc, so we need to use a function from
762        mass_simd.h. If this changes, then the #include <mass_simd.h> can
763        become conditional. */
764 #ifndef GMX_DOUBLE
765     return erfcf4(a);
766 #else
767     return erfcd4(a);
768 #endif
769 }
770
771 /* TODO: gmx_mm_erfc_p[sd] should be generalized using gmx_*_pr, so that it just works on BlueGene */
772
773 static gmx_inline int gmx_always_inline
774 gmx_anytrue_pb(gmx_mm_pb a)
775 {
776     /* The "anytrue" is done solely on the QPX AXU (which is the only
777        available FPU). This is awkward, because pretty much no
778        "horizontal" SIMD-vector operations exist, unlike x86 where
779        SSE4.1 added various kinds of horizontal operations. So we have
780        to make do with shifting vector elements and operating on the
781        results. This makes for lots of data dependency, but the main
782        alternative of storing to memory and reloading is not going to
783        help, either. OpenMP over 2 or 4 hardware threads per core will
784        hide much of the latency from the data dependency. The
785        vec_extract() lets the compiler correctly use a floating-point
786        comparison on the zeroth vector element, which avoids needing
787        memory at all.
788      */
789     gmx_mm_pb vec_shifted_left_0 = a;
790     gmx_mm_pb vec_shifted_left_1 = vec_sldw(a, a, 1);
791     gmx_mm_pb vec_shifted_left_2 = vec_sldw(a, a, 2);
792     gmx_mm_pb vec_shifted_left_3 = vec_sldw(a, a, 3);
793
794     gmx_mm_pb vec_return = vec_or(vec_or(vec_shifted_left_2, vec_shifted_left_3),
795                                   vec_or(vec_shifted_left_0, vec_shifted_left_1));
796     return (0.0 < vec_extract(vec_return, 0));
797 };
798
799 #undef gmx_always_inline
800
801 #endif /* GMX_CPU_ACCELERATION_IBM_QPX */
802
803 #ifdef __MIC__
804 #include "general_x86_mic.h"
805 #endif
806
807 #ifdef GMX_HAVE_SIMD_MACROS
808 /* Generic functions to extract a SIMD aligned pointer from a pointer x.
809  * x should have at least GMX_SIMD_WIDTH_HERE elements extra compared
810  * to how many you want to use, to avoid indexing outside the aligned region.
811  */
812
813 static gmx_inline real *
814 gmx_simd_align_real(const real *x)
815 {
816     return (real *)(((size_t)((x)+GMX_SIMD_WIDTH_HERE)) & (~((size_t)(GMX_SIMD_WIDTH_HERE*sizeof(real)-1))));
817 }
818
819 static gmx_inline int *
820 gmx_simd_align_int(const int *x)
821 {
822     return (int  *)(((size_t)((x)+GMX_SIMD_WIDTH_HERE)) & (~((size_t)(GMX_SIMD_WIDTH_HERE*sizeof(int )-1))));
823 }
824
825
826 /* Include the math functions which only need the above macros,
827  * generally these are the ones that don't need masking operations.
828  */
829 #ifdef GMX_DOUBLE
830 #include "math_double.h"
831 #else
832 #include "math_single.h"
833 #endif
834
835
836 #endif /* GMX_HAVE_SIMD_MACROS */
837
838 #endif