poulpy/spqlios/lib/spqlios/cplx/cplx_fftvec_avx2_fma.c
spqlios/lib/spqlios/cplx/cplx_fftvec_avx2_fma.cRaw
1#include <immintrin.h>
2#include <string.h>
3
4#include "cplx_fft_internal.h"
5#include "cplx_fft_private.h"
6
7typedef double D4MEM[4];
8
9EXPORT void cplx_fftvec_mul_fma(const CPLX_FFTVEC_MUL_PRECOMP* precomp, void* r, const void* a, const void* b) {
10 const uint32_t m = precomp->m;
11 const double(*aa)[4] = (double(*)[4])a;
12 const double(*bb)[4] = (double(*)[4])b;
13 double(*rr)[4] = (double(*)[4])r;
14 const double(*const aend)[4] = aa + (m >> 1);
15 do {
16 /*
17BEGIN_TEMPLATE
18const __m256d ari% = _mm256_loadu_pd(aa[%]);
19const __m256d bri% = _mm256_loadu_pd(bb[%]);
20const __m256d bir% = _mm256_shuffle_pd(bri%,bri%, 5); // conj of b
21const __m256d aii% = _mm256_shuffle_pd(ari%,ari%, 15); // im of a
22const __m256d pro% = _mm256_mul_pd(aii%,bir%);
23const __m256d arr% = _mm256_shuffle_pd(ari%,ari%, 0); // rr of a
24const __m256d res% = _mm256_fmaddsub_pd(arr%,bri%,pro%);
25_mm256_storeu_pd(rr[%],res%);
26rr += @; // ONCE
27aa += @; // ONCE
28bb += @; // ONCE
29END_TEMPLATE
30 */
31 // BEGIN_INTERLEAVE 4
32 // This block is automatically generated from the template above
33 // by the interleave.pl script. Please do not edit by hand
34 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
35 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
36 const __m256d ari2 = _mm256_loadu_pd(aa[2]);
37 const __m256d ari3 = _mm256_loadu_pd(aa[3]);
38 const __m256d bri0 = _mm256_loadu_pd(bb[0]);
39 const __m256d bri1 = _mm256_loadu_pd(bb[1]);
40 const __m256d bri2 = _mm256_loadu_pd(bb[2]);
41 const __m256d bri3 = _mm256_loadu_pd(bb[3]);
42 const __m256d bir0 = _mm256_shuffle_pd(bri0, bri0, 5); // conj of b
43 const __m256d bir1 = _mm256_shuffle_pd(bri1, bri1, 5); // conj of b
44 const __m256d bir2 = _mm256_shuffle_pd(bri2, bri2, 5); // conj of b
45 const __m256d bir3 = _mm256_shuffle_pd(bri3, bri3, 5); // conj of b
46 const __m256d aii0 = _mm256_shuffle_pd(ari0, ari0, 15); // im of a
47 const __m256d aii1 = _mm256_shuffle_pd(ari1, ari1, 15); // im of a
48 const __m256d aii2 = _mm256_shuffle_pd(ari2, ari2, 15); // im of a
49 const __m256d aii3 = _mm256_shuffle_pd(ari3, ari3, 15); // im of a
50 const __m256d pro0 = _mm256_mul_pd(aii0, bir0);
51 const __m256d pro1 = _mm256_mul_pd(aii1, bir1);
52 const __m256d pro2 = _mm256_mul_pd(aii2, bir2);
53 const __m256d pro3 = _mm256_mul_pd(aii3, bir3);
54 const __m256d arr0 = _mm256_shuffle_pd(ari0, ari0, 0); // rr of a
55 const __m256d arr1 = _mm256_shuffle_pd(ari1, ari1, 0); // rr of a
56 const __m256d arr2 = _mm256_shuffle_pd(ari2, ari2, 0); // rr of a
57 const __m256d arr3 = _mm256_shuffle_pd(ari3, ari3, 0); // rr of a
58 const __m256d res0 = _mm256_fmaddsub_pd(arr0, bri0, pro0);
59 const __m256d res1 = _mm256_fmaddsub_pd(arr1, bri1, pro1);
60 const __m256d res2 = _mm256_fmaddsub_pd(arr2, bri2, pro2);
61 const __m256d res3 = _mm256_fmaddsub_pd(arr3, bri3, pro3);
62 _mm256_storeu_pd(rr[0], res0);
63 _mm256_storeu_pd(rr[1], res1);
64 _mm256_storeu_pd(rr[2], res2);
65 _mm256_storeu_pd(rr[3], res3);
66 rr += 4; // ONCE
67 aa += 4; // ONCE
68 bb += 4; // ONCE
69 // END_INTERLEAVE
70 } while (aa < aend);
71}
72
73
74
75EXPORT void cplx_fftvec_addmul_fma(const CPLX_FFTVEC_ADDMUL_PRECOMP* precomp, void* r, const void* a, const void* b) {
76 const uint32_t m = precomp->m;
77 const double(*aa)[4] = (double(*)[4])a;
78 const double(*bb)[4] = (double(*)[4])b;
79 double(*rr)[4] = (double(*)[4])r;
80 const double(*const aend)[4] = aa + (m >> 1);
81 do {
82 /*
83BEGIN_TEMPLATE
84const __m256d ari% = _mm256_loadu_pd(aa[%]);
85const __m256d bri% = _mm256_loadu_pd(bb[%]);
86const __m256d rri% = _mm256_loadu_pd(rr[%]);
87const __m256d bir% = _mm256_shuffle_pd(bri%,bri%, 5);
88const __m256d aii% = _mm256_shuffle_pd(ari%,ari%, 15);
89const __m256d pro% = _mm256_fmaddsub_pd(aii%,bir%,rri%);
90const __m256d arr% = _mm256_shuffle_pd(ari%,ari%, 0);
91const __m256d res% = _mm256_fmaddsub_pd(arr%,bri%,pro%);
92_mm256_storeu_pd(rr[%],res%);
93rr += @; // ONCE
94aa += @; // ONCE
95bb += @; // ONCE
96END_TEMPLATE
97 */
98 // BEGIN_INTERLEAVE 2
99 // This block is automatically generated from the template above
100 // by the interleave.pl script. Please do not edit by hand
101 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
102 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
103 const __m256d bri0 = _mm256_loadu_pd(bb[0]);
104 const __m256d bri1 = _mm256_loadu_pd(bb[1]);
105 const __m256d rri0 = _mm256_loadu_pd(rr[0]);
106 const __m256d rri1 = _mm256_loadu_pd(rr[1]);
107 const __m256d bir0 = _mm256_shuffle_pd(bri0, bri0, 5);
108 const __m256d bir1 = _mm256_shuffle_pd(bri1, bri1, 5);
109 const __m256d aii0 = _mm256_shuffle_pd(ari0, ari0, 15);
110 const __m256d aii1 = _mm256_shuffle_pd(ari1, ari1, 15);
111 const __m256d pro0 = _mm256_fmaddsub_pd(aii0, bir0, rri0);
112 const __m256d pro1 = _mm256_fmaddsub_pd(aii1, bir1, rri1);
113 const __m256d arr0 = _mm256_shuffle_pd(ari0, ari0, 0);
114 const __m256d arr1 = _mm256_shuffle_pd(ari1, ari1, 0);
115 const __m256d res0 = _mm256_fmaddsub_pd(arr0, bri0, pro0);
116 const __m256d res1 = _mm256_fmaddsub_pd(arr1, bri1, pro1);
117 _mm256_storeu_pd(rr[0], res0);
118 _mm256_storeu_pd(rr[1], res1);
119 rr += 2; // ONCE
120 aa += 2; // ONCE
121 bb += 2; // ONCE
122 // END_INTERLEAVE
123 } while (aa < aend);
124}
125
126EXPORT void cplx_fftvec_bitwiddle_fma(const CPLX_FFTVEC_BITWIDDLE_PRECOMP* precomp, void* a, uint64_t slicea,
127 const void* omg) {
128 const uint32_t m = precomp->m;
129 const uint64_t OFFSET = slicea / sizeof(D4MEM);
130 D4MEM* aa = (D4MEM*)a;
131 const double(*const aend)[4] = aa + (m >> 1);
132 const __m256d om = _mm256_loadu_pd(omg);
133 const __m256d om1rr = _mm256_shuffle_pd(om, om, 0);
134 const __m256d om1ii = _mm256_shuffle_pd(om, om, 15);
135 const __m256d om2rr = _mm256_shuffle_pd(om, om, 0);
136 const __m256d om2ii = _mm256_shuffle_pd(om, om, 0);
137 const __m256d om3rr = _mm256_shuffle_pd(om, om, 15);
138 const __m256d om3ii = _mm256_shuffle_pd(om, om, 15);
139 do {
140 /*
141BEGIN_TEMPLATE
142__m256d ari% = _mm256_loadu_pd(aa[%]);
143__m256d bri% = _mm256_loadu_pd((aa+OFFSET)[%]);
144__m256d cri% = _mm256_loadu_pd((aa+2*OFFSET)[%]);
145__m256d dri% = _mm256_loadu_pd((aa+3*OFFSET)[%]);
146__m256d pa% = _mm256_shuffle_pd(cri%,cri%,5);
147__m256d pb% = _mm256_shuffle_pd(dri%,dri%,5);
148pa% = _mm256_mul_pd(pa%,om1ii);
149pb% = _mm256_mul_pd(pb%,om1ii);
150pa% = _mm256_fmaddsub_pd(cri%,om1rr,pa%);
151pb% = _mm256_fmaddsub_pd(dri%,om1rr,pb%);
152cri% = _mm256_sub_pd(ari%,pa%);
153dri% = _mm256_sub_pd(bri%,pb%);
154ari% = _mm256_add_pd(ari%,pa%);
155bri% = _mm256_add_pd(bri%,pb%);
156pa% = _mm256_shuffle_pd(bri%,bri%,5);
157pb% = _mm256_shuffle_pd(dri%,dri%,5);
158pa% = _mm256_mul_pd(pa%,om2ii);
159pb% = _mm256_mul_pd(pb%,om3ii);
160pa% = _mm256_fmaddsub_pd(bri%,om2rr,pa%);
161pb% = _mm256_fmaddsub_pd(dri%,om3rr,pb%);
162bri% = _mm256_sub_pd(ari%,pa%);
163dri% = _mm256_sub_pd(cri%,pb%);
164ari% = _mm256_add_pd(ari%,pa%);
165cri% = _mm256_add_pd(cri%,pb%);
166_mm256_storeu_pd(aa[%], ari%);
167_mm256_storeu_pd((aa+OFFSET)[%],bri%);
168_mm256_storeu_pd((aa+2*OFFSET)[%],cri%);
169_mm256_storeu_pd((aa+3*OFFSET)[%],dri%);
170aa += @; // ONCE
171END_TEMPLATE
172 */
173 // BEGIN_INTERLEAVE 1
174 // This block is automatically generated from the template above
175 // by the interleave.pl script. Please do not edit by hand
176 __m256d ari0 = _mm256_loadu_pd(aa[0]);
177 __m256d bri0 = _mm256_loadu_pd((aa + OFFSET)[0]);
178 __m256d cri0 = _mm256_loadu_pd((aa + 2 * OFFSET)[0]);
179 __m256d dri0 = _mm256_loadu_pd((aa + 3 * OFFSET)[0]);
180 __m256d pa0 = _mm256_shuffle_pd(cri0, cri0, 5);
181 __m256d pb0 = _mm256_shuffle_pd(dri0, dri0, 5);
182 pa0 = _mm256_mul_pd(pa0, om1ii);
183 pb0 = _mm256_mul_pd(pb0, om1ii);
184 pa0 = _mm256_fmaddsub_pd(cri0, om1rr, pa0);
185 pb0 = _mm256_fmaddsub_pd(dri0, om1rr, pb0);
186 cri0 = _mm256_sub_pd(ari0, pa0);
187 dri0 = _mm256_sub_pd(bri0, pb0);
188 ari0 = _mm256_add_pd(ari0, pa0);
189 bri0 = _mm256_add_pd(bri0, pb0);
190 pa0 = _mm256_shuffle_pd(bri0, bri0, 5);
191 pb0 = _mm256_shuffle_pd(dri0, dri0, 5);
192 pa0 = _mm256_mul_pd(pa0, om2ii);
193 pb0 = _mm256_mul_pd(pb0, om3ii);
194 pa0 = _mm256_fmaddsub_pd(bri0, om2rr, pa0);
195 pb0 = _mm256_fmaddsub_pd(dri0, om3rr, pb0);
196 bri0 = _mm256_sub_pd(ari0, pa0);
197 dri0 = _mm256_sub_pd(cri0, pb0);
198 ari0 = _mm256_add_pd(ari0, pa0);
199 cri0 = _mm256_add_pd(cri0, pb0);
200 _mm256_storeu_pd(aa[0], ari0);
201 _mm256_storeu_pd((aa + OFFSET)[0], bri0);
202 _mm256_storeu_pd((aa + 2 * OFFSET)[0], cri0);
203 _mm256_storeu_pd((aa + 3 * OFFSET)[0], dri0);
204 aa += 1; // ONCE
205 // END_INTERLEAVE
206 } while (aa < aend);
207}
208
209EXPORT void cplx_fftvec_sub2_to_fma(uint32_t m, void* r, const void* a, const void* b) {
210 const double(*aa)[4] = (double(*)[4])a;
211 const double(*bb)[4] = (double(*)[4])b;
212 double(*rr)[4] = (double(*)[4])r;
213 const double(*const aend)[4] = aa + (m >> 1);
214 do {
215 /*
216BEGIN_TEMPLATE
217const __m256d ari% = _mm256_loadu_pd(aa[%]);
218const __m256d bri% = _mm256_loadu_pd(bb[%]);
219const __m256d sum% = _mm256_add_pd(ari%,bri%);
220const __m256d rri% = _mm256_loadu_pd(rr[%]);
221const __m256d res% = _mm256_sub_pd(rri%,sum%);
222_mm256_storeu_pd(rr[%],res%);
223END_TEMPLATE
224 */
225 // BEGIN_INTERLEAVE 4
226 // This block is automatically generated from the template above
227 // by the interleave.pl script. Please do not edit by hand
228 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
229 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
230 const __m256d ari2 = _mm256_loadu_pd(aa[2]);
231 const __m256d ari3 = _mm256_loadu_pd(aa[3]);
232 const __m256d bri0 = _mm256_loadu_pd(bb[0]);
233 const __m256d bri1 = _mm256_loadu_pd(bb[1]);
234 const __m256d bri2 = _mm256_loadu_pd(bb[2]);
235 const __m256d bri3 = _mm256_loadu_pd(bb[3]);
236 const __m256d sum0 = _mm256_add_pd(ari0, bri0);
237 const __m256d sum1 = _mm256_add_pd(ari1, bri1);
238 const __m256d sum2 = _mm256_add_pd(ari2, bri2);
239 const __m256d sum3 = _mm256_add_pd(ari3, bri3);
240 const __m256d rri0 = _mm256_loadu_pd(rr[0]);
241 const __m256d rri1 = _mm256_loadu_pd(rr[1]);
242 const __m256d rri2 = _mm256_loadu_pd(rr[2]);
243 const __m256d rri3 = _mm256_loadu_pd(rr[3]);
244 const __m256d res0 = _mm256_sub_pd(rri0, sum0);
245 const __m256d res1 = _mm256_sub_pd(rri1, sum1);
246 const __m256d res2 = _mm256_sub_pd(rri2, sum2);
247 const __m256d res3 = _mm256_sub_pd(rri3, sum3);
248 _mm256_storeu_pd(rr[0], res0);
249 _mm256_storeu_pd(rr[1], res1);
250 _mm256_storeu_pd(rr[2], res2);
251 _mm256_storeu_pd(rr[3], res3);
252 // END_INTERLEAVE
253 rr += 4;
254 aa += 4;
255 bb += 4;
256 } while (aa < aend);
257}
258
259EXPORT void cplx_fftvec_add_fma(uint32_t m, void* r, const void* a, const void* b) {
260 const double(*aa)[4] = (double(*)[4])a;
261 const double(*bb)[4] = (double(*)[4])b;
262 double(*rr)[4] = (double(*)[4])r;
263 const double(*const aend)[4] = aa + (m >> 1);
264 do {
265 /*
266 BEGIN_TEMPLATE
267 const __m256d ari% = _mm256_loadu_pd(aa[%]);
268 const __m256d bri% = _mm256_loadu_pd(bb[%]);
269 const __m256d res% = _mm256_add_pd(ari%,bri%);
270 _mm256_storeu_pd(rr[%],res%);
271 rr += @; // ONCE
272 aa += @; // ONCE
273 bb += @; // ONCE
274 END_TEMPLATE
275 */
276 // BEGIN_INTERLEAVE 4
277 // This block is automatically generated from the template above
278 // by the interleave.pl script. Please do not edit by hand
279 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
280 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
281 const __m256d ari2 = _mm256_loadu_pd(aa[2]);
282 const __m256d ari3 = _mm256_loadu_pd(aa[3]);
283 const __m256d bri0 = _mm256_loadu_pd(bb[0]);
284 const __m256d bri1 = _mm256_loadu_pd(bb[1]);
285 const __m256d bri2 = _mm256_loadu_pd(bb[2]);
286 const __m256d bri3 = _mm256_loadu_pd(bb[3]);
287 const __m256d res0 = _mm256_add_pd(ari0, bri0);
288 const __m256d res1 = _mm256_add_pd(ari1, bri1);
289 const __m256d res2 = _mm256_add_pd(ari2, bri2);
290 const __m256d res3 = _mm256_add_pd(ari3, bri3);
291 _mm256_storeu_pd(rr[0], res0);
292 _mm256_storeu_pd(rr[1], res1);
293 _mm256_storeu_pd(rr[2], res2);
294 _mm256_storeu_pd(rr[3], res3);
295 rr += 4; // ONCE
296 aa += 4; // ONCE
297 bb += 4; // ONCE
298 // END_INTERLEAVE
299 } while (aa < aend);
300}
301
302EXPORT void cplx_fftvec_twiddle_fma(const CPLX_FFTVEC_TWIDDLE_PRECOMP* precomp, void* a, void* b, const void* omg) {
303 const uint32_t m = precomp->m;
304 double(*aa)[4] = (double(*)[4])a;
305 double(*bb)[4] = (double(*)[4])b;
306 const double(*const aend)[4] = aa + (m >> 1);
307 const __m256d om = _mm256_loadu_pd(omg);
308 const __m256d omrr = _mm256_shuffle_pd(om, om, 0);
309 const __m256d omii = _mm256_shuffle_pd(om, om, 15);
310 do {
311 /*
312BEGIN_TEMPLATE
313const __m256d bri% = _mm256_loadu_pd(bb[%]);
314const __m256d bir% = _mm256_shuffle_pd(bri%,bri%,5);
315__m256d p% = _mm256_mul_pd(bir%,omii);
316p% = _mm256_fmaddsub_pd(bri%,omrr,p%);
317const __m256d ari% = _mm256_loadu_pd(aa[%]);
318_mm256_storeu_pd(aa[%],_mm256_add_pd(ari%,p%));
319_mm256_storeu_pd(bb[%],_mm256_sub_pd(ari%,p%));
320bb += @; // ONCE
321aa += @; // ONCE
322END_TEMPLATE
323 */
324 // BEGIN_INTERLEAVE 4
325 // This block is automatically generated from the template above
326 // by the interleave.pl script. Please do not edit by hand
327 const __m256d bri0 = _mm256_loadu_pd(bb[0]);
328 const __m256d bri1 = _mm256_loadu_pd(bb[1]);
329 const __m256d bri2 = _mm256_loadu_pd(bb[2]);
330 const __m256d bri3 = _mm256_loadu_pd(bb[3]);
331 const __m256d bir0 = _mm256_shuffle_pd(bri0, bri0, 5);
332 const __m256d bir1 = _mm256_shuffle_pd(bri1, bri1, 5);
333 const __m256d bir2 = _mm256_shuffle_pd(bri2, bri2, 5);
334 const __m256d bir3 = _mm256_shuffle_pd(bri3, bri3, 5);
335 __m256d p0 = _mm256_mul_pd(bir0, omii);
336 __m256d p1 = _mm256_mul_pd(bir1, omii);
337 __m256d p2 = _mm256_mul_pd(bir2, omii);
338 __m256d p3 = _mm256_mul_pd(bir3, omii);
339 p0 = _mm256_fmaddsub_pd(bri0, omrr, p0);
340 p1 = _mm256_fmaddsub_pd(bri1, omrr, p1);
341 p2 = _mm256_fmaddsub_pd(bri2, omrr, p2);
342 p3 = _mm256_fmaddsub_pd(bri3, omrr, p3);
343 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
344 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
345 const __m256d ari2 = _mm256_loadu_pd(aa[2]);
346 const __m256d ari3 = _mm256_loadu_pd(aa[3]);
347 _mm256_storeu_pd(aa[0], _mm256_add_pd(ari0, p0));
348 _mm256_storeu_pd(aa[1], _mm256_add_pd(ari1, p1));
349 _mm256_storeu_pd(aa[2], _mm256_add_pd(ari2, p2));
350 _mm256_storeu_pd(aa[3], _mm256_add_pd(ari3, p3));
351 _mm256_storeu_pd(bb[0], _mm256_sub_pd(ari0, p0));
352 _mm256_storeu_pd(bb[1], _mm256_sub_pd(ari1, p1));
353 _mm256_storeu_pd(bb[2], _mm256_sub_pd(ari2, p2));
354 _mm256_storeu_pd(bb[3], _mm256_sub_pd(ari3, p3));
355 bb += 4; // ONCE
356 aa += 4; // ONCE
357 // END_INTERLEAVE
358 } while (aa < aend);
359}
360
361EXPORT void cplx_fftvec_copy_fma(uint32_t m, void* r, const void* a) {
362 const double(*aa)[4] = (double(*)[4])a;
363 double(*rr)[4] = (double(*)[4])r;
364 const double(*const aend)[4] = aa + (m >> 1);
365 do {
366 /*
367BEGIN_TEMPLATE
368const __m256d ari% = _mm256_loadu_pd(aa[%]);
369_mm256_storeu_pd(rr[%],ari%);
370rr += @; // ONCE
371aa += @; // ONCE
372END_TEMPLATE
373 */
374 // BEGIN_INTERLEAVE 4
375 // This block is automatically generated from the template above
376 // by the interleave.pl script. Please do not edit by hand
377 const __m256d ari0 = _mm256_loadu_pd(aa[0]);
378 const __m256d ari1 = _mm256_loadu_pd(aa[1]);
379 const __m256d ari2 = _mm256_loadu_pd(aa[2]);
380 const __m256d ari3 = _mm256_loadu_pd(aa[3]);
381 _mm256_storeu_pd(rr[0], ari0);
382 _mm256_storeu_pd(rr[1], ari1);
383 _mm256_storeu_pd(rr[2], ari2);
384 _mm256_storeu_pd(rr[3], ari3);
385 rr += 4; // ONCE
386 aa += 4; // ONCE
387 // END_INTERLEAVE
388 } while (aa < aend);
389}