1 /*
2 * Copyright (c) 2003, 2007-14 Matteo Frigo
3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
4 *
5 * Modifications by Romain Dolbeau & Erik Lindahl, derived from simd-avx.h
6 * Romain Dolbeau hereby places his modifications in the public domain.
7 * Erik Lindahl hereby places his modifications in the public domain.
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program 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
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 *
23 */
24
25 #if defined(FFTW_LDOUBLE) || defined(FFTW_QUAD)
26 #error "AVX2 only works in single or double precision"
27 #endif
28
29 #ifdef FFTW_SINGLE
30 # define DS(d,s) s /* single-precision option */
31 # define SUFF(name) name ## s
32 #else
33 # define DS(d,s) d /* double-precision option */
34 # define SUFF(name) name ## d
35 #endif
36
37 #define SIMD_SUFFIX _avx2 /* for renaming */
38 #define VL DS(2, 4) /* SIMD complex vector length */
39 #define SIMD_VSTRIDE_OKA(x) ((x) == 2)
40 #define SIMD_STRIDE_OKPAIR SIMD_STRIDE_OK
41
42 #if defined(__GNUC__) && !defined(__AVX2__) /* sanity check */
43 #error "compiling simd-avx2.h without avx2 support"
44 #endif
45
46 #ifdef _MSC_VER
47 #ifndef inline
48 #define inline __inline
49 #endif
50 #endif
51
52 #include <immintrin.h>
53
54 typedef DS(__m256d, __m256) V;
55 #define VADD SUFF(_mm256_add_p)
56 #define VSUB SUFF(_mm256_sub_p)
57 #define VMUL SUFF(_mm256_mul_p)
58 #define VXOR SUFF(_mm256_xor_p)
59 #define VSHUF SUFF(_mm256_shuffle_p)
60 #define VPERM1 SUFF(_mm256_permute_p)
61
62 #define SHUFVALD(fp0,fp1) \
63 (((fp1) << 3) | ((fp0) << 2) | ((fp1) << 1) | ((fp0)))
64 #define SHUFVALS(fp0,fp1,fp2,fp3) \
65 (((fp3) << 6) | ((fp2) << 4) | ((fp1) << 2) | ((fp0)))
66
67 #define VDUPL(x) DS(_mm256_movedup_pd(x), _mm256_moveldup_ps(x))
68 #define VDUPH(x) DS(_mm256_permute_pd(x,SHUFVALD(1,1)), _mm256_movehdup_ps(x))
69
70 #define VLIT(x0, x1) DS(_mm256_set_pd(x0, x1, x0, x1), _mm256_set_ps(x0, x1, x0, x1, x0, x1, x0, x1))
71 #define DVK(var, val) V var = VLIT(val, val)
72 #define LDK(x) x
73
LDA(const R * x,INT ivs,const R * aligned_like)74 static inline V LDA(const R *x, INT ivs, const R *aligned_like)
75 {
76 (void)aligned_like; /* UNUSED */
77 (void)ivs; /* UNUSED */
78 return SUFF(_mm256_loadu_p)(x);
79 }
80
STA(R * x,V v,INT ovs,const R * aligned_like)81 static inline void STA(R *x, V v, INT ovs, const R *aligned_like)
82 {
83 (void)aligned_like; /* UNUSED */
84 (void)ovs; /* UNUSED */
85 SUFF(_mm256_storeu_p)(x, v);
86 }
87
88 #if FFTW_SINGLE
89
90 # ifdef _MSC_VER
91 /* Temporarily disable the warning "uninitialized local variable
92 'name' used" and runtime checks for using a variable before it is
93 defined which is erroneously triggered by the LOADL0 / LOADH macros
94 as they only modify VAL partly each. */
95 # ifndef __INTEL_COMPILER
96 # pragma warning(disable : 4700)
97 # pragma runtime_checks("u", off)
98 # endif
99 # endif
100 # ifdef __INTEL_COMPILER
101 # pragma warning(disable : 592)
102 # endif
103
104 #define LOADH(addr, val) _mm_loadh_pi(val, (const __m64 *)(addr))
105 #define LOADL(addr, val) _mm_loadl_pi(val, (const __m64 *)(addr))
106 #define STOREH(addr, val) _mm_storeh_pi((__m64 *)(addr), val)
107 #define STOREL(addr, val) _mm_storel_pi((__m64 *)(addr), val)
108
LD(const R * x,INT ivs,const R * aligned_like)109 static inline V LD(const R *x, INT ivs, const R *aligned_like)
110 {
111 __m128 l0, l1, h0, h1;
112 (void)aligned_like; /* UNUSED */
113 #if defined(__ICC) || (__GNUC__ > 4) || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)
114 l0 = LOADL(x, SUFF(_mm_undefined_p)());
115 l1 = LOADL(x + ivs, SUFF(_mm_undefined_p)());
116 h0 = LOADL(x + 2*ivs, SUFF(_mm_undefined_p)());
117 h1 = LOADL(x + 3*ivs, SUFF(_mm_undefined_p)());
118 #else
119 l0 = LOADL(x, l0);
120 l1 = LOADL(x + ivs, l1);
121 h0 = LOADL(x + 2*ivs, h0);
122 h1 = LOADL(x + 3*ivs, h1);
123 #endif
124 l0 = SUFF(_mm_movelh_p)(l0,l1);
125 h0 = SUFF(_mm_movelh_p)(h0,h1);
126 return _mm256_insertf128_ps(_mm256_castps128_ps256(l0), h0, 1);
127 }
128
129 # ifdef _MSC_VER
130 # ifndef __INTEL_COMPILER
131 # pragma warning(default : 4700)
132 # pragma runtime_checks("u", restore)
133 # endif
134 # endif
135 # ifdef __INTEL_COMPILER
136 # pragma warning(default : 592)
137 # endif
138
ST(R * x,V v,INT ovs,const R * aligned_like)139 static inline void ST(R *x, V v, INT ovs, const R *aligned_like)
140 {
141 __m128 h = _mm256_extractf128_ps(v, 1);
142 __m128 l = _mm256_castps256_ps128(v);
143 (void)aligned_like; /* UNUSED */
144 /* WARNING: the extra_iter hack depends upon STOREL occurring
145 after STOREH */
146 STOREH(x + 3*ovs, h);
147 STOREL(x + 2*ovs, h);
148 STOREH(x + ovs, l);
149 STOREL(x, l);
150 }
151
152 #define STM2(x, v, ovs, aligned_like) /* no-op */
STN2(R * x,V v0,V v1,INT ovs)153 static inline void STN2(R *x, V v0, V v1, INT ovs)
154 {
155 V x0 = VSHUF(v0, v1, SHUFVALS(0, 1, 0, 1));
156 V x1 = VSHUF(v0, v1, SHUFVALS(2, 3, 2, 3));
157 __m128 h0 = _mm256_extractf128_ps(x0, 1);
158 __m128 l0 = _mm256_castps256_ps128(x0);
159 __m128 h1 = _mm256_extractf128_ps(x1, 1);
160 __m128 l1 = _mm256_castps256_ps128(x1);
161 *(__m128 *)(x + 3*ovs) = h1;
162 *(__m128 *)(x + 2*ovs) = h0;
163 *(__m128 *)(x + 1*ovs) = l1;
164 *(__m128 *)(x + 0*ovs) = l0;
165 }
166
167 #define STM4(x, v, ovs, aligned_like) /* no-op */
168 #define STN4(x, v0, v1, v2, v3, ovs) \
169 { \
170 V xxx0, xxx1, xxx2, xxx3; \
171 V yyy0, yyy1, yyy2, yyy3; \
172 xxx0 = _mm256_unpacklo_ps(v0, v2); \
173 xxx1 = _mm256_unpackhi_ps(v0, v2); \
174 xxx2 = _mm256_unpacklo_ps(v1, v3); \
175 xxx3 = _mm256_unpackhi_ps(v1, v3); \
176 yyy0 = _mm256_unpacklo_ps(xxx0, xxx2); \
177 yyy1 = _mm256_unpackhi_ps(xxx0, xxx2); \
178 yyy2 = _mm256_unpacklo_ps(xxx1, xxx3); \
179 yyy3 = _mm256_unpackhi_ps(xxx1, xxx3); \
180 *(__m128 *)(x + 0 * ovs) = _mm256_castps256_ps128(yyy0); \
181 *(__m128 *)(x + 4 * ovs) = _mm256_extractf128_ps(yyy0, 1); \
182 *(__m128 *)(x + 1 * ovs) = _mm256_castps256_ps128(yyy1); \
183 *(__m128 *)(x + 5 * ovs) = _mm256_extractf128_ps(yyy1, 1); \
184 *(__m128 *)(x + 2 * ovs) = _mm256_castps256_ps128(yyy2); \
185 *(__m128 *)(x + 6 * ovs) = _mm256_extractf128_ps(yyy2, 1); \
186 *(__m128 *)(x + 3 * ovs) = _mm256_castps256_ps128(yyy3); \
187 *(__m128 *)(x + 7 * ovs) = _mm256_extractf128_ps(yyy3, 1); \
188 }
189
190 #else
VMOVAPD_LD(const R * x)191 static inline __m128d VMOVAPD_LD(const R *x)
192 {
193 /* gcc-4.6 miscompiles the combination _mm256_castpd128_pd256(VMOVAPD_LD(x))
194 into a 256-bit vmovapd, which requires 32-byte aligment instead of
195 16-byte alignment.
196
197 Force the use of vmovapd via asm until compilers stabilize.
198 */
199 #if defined(__GNUC__)
200 __m128d var;
201 __asm__("vmovapd %1, %0\n" : "=x"(var) : "m"(x[0]));
202 return var;
203 #else
204 return *(const __m128d *)x;
205 #endif
206 }
207
LD(const R * x,INT ivs,const R * aligned_like)208 static inline V LD(const R *x, INT ivs, const R *aligned_like)
209 {
210 V var;
211 (void)aligned_like; /* UNUSED */
212 var = _mm256_castpd128_pd256(VMOVAPD_LD(x));
213 var = _mm256_insertf128_pd(var, *(const __m128d *)(x+ivs), 1);
214 return var;
215 }
216
ST(R * x,V v,INT ovs,const R * aligned_like)217 static inline void ST(R *x, V v, INT ovs, const R *aligned_like)
218 {
219 (void)aligned_like; /* UNUSED */
220 /* WARNING: the extra_iter hack depends upon the store of the low
221 part occurring after the store of the high part */
222 *(__m128d *)(x + ovs) = _mm256_extractf128_pd(v, 1);
223 *(__m128d *)x = _mm256_castpd256_pd128(v);
224 }
225
226
227 #define STM2 ST
228 #define STN2(x, v0, v1, ovs) /* nop */
229 #define STM4(x, v, ovs, aligned_like) /* no-op */
230
231 /* STN4 is a macro, not a function, thanks to Visual C++ developers
232 deciding "it would be infrequent that people would want to pass more
233 than 3 [__m128 parameters] by value." Even though the comment
234 was made about __m128 parameters, it appears to apply to __m256
235 parameters as well. */
236 #define STN4(x, v0, v1, v2, v3, ovs) \
237 { \
238 V xxx0, xxx1, xxx2, xxx3; \
239 xxx0 = _mm256_unpacklo_pd(v0, v1); \
240 xxx1 = _mm256_unpackhi_pd(v0, v1); \
241 xxx2 = _mm256_unpacklo_pd(v2, v3); \
242 xxx3 = _mm256_unpackhi_pd(v2, v3); \
243 STA(x, _mm256_permute2f128_pd(xxx0, xxx2, 0x20), 0, 0); \
244 STA(x + ovs, _mm256_permute2f128_pd(xxx1, xxx3, 0x20), 0, 0); \
245 STA(x + 2 * ovs, _mm256_permute2f128_pd(xxx0, xxx2, 0x31), 0, 0); \
246 STA(x + 3 * ovs, _mm256_permute2f128_pd(xxx1, xxx3, 0x31), 0, 0); \
247 }
248 #endif
249
FLIP_RI(V x)250 static inline V FLIP_RI(V x)
251 {
252 return VPERM1(x, DS(SHUFVALD(1, 0), SHUFVALS(1, 0, 3, 2)));
253 }
254
VCONJ(V x)255 static inline V VCONJ(V x)
256 {
257 /* Produce a SIMD vector[VL] of (0 + -0i).
258
259 We really want to write this:
260
261 V pmpm = VLIT(-0.0, 0.0);
262
263 but historically some compilers have ignored the distiction
264 between +0 and -0. It looks like 'gcc-8 -fast-math' treats -0
265 as 0 too.
266 */
267 union uvec {
268 unsigned u[8];
269 V v;
270 };
271 static const union uvec pmpm = {
272 #ifdef FFTW_SINGLE
273 { 0x00000000, 0x80000000, 0x00000000, 0x80000000,
274 0x00000000, 0x80000000, 0x00000000, 0x80000000 }
275 #else
276 { 0x00000000, 0x00000000, 0x00000000, 0x80000000,
277 0x00000000, 0x00000000, 0x00000000, 0x80000000 }
278 #endif
279 };
280 return VXOR(pmpm.v, x);
281 }
282
VBYI(V x)283 static inline V VBYI(V x)
284 {
285 return FLIP_RI(VCONJ(x));
286 }
287
288 /* FMA support */
289 #define VFMA SUFF(_mm256_fmadd_p)
290 #define VFNMS SUFF(_mm256_fnmadd_p)
291 #define VFMS SUFF(_mm256_fmsub_p)
292 #define VFMAI(b, c) SUFF(_mm256_addsub_p)(c, FLIP_RI(b)) /* VADD(c, VBYI(b)) */
293 #define VFNMSI(b, c) VSUB(c, VBYI(b))
294 #define VFMACONJ(b,c) VADD(VCONJ(b),c)
295 #define VFMSCONJ(b,c) VSUB(VCONJ(b),c)
296 #define VFNMSCONJ(b,c) SUFF(_mm256_addsub_p)(c, b) /* VSUB(c, VCONJ(b)) */
297
VZMUL(V tx,V sr)298 static inline V VZMUL(V tx, V sr)
299 {
300 /* V tr = VDUPL(tx); */
301 /* V ti = VDUPH(tx); */
302 /* tr = VMUL(sr, tr); */
303 /* sr = VBYI(sr); */
304 /* return VFMA(ti, sr, tr); */
305 return SUFF(_mm256_fmaddsub_p)(sr, VDUPL(tx), VMUL(FLIP_RI(sr), VDUPH(tx)));
306 }
307
VZMULJ(V tx,V sr)308 static inline V VZMULJ(V tx, V sr)
309 {
310 /* V tr = VDUPL(tx); */
311 /* V ti = VDUPH(tx); */
312 /* tr = VMUL(sr, tr); */
313 /* sr = VBYI(sr); */
314 /* return VFNMS(ti, sr, tr); */
315 return SUFF(_mm256_fmsubadd_p)(sr, VDUPL(tx), VMUL(FLIP_RI(sr), VDUPH(tx)));
316 }
317
VZMULI(V tx,V sr)318 static inline V VZMULI(V tx, V sr)
319 {
320 V tr = VDUPL(tx);
321 V ti = VDUPH(tx);
322 ti = VMUL(ti, sr);
323 sr = VBYI(sr);
324 return VFMS(tr, sr, ti);
325 /*
326 * Keep the old version
327 * (2 permute, 1 shuffle, 1 constant load (L1), 1 xor, 2 fp), since the below FMA one
328 * would be 2 permute, 1 shuffle, 1 xor (setzero), 3 fp), but with a longer pipeline.
329 *
330 * Alternative new fma version:
331 * return SUFF(_mm256_addsub_p)(SUFF(_mm256_fnmadd_p)(sr, VDUPH(tx), SUFF(_mm256_setzero_p)()),
332 * VMUL(FLIP_RI(sr), VDUPL(tx)));
333 */
334 }
335
VZMULIJ(V tx,V sr)336 static inline V VZMULIJ(V tx, V sr)
337 {
338 /* V tr = VDUPL(tx); */
339 /* V ti = VDUPH(tx); */
340 /* ti = VMUL(ti, sr); */
341 /* sr = VBYI(sr); */
342 /* return VFMA(tr, sr, ti); */
343 return SUFF(_mm256_fmaddsub_p)(sr, VDUPH(tx), VMUL(FLIP_RI(sr), VDUPL(tx)));
344 }
345
346 /* twiddle storage #1: compact, slower */
347 #ifdef FFTW_SINGLE
348 # define VTW1(v,x) {TW_CEXP, v, x}, {TW_CEXP, v+1, x}, {TW_CEXP, v+2, x}, {TW_CEXP, v+3, x}
349 #else
350 # define VTW1(v,x) {TW_CEXP, v, x}, {TW_CEXP, v+1, x}
351 #endif
352 #define TWVL1 (VL)
353
BYTW1(const R * t,V sr)354 static inline V BYTW1(const R *t, V sr)
355 {
356 return VZMUL(LDA(t, 2, t), sr);
357 }
358
BYTWJ1(const R * t,V sr)359 static inline V BYTWJ1(const R *t, V sr)
360 {
361 return VZMULJ(LDA(t, 2, t), sr);
362 }
363
364 /* twiddle storage #2: twice the space, faster (when in cache) */
365 #ifdef FFTW_SINGLE
366 # define VTW2(v,x) \
367 {TW_COS, v, x}, {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+1, x}, \
368 {TW_COS, v+2, x}, {TW_COS, v+2, x}, {TW_COS, v+3, x}, {TW_COS, v+3, x}, \
369 {TW_SIN, v, -x}, {TW_SIN, v, x}, {TW_SIN, v+1, -x}, {TW_SIN, v+1, x}, \
370 {TW_SIN, v+2, -x}, {TW_SIN, v+2, x}, {TW_SIN, v+3, -x}, {TW_SIN, v+3, x}
371 #else
372 # define VTW2(v,x) \
373 {TW_COS, v, x}, {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+1, x}, \
374 {TW_SIN, v, -x}, {TW_SIN, v, x}, {TW_SIN, v+1, -x}, {TW_SIN, v+1, x}
375 #endif
376 #define TWVL2 (2 * VL)
377
BYTW2(const R * t,V sr)378 static inline V BYTW2(const R *t, V sr)
379 {
380 const V *twp = (const V *)t;
381 V si = FLIP_RI(sr);
382 V tr = twp[0], ti = twp[1];
383 return VFMA(tr, sr, VMUL(ti, si));
384 }
385
BYTWJ2(const R * t,V sr)386 static inline V BYTWJ2(const R *t, V sr)
387 {
388 const V *twp = (const V *)t;
389 V si = FLIP_RI(sr);
390 V tr = twp[0], ti = twp[1];
391 return VFNMS(ti, si, VMUL(tr, sr));
392 }
393
394 /* twiddle storage #3 */
395 #define VTW3 VTW1
396 #define TWVL3 TWVL1
397
398 /* twiddle storage for split arrays */
399 #ifdef FFTW_SINGLE
400 # define VTWS(v,x) \
401 {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+2, x}, {TW_COS, v+3, x}, \
402 {TW_COS, v+4, x}, {TW_COS, v+5, x}, {TW_COS, v+6, x}, {TW_COS, v+7, x}, \
403 {TW_SIN, v, x}, {TW_SIN, v+1, x}, {TW_SIN, v+2, x}, {TW_SIN, v+3, x}, \
404 {TW_SIN, v+4, x}, {TW_SIN, v+5, x}, {TW_SIN, v+6, x}, {TW_SIN, v+7, x}
405 #else
406 # define VTWS(v,x) \
407 {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+2, x}, {TW_COS, v+3, x}, \
408 {TW_SIN, v, x}, {TW_SIN, v+1, x}, {TW_SIN, v+2, x}, {TW_SIN, v+3, x}
409 #endif
410 #define TWVLS (2 * VL)
411
412 #define VLEAVE _mm256_zeroupper
413
414 #include "simd-common.h"
415