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