Add Power/PowerPC VMX SIMD support
[alexxy/gromacs.git] / src / gromacs / simd / impl_ibm_vmx / impl_ibm_vmx.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 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.
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 #ifndef GMX_SIMD_IMPLEMENTATION_IBM_VMX_H
37 #define GMX_SIMD_IMPLEMENTATION_IBM_VMX_H
38
39 #include <math.h>
40
41 #include <altivec.h>
42
43 /* Make sure we do not screw up c++ - undefine vector/bool, and rely on __vector and __bool */
44 #undef vector
45 #undef bool
46
47 /* IBM VMX SIMD instruction wrappers. Power6 and later.
48  *
49  * Please see documentation in gromacs/simd/simd.h for the available
50  * defines.
51  */
52 /* Capability definitions for IBM VMX */
53 #define GMX_SIMD_HAVE_FLOAT
54 #undef  GMX_SIMD_HAVE_DOUBLE
55 #define GMX_SIMD_HAVE_HARDWARE
56 #undef  GMX_SIMD_HAVE_LOADU
57 #undef  GMX_SIMD_HAVE_STOREU
58 #define GMX_SIMD_HAVE_LOGICAL
59 /* VMX only provides fmadd/fnmadd (our definitions), but not fmsub/fnmsub.
60  * However, fnmadd is what we need for 1/sqrt(x).
61  */
62 #define GMX_SIMD_HAVE_FMA
63 #undef  GMX_SIMD_HAVE_FRACTION
64 #define GMX_SIMD_HAVE_FINT32
65 #undef  GMX_SIMD_HAVE_FINT32_EXTRACT
66 #define GMX_SIMD_HAVE_FINT32_LOGICAL
67 #define GMX_SIMD_HAVE_FINT32_ARITHMETICS
68 #undef  GMX_SIMD_HAVE_DINT32
69 #undef  GMX_SIMD_HAVE_DINT32_EXTRACT
70 #undef  GMX_SIMD_HAVE_DINT32_LOGICAL
71 #undef  GMX_SIMD_HAVE_DINT32_ARITHMETICS
72 #define GMX_SIMD4_HAVE_FLOAT
73 #undef  GMX_SIMD4_HAVE_DOUBLE
74
75 /* Implementation details */
76 #define GMX_SIMD_FLOAT_WIDTH         4
77 #undef  GMX_SIMD_DOUBLE_WIDTH
78 #define GMX_SIMD_FINT32_WIDTH        4
79 #undef  GMX_SIMD_DINT32_WIDTH
80 #define GMX_SIMD_RSQRT_BITS         14
81 #define GMX_SIMD_RCP_BITS           14
82
83 /****************************************************
84  *      SINGLE PRECISION SIMD IMPLEMENTATION        *
85  ****************************************************/
86 #define gmx_simd_float_t           __vector float
87 #define gmx_simd_load_f(m)         vec_ld(0, (const __vector float *)m)
88 #define gmx_simd_load1_f(m)        gmx_simd_load1_f_ibm_vmx(m)
89 #define gmx_simd_set1_f(x)         gmx_simd_set1_f_ibm_vmx(x)
90 #define gmx_simd_store_f(m, x)     vec_st(x, 0, (__vector float *)m)
91 #undef  gmx_simd_loadu_f
92 #undef  gmx_simd_storeu_f
93 #define gmx_simd_setzero_f()       ((__vector float)vec_splat_u32(0))
94 #define gmx_simd_add_f(a, b)       vec_add(a, b)
95 #define gmx_simd_sub_f(a, b)       vec_sub(a, b)
96 #define gmx_simd_mul_f(a, b)       vec_mul(a, b)
97 #define gmx_simd_fmadd_f(a, b, c)  vec_madd(a, b, c)
98 #define gmx_simd_fmsub_f(a, b, c)  vec_sub(vec_mul(a, b), c)
99 /* IBM uses an alternative FMA definition, so -a*b+c=-(a*b-c) is "nmsub" */
100 #define gmx_simd_fnmadd_f(a, b, c) vec_nmsub(a, b, c)
101 /* IBM uses an alternative FMA definition, so -a*b-c=-(a*b+c) is "nmadd" */
102 #define gmx_simd_fnmsub_f(a, b, c) vec_sub(gmx_simd_setzero_f(), vec_madd(a, b, c))
103 #define gmx_simd_and_f(a, b)       vec_and(a, b)
104 #define gmx_simd_andnot_f(a, b)    vec_andc(b, a)
105 #define gmx_simd_or_f(a, b)        vec_or(a, b)
106 #define gmx_simd_xor_f(a, b)       vec_xor(a, b)
107 #define gmx_simd_rsqrt_f(a)        vec_rsqrte(a)
108 #define gmx_simd_rcp_f(a)          vec_re(a)
109 #define gmx_simd_fabs_f(a)         vec_abs(a)
110 #define gmx_simd_fneg_f(a)         vec_xor(a, (__vector float)vec_sl(vec_splat_u32(-1), vec_splat_u32(-1)))
111 #define gmx_simd_max_f(a, b)       vec_max(a, b)
112 #define gmx_simd_min_f(a, b)       vec_min(a, b)
113 #define gmx_simd_round_f(a)        vec_round(a)
114 #define gmx_simd_trunc_f(a)        vec_trunc(a)
115 #define gmx_simd_fraction_f(x)     vec_sub(x, vec_trunc(x))
116 #define gmx_simd_get_exponent_f(a) gmx_simd_get_exponent_f_ibm_vmx(a)
117 #define gmx_simd_get_mantissa_f(a) gmx_simd_get_mantissa_f_ibm_vmx(a)
118 #define gmx_simd_set_exponent_f(a) gmx_simd_set_exponent_f_ibm_vmx(a)
119 /* integer datatype corresponding to float: gmx_simd_fint32_t */
120 #define gmx_simd_fint32_t          __vector int
121 #define gmx_simd_load_fi(m)        vec_ld(0, (const __vector int *)m)
122 #define gmx_simd_set1_fi(i)        gmx_simd_set1_fi_ibm_vmx((int)i)
123 #define gmx_simd_store_fi(m, x)    vec_st(x, 0, (__vector int *)m)
124 #undef  gmx_simd_loadu_fi
125 #undef  gmx_simd_storeu_fi
126 #define gmx_simd_setzero_fi()      vec_splat_s32(0)
127 #define gmx_simd_cvt_f2i(a)        vec_cts(vec_round(a), 0)
128 #define gmx_simd_cvtt_f2i(a)       vec_cts(a, 0)
129 #define gmx_simd_cvt_i2f(a)        vec_ctf(a, 0)
130 #undef  gmx_simd_extract_fi
131 /* Integer logical ops on gmx_simd_fint32_t */
132 /* The shift constant magic requires an explanation:
133  * VMX only allows literals up to 15 to be created directly with vec_splat_u32,
134  * and we need to be able to shift up to 31 bits. The code on the right hand
135  * side splits the constant in three parts with values in the range 0..15.
136  * Since the argument has to be a constant (but our and VMX requirement), these
137  * constants will be evaluated at compile-time, and if one or two parts evaluate
138  * to zero they will be removed with -O2 or higher optimization (checked).
139  */
140 #define gmx_simd_slli_fi(a, i)      vec_sl(a, vec_add(vec_add(vec_splat_u32( (((i&0xF)+(i/16))&0xF)+i/31 ), vec_splat_u32( (i/16)*15 )), vec_splat_u32( (i/31)*15 )))
141 #define gmx_simd_srli_fi(a, i)      vec_sr(a, vec_add(vec_add(vec_splat_u32( (((i&0xF)+(i/16))&0xF)+i/31 ), vec_splat_u32( (i/16)*15 )), vec_splat_u32( (i/31)*15 )))
142 #define gmx_simd_and_fi(a, b)       vec_and(a, b)
143 #define gmx_simd_andnot_fi(a, b)   vec_andc(b, a)
144 #define gmx_simd_or_fi(a, b)        vec_or(a, b)
145 #define gmx_simd_xor_fi(a, b)       vec_xor(a, b)
146 /* Integer arithmetic ops on gmx_simd_fint32_t */
147 #define gmx_simd_add_fi(a, b)       vec_add(a, b)
148 #define gmx_simd_sub_fi(a, b)       vec_sub(a, b)
149 #define gmx_simd_mul_fi(a, b)       vec_mule((__vector short)a, (__vector short)b)
150 /* Boolean & comparison operations on gmx_simd_float_t */
151 #define gmx_simd_fbool_t           __vector __bool int
152 #define gmx_simd_cmpeq_f(a, b)     vec_cmpeq(a, b)
153 #define gmx_simd_cmplt_f(a, b)     vec_cmplt(a, b)
154 #define gmx_simd_cmple_f(a, b)     vec_cmple(a, b)
155 #define gmx_simd_and_fb(a, b)      vec_and(a, b)
156 #define gmx_simd_or_fb(a, b)       vec_or(a, b)
157 #define gmx_simd_anytrue_fb(a)     vec_any_ne(a, (__vector __bool int)vec_splat_u32(0))
158 #define gmx_simd_blendzero_f(a, sel)    vec_and(a, (__vector float)sel)
159 #define gmx_simd_blendnotzero_f(a, sel) vec_andc(a, (__vector float)sel)
160 #define gmx_simd_blendv_f(a, b, sel)    vec_sel(a, b, sel)
161 #define gmx_simd_reduce_f(a)       gmx_simd_reduce_f_ibm_vmx(a)
162 /* Boolean & comparison operations on gmx_simd_fint32_t */
163 #define gmx_simd_fibool_t          __vector __bool int
164 #define gmx_simd_cmpeq_fi(a, b)     vec_cmpeq(a, b)
165 #define gmx_simd_cmplt_fi(a, b)     vec_cmplt(a, b)
166 #define gmx_simd_and_fib(a, b)      vec_and(a, b)
167 #define gmx_simd_or_fib(a, b)       vec_or(a, b)
168 #define gmx_simd_anytrue_fib(a)          vec_any_ne(a, (__vector __bool int)vec_splat_u32(0))
169 #define gmx_simd_blendzero_fi(a, sel)    vec_and(a, (__vector int)sel)
170 #define gmx_simd_blendnotzero_fi(a, sel) vec_andc(a, (__vector int)sel)
171 #define gmx_simd_blendv_fi(a, b, sel)    vec_sel(a, b, sel)
172 /* Conversions between different booleans */
173 #define gmx_simd_cvt_fb2fib(x)     (x)
174 #define gmx_simd_cvt_fib2fb(x)     (x)
175
176 /* Double is not available with VMX SIMD */
177
178 /****************************************************
179  * IMPLEMENTATION HELPER FUNCTIONS                  *
180  ****************************************************/
181 static gmx_inline gmx_simd_float_t
182 gmx_simd_set1_f_ibm_vmx(const float x)
183 {
184     /* In the old days when PPC was all big endian we could
185      * use the vec_lvsl() instruction to permute bytes based on
186      * a load adress. However, at least with gcc-4.8.2 the bytes
187      * end up reversed on Power8 running little endian (Linux).
188      * Since this is not a critical instruction we work around
189      * it by first putting the data in an aligned position before
190      * loading, so we can avoid vec_lvsl() entirely. We can
191      * do this slightly faster on GCC with alignment attributes.
192      */
193     __vector float vx;
194 #ifdef __GNUC__
195     float alignedx __attribute ((aligned (16)));
196     alignedx = x;
197     vx       = vec_lde(0, &alignedx);
198 #else
199     struct {
200         vector float vx; float x[4];
201     } conv;
202     conv.x[0] = x;
203     vx        = vec_lde(0, conv.x);
204 #endif
205     return vec_splat(vx, 0);
206 }
207
208 static gmx_inline gmx_simd_float_t
209 gmx_simd_load1_f_ibm_vmx(const float * m)
210 {
211     return gmx_simd_set1_f_ibm_vmx(*m);
212 }
213
214 static gmx_inline gmx_simd_fint32_t
215 gmx_simd_set1_fi_ibm_vmx(const int x)
216 {
217     /* See comment in gmx_simd_set1_f_ibm_vmx why we
218      * cannot use vec_lvsl().
219      */
220     __vector int vx;
221 #ifdef __GNUC__
222     int alignedx __attribute ((aligned (16)));
223     alignedx = x;
224     vx       = vec_lde(0, &alignedx);
225 #else
226     struct {
227         vector int vx; int x[4];
228     } conv;
229     conv.x[0] = x;
230     vx        = vec_lde(0, conv.x);
231 #endif
232     return vec_splat(vx, 0);
233 }
234
235
236 static gmx_inline gmx_simd_float_t
237 gmx_simd_get_exponent_f_ibm_vmx(gmx_simd_float_t x)
238 {
239     /* Generate 0x7F800000 without memory operations */
240     gmx_simd_float_t     expmask = (__vector float)gmx_simd_slli_fi(vec_add(vec_splat_s32(15), vec_sl(vec_splat_s32(15), vec_splat_u32(4))), 23);
241     gmx_simd_fint32_t    i127    = vec_sub(vec_sl(vec_splat_s32(1), vec_splat_u32(7)), vec_splat_s32(1));
242     gmx_simd_fint32_t    iexp;
243
244     iexp = (__vector int)gmx_simd_and_f(x, expmask);
245     iexp = vec_sub(gmx_simd_srli_fi(iexp, 23), i127);
246     return vec_ctf(iexp, 0);
247 }
248
249 static gmx_inline gmx_simd_float_t
250 gmx_simd_get_mantissa_f_ibm_vmx(gmx_simd_float_t x)
251 {
252     gmx_simd_float_t     expmask = (__vector float)gmx_simd_slli_fi(vec_add(vec_splat_s32(15), vec_sl(vec_splat_s32(15), vec_splat_u32(4))), 23);
253
254     /* Get mantissa. By taking the absolute value (to get rid of the sign bit) we can
255      * use the same mask as for gmx_simd_get_exponent_f() (but complement it). Since
256      * these two routines are typically called together, this will save a few operations.
257      */
258     x = gmx_simd_andnot_f(expmask, vec_abs(x));
259     /* Reset zero (but correctly biased) exponent */
260     return gmx_simd_or_f(x, vec_ctf(vec_splat_s32(1), 0));
261 }
262
263 static gmx_inline gmx_simd_float_t
264 gmx_simd_set_exponent_f_ibm_vmx(gmx_simd_float_t x)
265 {
266     gmx_simd_fint32_t  iexp = gmx_simd_cvt_f2i(x);
267     gmx_simd_fint32_t  i127 = vec_sub(vec_sl(vec_splat_s32(1), vec_splat_u32(7)), vec_splat_s32(1));
268
269     iexp = gmx_simd_slli_fi(vec_add(iexp, i127), 23);
270     return (__vector float)iexp;
271 }
272
273 static gmx_inline float
274 gmx_simd_reduce_f_ibm_vmx(gmx_simd_float_t x)
275 {
276     float res;
277     x = vec_add(x, vec_sld(x, x, 8));
278     x = vec_add(x, vec_sld(x, x, 4));
279     vec_ste(x, 0, &res);
280     return res;
281 }
282
283
284
285 /* SINGLE */
286 #define gmx_simd4_float_t                gmx_simd_float_t
287 #define gmx_simd4_load_f                 gmx_simd_load_f
288 #define gmx_simd4_load1_f                gmx_simd_load1_f
289 #define gmx_simd4_set1_f                 gmx_simd_set1_f
290 #define gmx_simd4_store_f                gmx_simd_store_f
291 #define gmx_simd4_loadu_f                gmx_simd_loadu_f
292 #define gmx_simd4_storeu_f               gmx_simd_storeu_f
293 #define gmx_simd4_setzero_f              gmx_simd_setzero_f
294 #define gmx_simd4_add_f                  gmx_simd_add_f
295 #define gmx_simd4_sub_f                  gmx_simd_sub_f
296 #define gmx_simd4_mul_f                  gmx_simd_mul_f
297 #define gmx_simd4_fmadd_f                gmx_simd_fmadd_f
298 #define gmx_simd4_fmsub_f                gmx_simd_fmsub_f
299 #define gmx_simd4_fnmadd_f               gmx_simd_fnmadd_f
300 #define gmx_simd4_fnmsub_f               gmx_simd_fnmsub_f
301 #define gmx_simd4_and_f                  gmx_simd_and_f
302 #define gmx_simd4_andnot_f               gmx_simd_andnot_f
303 #define gmx_simd4_or_f                   gmx_simd_or_f
304 #define gmx_simd4_xor_f                  gmx_simd_xor_f
305 #define gmx_simd4_rsqrt_f                gmx_simd_rsqrt_f
306 #define gmx_simd4_rcp_f                  gmx_simd_rcp_f
307 #define gmx_simd4_fabs_f                 gmx_simd_fabs_f
308 #define gmx_simd4_fneg_f                 gmx_simd_fneg_f
309 #define gmx_simd4_max_f                  gmx_simd_max_f
310 #define gmx_simd4_min_f                  gmx_simd_min_f
311 #define gmx_simd4_round_f                gmx_simd_round_f
312 #define gmx_simd4_trunc_f                gmx_simd_trunc_f
313 #define gmx_simd4_fraction_f             gmx_simd_fraction_f
314 #define gmx_simd4_get_exponent_f         gmx_simd_get_exponent_f
315 #define gmx_simd4_get_mantissa_f         gmx_simd_get_mantissa_f
316 #define gmx_simd4_set_exponent_f         gmx_simd_set_exponent_f
317 #define gmx_simd4_dotproduct3_f          gmx_simd4_dotproduct3_f_ibm_vmx
318 #define gmx_simd4_fint32_t               gmx_simd_fint32_t
319 #define gmx_simd4_load_fi                gmx_simd_load_fi
320 #define gmx_simd4_load1_fi               gmx_simd_load1_fi
321 #define gmx_simd4_set1_fi                gmx_simd_set1_fi
322 #define gmx_simd4_store_fi               gmx_simd_store_fi
323 #define gmx_simd4_loadu_fi               gmx_simd_loadu_fi
324 #define gmx_simd4_storeu_fi              gmx_simd_storeu_fi
325 #define gmx_simd4_setzero_fi             gmx_simd_setzero_fi
326 #define gmx_simd4_cvt_f2i                gmx_simd_cvt_f2i
327 #define gmx_simd4_cvtt_f2i               gmx_simd_cvtt_f2i
328 #define gmx_simd4_cvt_i2f                gmx_simd_cvt_i2f
329 #define gmx_simd4_fbool_t                gmx_simd_fbool_t
330 #define gmx_simd4_cmpeq_f                gmx_simd_cmpeq_f
331 #define gmx_simd4_cmplt_f                gmx_simd_cmplt_f
332 #define gmx_simd4_cmple_f                gmx_simd_cmple_f
333 #define gmx_simd4_and_fb                 gmx_simd_and_fb
334 #define gmx_simd4_or_fb                  gmx_simd_or_fb
335 #define gmx_simd4_anytrue_fb             gmx_simd_anytrue_fb
336 #define gmx_simd4_blendzero_f            gmx_simd_blendzero_f
337 #define gmx_simd4_blendnotzero_f         gmx_simd_blendnotzero_f
338 #define gmx_simd4_blendv_f               gmx_simd_blendv_f
339 #define gmx_simd4_reduce_f               gmx_simd_reduce_f
340
341 static gmx_inline float
342 gmx_simd4_dotproduct3_f_ibm_vmx(gmx_simd4_float_t a, gmx_simd4_float_t b)
343 {
344     gmx_simd4_float_t c = vec_mul(a, b);
345     /* Keep only elements 0,1,2 by shifting in zero from right */
346     c = vec_sld(c, gmx_simd_setzero_f(), 4);
347     /* calculate sum */
348     return gmx_simd_reduce_f_ibm_vmx(c);
349 }
350
351 /* Function to check whether SIMD operations have resulted in overflow.
352  * For now, this is unfortunately a dummy for this architecture.
353  */
354 static int
355 gmx_simd_check_and_reset_overflow(void)
356 {
357     return 0;
358 }
359
360 #endif /* GMX_SIMD_IMPLEMENTATION_IBM_VMX_H */