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