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