elpa2_kernels_complex_avx-avx2_2hv_double_precision.c 50.5 KB
Newer Older
1
2
//    This file is part of ELPA.
//
Andreas Marek's avatar
Andreas Marek committed
3
//    The ELPA library was originally created by the ELPA consortium,
4
5
//    consisting of the following organizations:
//
6
7
//    - Max Planck Computing and Data Facility (MPCDF), formerly known as
//      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
8
9
10
//    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
//      Informatik,
//    - Technische Universität München, Lehrstuhl für Informatik mit
Andreas Marek's avatar
Andreas Marek committed
11
12
13
14
15
//      Schwerpunkt Wissenschaftliches Rechnen ,
//    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
//    - Max-Plack-Institut für Mathematik in den Naturwissenschaftrn,
//      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
//      and
16
17
//    - IBM Deutschland GmbH
//
18
//    This particular source code file contains additions, changes and
Andreas Marek's avatar
Andreas Marek committed
19
//    enhancements authored by Intel Corporation which is not part of
20
//    the ELPA consortium.
21
22
//
//    More information can be found here:
23
//    http://elpa.mpcdf.mpg.de/
24
25
//
//    ELPA is free software: you can redistribute it and/or modify
Andreas Marek's avatar
Andreas Marek committed
26
27
//    it under the terms of the version 3 of the license of the
//    GNU Lesser General Public License as published by the Free
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
//    Software Foundation.
//
//    ELPA is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU Lesser General Public License for more details.
//
//    You should have received a copy of the GNU Lesser General Public License
//    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
//
//    ELPA reflects a substantial effort on the part of the original
//    ELPA consortium, and we ask you to respect the spirit of the
//    license that we chose: i.e., please contribute any changes you
//    may have back to the original ELPA library distribution, and keep
//    any derivatives of ELPA under the same license that we chose for
//    the original distribution, the GNU Lesser General Public License.
//
//
// --------------------------------------------------------------------------------------------------
//
// This file contains the compute intensive kernels for the Householder transformations.
// It should be compiled with the highest possible optimization level.
//
// On Intel Nehalem or Intel Westmere or AMD Magny Cours use -O3 -msse3
// On Intel Sandy Bridge use -O3 -mavx
//
// Copyright of the original code rests with the authors inside the ELPA
// consortium. The copyright of any additional modifications shall rest
// with their original authors, but shall adhere to the licensing terms
// distributed along with the original code in the file "COPYING".
//
// Author: Alexander Heinecke (alexander.heinecke@mytum.de)
60
// Adapted for building a shared-library by Andreas Marek, MPCDF (andreas.marek@mpcdf.mpg.de)
61
// --------------------------------------------------------------------------------------------------
62
#include "config-f90.h"
63

Andreas Marek's avatar
Andreas Marek committed
64
#include <complex.h>
65
66
67
68
#include <x86intrin.h>

#define __forceinline __attribute__((always_inline))

69
70
#ifdef HAVE_AVX2

71
72
73
74
75
76
77
78
79
80
81
#ifdef __FMA4__
#define __ELPA_USE_FMA__
#define _mm256_FMADDSUB_pd(a,b,c) _mm256_maddsub_pd(a,b,c)
#define _mm256_FMSUBADD_pd(a,b,c) _mm256_msubadd_pd(a,b,c)
#endif

#ifdef __AVX2__
#define __ELPA_USE_FMA__
#define _mm256_FMADDSUB_pd(a,b,c) _mm256_fmaddsub_pd(a,b,c)
#define _mm256_FMSUBADD_pd(a,b,c) _mm256_fmsubadd_pd(a,b,c)
#endif
82

83
84
#endif

85
//Forward declaration
86
87
88
89
static __forceinline void hh_trafo_complex_kernel_8_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s);
static __forceinline void hh_trafo_complex_kernel_6_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s);
static __forceinline void hh_trafo_complex_kernel_4_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s);
static __forceinline void hh_trafo_complex_kernel_2_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s);
90

91
/*
92
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
93
!f> interface
94
95
!f>   subroutine double_hh_trafo_complex_avx_avx2_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                             bind(C, name="double_hh_trafo_complex_avx_avx2_2hv_double")
96
97
98
99
100
101
102
103
104
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int)     :: pnb, pnq, pldq, pldh
!f>     complex(kind=c_double)     :: q(*)
!f>     complex(kind=c_double)     :: hh(pnb,2)
!f>   end subroutine
!f> end interface
!f>#endif
*/

105
void double_hh_trafo_complex_avx_avx2_2hv_double(double complex* q, double complex* hh, int* pnb, int* pnq, int* pldq, int* pldh)
106
107
108
109
110
111
112
{
	int i;
	int nb = *pnb;
	int nq = *pldq;
	int ldq = *pldq;
	int ldh = *pldh;

Andreas Marek's avatar
Andreas Marek committed
113
	double complex s = conj(hh[(ldh)+1])*1.0;
114
115
116
117
118
119
120
121
	for (i = 2; i < nb; i++)
	{
		s += hh[i-1] * conj(hh[(i+ldh)]);
	}

#if 1
	for (i = 0; i < nq-4; i+=8)
	{
122
		hh_trafo_complex_kernel_8_AVX_2hv_double(&q[i], hh, nb, ldq, ldh, s);
123
124
125
	}
	if (nq-i > 0)
	{
126
		hh_trafo_complex_kernel_4_AVX_2hv_double(&q[i], hh, nb, ldq, ldh, s);
127
128
129
130
	}
#else
	for (i = 0; i < nq-4; i+=6)
	{
131
		hh_trafo_complex_kernel_6_AVX_2hv_double(&q[i], hh, nb, ldq, ldh, s);
132
133
134
	}
	if (nq-i > 2)
	{
135
		hh_trafo_complex_kernel_4_AVX_2hv_double(&q[i], hh, nb, ldq, ldh, s);
136
137
138
	}
	else if (nq-i > 0)
	{
139
		hh_trafo_complex_kernel_2_AVX_2hv_double(&q[i], hh, nb, ldq, ldh, s);
140
141
142
143
	}
#endif
}

144
static __forceinline void hh_trafo_complex_kernel_8_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s)
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
{
	double* q_dbl = (double*)q;
	double* hh_dbl = (double*)hh;
	double* s_dbl = (double*)(&s);

	__m256d x1, x2, x3, x4;
	__m256d y1, y2, y3, y4;
	__m256d q1, q2, q3, q4;
	__m256d h1_real, h1_imag, h2_real, h2_imag;
	__m256d tmp1, tmp2, tmp3, tmp4;
	int i=0;

	__m256d sign = (__m256d)_mm256_set_epi64x(0x8000000000000000, 0x8000000000000000, 0x8000000000000000, 0x8000000000000000);

	x1 = _mm256_load_pd(&q_dbl[(2*ldq)+0]);
	x2 = _mm256_load_pd(&q_dbl[(2*ldq)+4]);
	x3 = _mm256_load_pd(&q_dbl[(2*ldq)+8]);
	x4 = _mm256_load_pd(&q_dbl[(2*ldq)+12]);

	h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+1)*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+1)*2)+1]);
166
#ifndef __ELPA_USE_FMA__
167
168
169
170
171
172
173
174
175
176
	// conjugate
	h2_imag = _mm256_xor_pd(h2_imag, sign);
#endif

	y1 = _mm256_load_pd(&q_dbl[0]);
	y2 = _mm256_load_pd(&q_dbl[4]);
	y3 = _mm256_load_pd(&q_dbl[8]);
	y4 = _mm256_load_pd(&q_dbl[12]);

	tmp1 = _mm256_mul_pd(h2_imag, x1);
177
178
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_add_pd(y1, _mm256_FMSUBADD_pd(h2_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
179
180
181
182
#else
	y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, x2);
183
184
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_add_pd(y2, _mm256_FMSUBADD_pd(h2_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
185
186
187
188
#else
	y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, x3);
189
190
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_add_pd(y3, _mm256_FMSUBADD_pd(h2_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
191
192
193
194
#else
	y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	tmp4 = _mm256_mul_pd(h2_imag, x4);
195
196
#ifdef __ELPA_USE_FMA__
	y4 = _mm256_add_pd(y4, _mm256_FMSUBADD_pd(h2_real, x4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
197
198
199
200
201
202
203
204
205
206
207
208
209
#else
	y4 = _mm256_add_pd(y4, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

	for (i = 2; i < nb; i++)
	{
		q1 = _mm256_load_pd(&q_dbl[(2*i*ldq)+0]);
		q2 = _mm256_load_pd(&q_dbl[(2*i*ldq)+4]);
		q3 = _mm256_load_pd(&q_dbl[(2*i*ldq)+8]);
		q4 = _mm256_load_pd(&q_dbl[(2*i*ldq)+12]);

		h1_real = _mm256_broadcast_sd(&hh_dbl[(i-1)*2]);
		h1_imag = _mm256_broadcast_sd(&hh_dbl[((i-1)*2)+1]);
210
#ifndef __ELPA_USE_FMA__
211
212
213
214
215
		// conjugate
		h1_imag = _mm256_xor_pd(h1_imag, sign);
#endif

		tmp1 = _mm256_mul_pd(h1_imag, q1);
216
217
#ifdef __ELPA_USE_FMA__
		x1 = _mm256_add_pd(x1, _mm256_FMSUBADD_pd(h1_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
218
219
220
221
#else
		x1 = _mm256_add_pd(x1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h1_imag, q2);
222
223
#ifdef __ELPA_USE_FMA__
		x2 = _mm256_add_pd(x2, _mm256_FMSUBADD_pd(h1_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
224
225
226
227
#else
		x2 = _mm256_add_pd(x2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h1_imag, q3);
228
229
#ifdef __ELPA_USE_FMA__
		x3 = _mm256_add_pd(x3, _mm256_FMSUBADD_pd(h1_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
230
231
232
233
#else
		x3 = _mm256_add_pd(x3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
		tmp4 = _mm256_mul_pd(h1_imag, q4);
234
235
#ifdef __ELPA_USE_FMA__
		x4 = _mm256_add_pd(x4, _mm256_FMSUBADD_pd(h1_real, q4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
236
237
238
239
240
241
#else
		x4 = _mm256_add_pd(x4, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

		h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+i)*2]);
		h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+i)*2)+1]);
242
#ifndef __ELPA_USE_FMA__
243
244
245
246
247
		// conjugate
		h2_imag = _mm256_xor_pd(h2_imag, sign);
#endif

		tmp1 = _mm256_mul_pd(h2_imag, q1);
248
249
#ifdef __ELPA_USE_FMA__
		y1 = _mm256_add_pd(y1, _mm256_FMSUBADD_pd(h2_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
250
251
252
253
#else
		y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h2_imag, q2);
254
255
#ifdef __ELPA_USE_FMA__
		y2 = _mm256_add_pd(y2, _mm256_FMSUBADD_pd(h2_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
256
257
258
259
#else
		y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h2_imag, q3);
260
261
#ifdef __ELPA_USE_FMA__
		y3 = _mm256_add_pd(y3, _mm256_FMSUBADD_pd(h2_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
262
263
264
265
#else
		y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
		tmp4 = _mm256_mul_pd(h2_imag, q4);
266
267
#ifdef __ELPA_USE_FMA__
		y4 = _mm256_add_pd(y4, _mm256_FMSUBADD_pd(h2_real, q4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
268
269
270
271
272
273
274
#else
		y4 = _mm256_add_pd(y4, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif
	}

	h1_real = _mm256_broadcast_sd(&hh_dbl[(nb-1)*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[((nb-1)*2)+1]);
275
#ifndef __ELPA_USE_FMA__
276
277
278
279
280
281
282
283
284
285
	// conjugate
	h1_imag = _mm256_xor_pd(h1_imag, sign);
#endif

	q1 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+0]);
	q2 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+4]);
	q3 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+8]);
	q4 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+12]);

	tmp1 = _mm256_mul_pd(h1_imag, q1);
286
287
#ifdef __ELPA_USE_FMA__
	x1 = _mm256_add_pd(x1, _mm256_FMSUBADD_pd(h1_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
288
289
290
291
#else
	x1 = _mm256_add_pd(x1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, q2);
292
293
#ifdef __ELPA_USE_FMA__
	x2 = _mm256_add_pd(x2, _mm256_FMSUBADD_pd(h1_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
294
295
296
297
#else
	x2 = _mm256_add_pd(x2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, q3);
298
299
#ifdef __ELPA_USE_FMA__
	x3 = _mm256_add_pd(x3, _mm256_FMSUBADD_pd(h1_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
300
301
302
303
#else
	x3 = _mm256_add_pd(x3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	tmp4 = _mm256_mul_pd(h1_imag, q4);
304
305
#ifdef __ELPA_USE_FMA__
	x4 = _mm256_add_pd(x4, _mm256_FMSUBADD_pd(h1_real, q4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
306
307
308
309
310
311
312
313
314
315
#else
	x4 = _mm256_add_pd(x4, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

	h1_real = _mm256_broadcast_sd(&hh_dbl[0]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[1]);
	h1_real = _mm256_xor_pd(h1_real, sign);
	h1_imag = _mm256_xor_pd(h1_imag, sign);

	tmp1 = _mm256_mul_pd(h1_imag, x1);
316
317
#ifdef __ELPA_USE_FMA__
	x1 = _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
318
319
320
321
#else
	x1 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, x2);
322
323
#ifdef __ELPA_USE_FMA__
	x2 = _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5));
324
325
326
327
#else
	x2 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, x3);
328
329
#ifdef __ELPA_USE_FMA__
	x3 = _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5));
330
331
332
333
#else
	x3 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5));
#endif
	tmp4 = _mm256_mul_pd(h1_imag, x4);
334
335
#ifdef __ELPA_USE_FMA__
	x4 = _mm256_FMADDSUB_pd(h1_real, x4, _mm256_shuffle_pd(tmp4, tmp4, 0x5));
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
#else
	x4 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x4), _mm256_shuffle_pd(tmp4, tmp4, 0x5));
#endif

	h1_real = _mm256_broadcast_sd(&hh_dbl[ldh*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[(ldh*2)+1]);
	h2_real = _mm256_broadcast_sd(&hh_dbl[ldh*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[(ldh*2)+1]);

	h1_real = _mm256_xor_pd(h1_real, sign);
	h1_imag = _mm256_xor_pd(h1_imag, sign);
	h2_real = _mm256_xor_pd(h2_real, sign);
	h2_imag = _mm256_xor_pd(h2_imag, sign);

	__m128d tmp_s_128 = _mm_loadu_pd(s_dbl);
	tmp2 = _mm256_broadcast_pd(&tmp_s_128);
	tmp1 = _mm256_mul_pd(h2_imag, tmp2);
353
354
#ifdef __ELPA_USE_FMA__
	tmp2 = _mm256_FMADDSUB_pd(h2_real, tmp2, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
355
356
357
358
359
360
361
362
#else
	tmp2 = _mm256_addsub_pd( _mm256_mul_pd(h2_real, tmp2), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	_mm_storeu_pd(s_dbl, _mm256_castpd256_pd128(tmp2));
	h2_real = _mm256_broadcast_sd(&s_dbl[0]);
	h2_imag = _mm256_broadcast_sd(&s_dbl[1]);

	tmp1 = _mm256_mul_pd(h1_imag, y1);
363
364
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_FMADDSUB_pd(h1_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
365
366
367
368
#else
	y1 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, y2);
369
370
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_FMADDSUB_pd(h1_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5));
371
372
373
374
#else
	y2 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, y3);
375
376
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_FMADDSUB_pd(h1_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5));
377
378
379
380
#else
	y3 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5));
#endif
	tmp4 = _mm256_mul_pd(h1_imag, y4);
381
382
#ifdef __ELPA_USE_FMA__
	y4 = _mm256_FMADDSUB_pd(h1_real, y4, _mm256_shuffle_pd(tmp4, tmp4, 0x5));
383
384
385
386
387
#else
	y4 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y4), _mm256_shuffle_pd(tmp4, tmp4, 0x5));
#endif

	tmp1 = _mm256_mul_pd(h2_imag, x1);
388
389
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_add_pd(y1, _mm256_FMADDSUB_pd(h2_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
390
391
392
393
#else
	y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, x2);
394
395
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_add_pd(y2, _mm256_FMADDSUB_pd(h2_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
396
397
398
399
#else
	y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, x3);
400
401
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_add_pd(y3, _mm256_FMADDSUB_pd(h2_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
402
403
404
405
#else
	y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	tmp4 = _mm256_mul_pd(h2_imag, x4);
406
407
#ifdef __ELPA_USE_FMA__
	y4 = _mm256_add_pd(y4, _mm256_FMADDSUB_pd(h2_real, x4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
#else
	y4 = _mm256_add_pd(y4, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

	q1 = _mm256_load_pd(&q_dbl[0]);
	q2 = _mm256_load_pd(&q_dbl[4]);
	q3 = _mm256_load_pd(&q_dbl[8]);
	q4 = _mm256_load_pd(&q_dbl[12]);

	q1 = _mm256_add_pd(q1, y1);
	q2 = _mm256_add_pd(q2, y2);
	q3 = _mm256_add_pd(q3, y3);
	q4 = _mm256_add_pd(q4, y4);

	_mm256_store_pd(&q_dbl[0], q1);
	_mm256_store_pd(&q_dbl[4], q2);
	_mm256_store_pd(&q_dbl[8], q3);
	_mm256_store_pd(&q_dbl[12], q4);

	h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+1)*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+1)*2)+1]);

	q1 = _mm256_load_pd(&q_dbl[(ldq*2)+0]);
	q2 = _mm256_load_pd(&q_dbl[(ldq*2)+4]);
	q3 = _mm256_load_pd(&q_dbl[(ldq*2)+8]);
	q4 = _mm256_load_pd(&q_dbl[(ldq*2)+12]);

	q1 = _mm256_add_pd(q1, x1);
	q2 = _mm256_add_pd(q2, x2);
	q3 = _mm256_add_pd(q3, x3);
	q4 = _mm256_add_pd(q4, x4);

	tmp1 = _mm256_mul_pd(h2_imag, y1);
441
442
#ifdef __ELPA_USE_FMA__
	q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h2_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
443
444
445
446
447
#else
	q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, y2);
#ifdef __FMA4_
448
	q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h2_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
449
450
451
452
#else
	q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, y3);
453
454
#ifdef __ELPA_USE_FMA__
	q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h2_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
455
456
457
458
#else
	q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	tmp4 = _mm256_mul_pd(h2_imag, y4);
459
460
#ifdef __ELPA_USE_FMA__
	q4 = _mm256_add_pd(q4, _mm256_FMADDSUB_pd(h2_real, y4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
#else
	q4 = _mm256_add_pd(q4, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

	_mm256_store_pd(&q_dbl[(ldq*2)+0], q1);
	_mm256_store_pd(&q_dbl[(ldq*2)+4], q2);
	_mm256_store_pd(&q_dbl[(ldq*2)+8], q3);
	_mm256_store_pd(&q_dbl[(ldq*2)+12], q4);

	for (i = 2; i < nb; i++)
	{
		q1 = _mm256_load_pd(&q_dbl[(2*i*ldq)+0]);
		q2 = _mm256_load_pd(&q_dbl[(2*i*ldq)+4]);
		q3 = _mm256_load_pd(&q_dbl[(2*i*ldq)+8]);
		q4 = _mm256_load_pd(&q_dbl[(2*i*ldq)+12]);

		h1_real = _mm256_broadcast_sd(&hh_dbl[(i-1)*2]);
		h1_imag = _mm256_broadcast_sd(&hh_dbl[((i-1)*2)+1]);

		tmp1 = _mm256_mul_pd(h1_imag, x1);
481
482
#ifdef __ELPA_USE_FMA__
		q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
483
484
485
486
#else
		q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h1_imag, x2);
487
488
#ifdef __ELPA_USE_FMA__
		q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
489
490
491
492
#else
		q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h1_imag, x3);
493
494
#ifdef __ELPA_USE_FMA__
		q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
495
496
497
498
#else
		q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
		tmp4 = _mm256_mul_pd(h1_imag, x4);
499
500
#ifdef __ELPA_USE_FMA__
		q4 = _mm256_add_pd(q4, _mm256_FMADDSUB_pd(h1_real, x4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
501
502
503
504
505
506
507
508
#else
		q4 = _mm256_add_pd(q4, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

		h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+i)*2]);
		h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+i)*2)+1]);

		tmp1 = _mm256_mul_pd(h2_imag, y1);
509
510
#ifdef __ELPA_USE_FMA__
		q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h2_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
511
512
513
514
#else
		q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h2_imag, y2);
515
516
#ifdef __ELPA_USE_FMA__
		q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h2_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
517
518
519
520
#else
		q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h2_imag, y3);
521
522
#ifdef __ELPA_USE_FMA__
		q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h2_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
523
524
525
526
#else
		q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
		tmp4 = _mm256_mul_pd(h2_imag, y4);
527
528
#ifdef __ELPA_USE_FMA__
		q4 = _mm256_add_pd(q4, _mm256_FMADDSUB_pd(h2_real, y4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
#else
		q4 = _mm256_add_pd(q4, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

		_mm256_store_pd(&q_dbl[(2*i*ldq)+0], q1);
		_mm256_store_pd(&q_dbl[(2*i*ldq)+4], q2);
		_mm256_store_pd(&q_dbl[(2*i*ldq)+8], q3);
		_mm256_store_pd(&q_dbl[(2*i*ldq)+12], q4);
	}
	h1_real = _mm256_broadcast_sd(&hh_dbl[(nb-1)*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[((nb-1)*2)+1]);

	q1 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+0]);
	q2 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+4]);
	q3 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+8]);
	q4 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+12]);

	tmp1 = _mm256_mul_pd(h1_imag, x1);
547
548
#ifdef __ELPA_USE_FMA__
	q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
549
550
551
552
#else
	q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, x2);
553
554
#ifdef __ELPA_USE_FMA__
	q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
555
556
557
558
#else
	q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, x3);
559
560
#ifdef __ELPA_USE_FMA__
	q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
561
562
563
564
#else
	q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	tmp4 = _mm256_mul_pd(h1_imag, x4);
565
566
#ifdef __ELPA_USE_FMA__
	q4 = _mm256_add_pd(q4, _mm256_FMADDSUB_pd(h1_real, x4, _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
567
568
569
570
571
572
573
574
575
576
#else
	q4 = _mm256_add_pd(q4, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x4), _mm256_shuffle_pd(tmp4, tmp4, 0x5)));
#endif

	_mm256_store_pd(&q_dbl[(2*nb*ldq)+0], q1);
	_mm256_store_pd(&q_dbl[(2*nb*ldq)+4], q2);
	_mm256_store_pd(&q_dbl[(2*nb*ldq)+8], q3);
	_mm256_store_pd(&q_dbl[(2*nb*ldq)+12], q4);
}

577
static __forceinline void hh_trafo_complex_kernel_6_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s)
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
{
	double* q_dbl = (double*)q;
	double* hh_dbl = (double*)hh;
	double* s_dbl = (double*)(&s);

	__m256d x1, x2, x3;
	__m256d y1, y2, y3;
	__m256d q1, q2, q3;
	__m256d h1_real, h1_imag, h2_real, h2_imag;
	__m256d tmp1, tmp2, tmp3;
	int i=0;

	__m256d sign = (__m256d)_mm256_set_epi64x(0x8000000000000000, 0x8000000000000000, 0x8000000000000000, 0x8000000000000000);

	x1 = _mm256_load_pd(&q_dbl[(2*ldq)+0]);
	x2 = _mm256_load_pd(&q_dbl[(2*ldq)+4]);
	x3 = _mm256_load_pd(&q_dbl[(2*ldq)+8]);

	h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+1)*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+1)*2)+1]);
598
#ifndef __ELPA_USE_FMA__
599
600
601
602
603
604
605
606
607
	// conjugate
	h2_imag = _mm256_xor_pd(h2_imag, sign);
#endif

	y1 = _mm256_load_pd(&q_dbl[0]);
	y2 = _mm256_load_pd(&q_dbl[4]);
	y3 = _mm256_load_pd(&q_dbl[8]);

	tmp1 = _mm256_mul_pd(h2_imag, x1);
608
609
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_add_pd(y1, _mm256_FMSUBADD_pd(h2_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
610
611
612
613
#else
	y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, x2);
614
615
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_add_pd(y2, _mm256_FMSUBADD_pd(h2_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
616
617
618
619
#else
	y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, x3);
620
621
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_add_pd(y3, _mm256_FMSUBADD_pd(h2_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
622
623
624
625
626
627
628
629
630
631
632
633
#else
	y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

	for (i = 2; i < nb; i++)
	{
		q1 = _mm256_load_pd(&q_dbl[(2*i*ldq)+0]);
		q2 = _mm256_load_pd(&q_dbl[(2*i*ldq)+4]);
		q3 = _mm256_load_pd(&q_dbl[(2*i*ldq)+8]);

		h1_real = _mm256_broadcast_sd(&hh_dbl[(i-1)*2]);
		h1_imag = _mm256_broadcast_sd(&hh_dbl[((i-1)*2)+1]);
634
#ifndef __ELPA_USE_FMA__
635
636
637
638
639
		// conjugate
		h1_imag = _mm256_xor_pd(h1_imag, sign);
#endif

		tmp1 = _mm256_mul_pd(h1_imag, q1);
640
641
#ifdef __ELPA_USE_FMA__
		x1 = _mm256_add_pd(x1, _mm256_FMSUBADD_pd(h1_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
642
643
644
645
#else
		x1 = _mm256_add_pd(x1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h1_imag, q2);
646
647
#ifdef __ELPA_USE_FMA__
		x2 = _mm256_add_pd(x2, _mm256_FMSUBADD_pd(h1_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
648
649
650
651
#else
		x2 = _mm256_add_pd(x2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h1_imag, q3);
652
653
#ifdef __ELPA_USE_FMA__
		x3 = _mm256_add_pd(x3, _mm256_FMSUBADD_pd(h1_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
654
655
656
657
658
659
#else
		x3 = _mm256_add_pd(x3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

		h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+i)*2]);
		h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+i)*2)+1]);
660
#ifndef __ELPA_USE_FMA__
661
662
663
664
665
		// conjugate
		h2_imag = _mm256_xor_pd(h2_imag, sign);
#endif

		tmp1 = _mm256_mul_pd(h2_imag, q1);
666
667
#ifdef __ELPA_USE_FMA__
		y1 = _mm256_add_pd(y1, _mm256_FMSUBADD_pd(h2_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
668
669
670
671
#else
		y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h2_imag, q2);
672
673
#ifdef __ELPA_USE_FMA__
		y2 = _mm256_add_pd(y2, _mm256_FMSUBADD_pd(h2_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
674
675
676
677
#else
		y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h2_imag, q3);
678
679
#ifdef __ELPA_USE_FMA__
		y3 = _mm256_add_pd(y3, _mm256_FMSUBADD_pd(h2_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
680
681
682
683
684
685
686
#else
		y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif
	}

	h1_real = _mm256_broadcast_sd(&hh_dbl[(nb-1)*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[((nb-1)*2)+1]);
687
#ifndef __ELPA_USE_FMA__
688
689
690
691
692
693
694
695
696
	// conjugate
	h1_imag = _mm256_xor_pd(h1_imag, sign);
#endif

	q1 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+0]);
	q2 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+4]);
	q3 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+8]);

	tmp1 = _mm256_mul_pd(h1_imag, q1);
697
698
#ifdef __ELPA_USE_FMA__
	x1 = _mm256_add_pd(x1, _mm256_FMSUBADD_pd(h1_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
699
700
701
702
#else
	x1 = _mm256_add_pd(x1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, q2);
703
704
#ifdef __ELPA_USE_FMA__
	x2 = _mm256_add_pd(x2, _mm256_FMSUBADD_pd(h1_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
705
706
707
708
#else
	x2 = _mm256_add_pd(x2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, q3);
709
710
#ifdef __ELPA_USE_FMA__
	x3 = _mm256_add_pd(x3, _mm256_FMSUBADD_pd(h1_real, q3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
711
712
713
714
715
716
717
718
719
720
#else
	x3 = _mm256_add_pd(x3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

	h1_real = _mm256_broadcast_sd(&hh_dbl[0]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[1]);
	h1_real = _mm256_xor_pd(h1_real, sign);
	h1_imag = _mm256_xor_pd(h1_imag, sign);

	tmp1 = _mm256_mul_pd(h1_imag, x1);
721
722
#ifdef __ELPA_USE_FMA__
	x1 = _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
723
724
725
726
#else
	x1 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, x2);
727
728
#ifdef __ELPA_USE_FMA__
	x2 = _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5));
729
730
731
732
#else
	x2 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, x3);
733
734
#ifdef __ELPA_USE_FMA__
	x3 = _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5));
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
#else
	x3 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5));
#endif

	h1_real = _mm256_broadcast_sd(&hh_dbl[ldh*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[(ldh*2)+1]);
	h2_real = _mm256_broadcast_sd(&hh_dbl[ldh*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[(ldh*2)+1]);

	h1_real = _mm256_xor_pd(h1_real, sign);
	h1_imag = _mm256_xor_pd(h1_imag, sign);
	h2_real = _mm256_xor_pd(h2_real, sign);
	h2_imag = _mm256_xor_pd(h2_imag, sign);

	__m128d tmp_s_128 = _mm_loadu_pd(s_dbl);
	tmp2 = _mm256_broadcast_pd(&tmp_s_128);
	tmp1 = _mm256_mul_pd(h2_imag, tmp2);
752
753
#ifdef __ELPA_USE_FMA__
	tmp2 = _mm256_FMADDSUB_pd(h2_real, tmp2, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
754
755
756
757
758
759
760
761
#else
	tmp2 = _mm256_addsub_pd( _mm256_mul_pd(h2_real, tmp2), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	_mm_storeu_pd(s_dbl, _mm256_castpd256_pd128(tmp2));
	h2_real = _mm256_broadcast_sd(&s_dbl[0]);
	h2_imag = _mm256_broadcast_sd(&s_dbl[1]);

	tmp1 = _mm256_mul_pd(h1_imag, y1);
762
763
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_FMADDSUB_pd(h1_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5));
764
765
766
767
#else
	y1 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, y2);
768
769
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_FMADDSUB_pd(h1_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5));
770
771
772
773
#else
	y2 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, y3);
774
775
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_FMADDSUB_pd(h1_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5));
776
777
778
779
780
#else
	y3 = _mm256_addsub_pd( _mm256_mul_pd(h1_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5));
#endif

	tmp1 = _mm256_mul_pd(h2_imag, x1);
781
782
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_add_pd(y1, _mm256_FMADDSUB_pd(h2_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
783
784
785
786
#else
	y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, x2);
787
788
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_add_pd(y2, _mm256_FMADDSUB_pd(h2_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
789
790
791
792
#else
	y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, x3);
793
794
#ifdef __ELPA_USE_FMA__
	y3 = _mm256_add_pd(y3, _mm256_FMADDSUB_pd(h2_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
#else
	y3 = _mm256_add_pd(y3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

	q1 = _mm256_load_pd(&q_dbl[0]);
	q2 = _mm256_load_pd(&q_dbl[4]);
	q3 = _mm256_load_pd(&q_dbl[8]);

	q1 = _mm256_add_pd(q1, y1);
	q2 = _mm256_add_pd(q2, y2);
	q3 = _mm256_add_pd(q3, y3);

	_mm256_store_pd(&q_dbl[0], q1);
	_mm256_store_pd(&q_dbl[4], q2);
	_mm256_store_pd(&q_dbl[8], q3);

	h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+1)*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+1)*2)+1]);

	q1 = _mm256_load_pd(&q_dbl[(ldq*2)+0]);
	q2 = _mm256_load_pd(&q_dbl[(ldq*2)+4]);
	q3 = _mm256_load_pd(&q_dbl[(ldq*2)+8]);

	q1 = _mm256_add_pd(q1, x1);
	q2 = _mm256_add_pd(q2, x2);
	q3 = _mm256_add_pd(q3, x3);

	tmp1 = _mm256_mul_pd(h2_imag, y1);
823
824
#ifdef __ELPA_USE_FMA__
	q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h2_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
825
826
827
828
829
#else
	q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, y2);
#ifdef __FMA4_
830
	q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h2_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
831
832
833
834
#else
	q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h2_imag, y3);
835
836
#ifdef __ELPA_USE_FMA__
	q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h2_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
#else
	q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

	_mm256_store_pd(&q_dbl[(ldq*2)+0], q1);
	_mm256_store_pd(&q_dbl[(ldq*2)+4], q2);
	_mm256_store_pd(&q_dbl[(ldq*2)+8], q3);

	for (i = 2; i < nb; i++)
	{
		q1 = _mm256_load_pd(&q_dbl[(2*i*ldq)+0]);
		q2 = _mm256_load_pd(&q_dbl[(2*i*ldq)+4]);
		q3 = _mm256_load_pd(&q_dbl[(2*i*ldq)+8]);

		h1_real = _mm256_broadcast_sd(&hh_dbl[(i-1)*2]);
		h1_imag = _mm256_broadcast_sd(&hh_dbl[((i-1)*2)+1]);

		tmp1 = _mm256_mul_pd(h1_imag, x1);
855
856
#ifdef __ELPA_USE_FMA__
		q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
857
858
859
860
#else
		q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h1_imag, x2);
861
862
#ifdef __ELPA_USE_FMA__
		q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
863
864
865
866
#else
		q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h1_imag, x3);
867
868
#ifdef __ELPA_USE_FMA__
		q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
869
870
871
872
873
874
875
876
#else
		q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

		h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+i)*2]);
		h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+i)*2)+1]);

		tmp1 = _mm256_mul_pd(h2_imag, y1);
877
878
#ifdef __ELPA_USE_FMA__
		q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h2_real, y1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
879
880
881
882
#else
		q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h2_imag, y2);
883
884
#ifdef __ELPA_USE_FMA__
		q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h2_real, y2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
885
886
887
888
#else
		q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
		tmp3 = _mm256_mul_pd(h2_imag, y3);
889
890
#ifdef __ELPA_USE_FMA__
		q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h2_real, y3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
#else
		q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h2_real, y3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

		_mm256_store_pd(&q_dbl[(2*i*ldq)+0], q1);
		_mm256_store_pd(&q_dbl[(2*i*ldq)+4], q2);
		_mm256_store_pd(&q_dbl[(2*i*ldq)+8], q3);
	}
	h1_real = _mm256_broadcast_sd(&hh_dbl[(nb-1)*2]);
	h1_imag = _mm256_broadcast_sd(&hh_dbl[((nb-1)*2)+1]);

	q1 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+0]);
	q2 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+4]);
	q3 = _mm256_load_pd(&q_dbl[(2*nb*ldq)+8]);

	tmp1 = _mm256_mul_pd(h1_imag, x1);
907
908
#ifdef __ELPA_USE_FMA__
	q1 = _mm256_add_pd(q1, _mm256_FMADDSUB_pd(h1_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
909
910
911
912
#else
	q1 = _mm256_add_pd(q1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h1_imag, x2);
913
914
#ifdef __ELPA_USE_FMA__
	q2 = _mm256_add_pd(q2, _mm256_FMADDSUB_pd(h1_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
915
916
917
918
#else
	q2 = _mm256_add_pd(q2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif
	tmp3 = _mm256_mul_pd(h1_imag, x3);
919
920
#ifdef __ELPA_USE_FMA__
	q3 = _mm256_add_pd(q3, _mm256_FMADDSUB_pd(h1_real, x3, _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
921
922
923
924
925
926
927
928
929
#else
	q3 = _mm256_add_pd(q3, _mm256_addsub_pd( _mm256_mul_pd(h1_real, x3), _mm256_shuffle_pd(tmp3, tmp3, 0x5)));
#endif

	_mm256_store_pd(&q_dbl[(2*nb*ldq)+0], q1);
	_mm256_store_pd(&q_dbl[(2*nb*ldq)+4], q2);
	_mm256_store_pd(&q_dbl[(2*nb*ldq)+8], q3);
}

930
static __forceinline void hh_trafo_complex_kernel_4_AVX_2hv_double(double complex* q, double complex* hh, int nb, int ldq, int ldh, double complex s)
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
{
	double* q_dbl = (double*)q;
	double* hh_dbl = (double*)hh;
	double* s_dbl = (double*)(&s);

	__m256d x1, x2;
	__m256d y1, y2;
	__m256d q1, q2;
	__m256d h1_real, h1_imag, h2_real, h2_imag;
	__m256d tmp1, tmp2;
	int i=0;

	__m256d sign = (__m256d)_mm256_set_epi64x(0x8000000000000000, 0x8000000000000000, 0x8000000000000000, 0x8000000000000000);

	x1 = _mm256_load_pd(&q_dbl[(2*ldq)+0]);
	x2 = _mm256_load_pd(&q_dbl[(2*ldq)+4]);

	h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+1)*2]);
	h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+1)*2)+1]);
950
#ifndef __ELPA_USE_FMA__
951
952
953
954
955
956
957
958
	// conjugate
	h2_imag = _mm256_xor_pd(h2_imag, sign);
#endif

	y1 = _mm256_load_pd(&q_dbl[0]);
	y2 = _mm256_load_pd(&q_dbl[4]);

	tmp1 = _mm256_mul_pd(h2_imag, x1);
959
960
#ifdef __ELPA_USE_FMA__
	y1 = _mm256_add_pd(y1, _mm256_FMSUBADD_pd(h2_real, x1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
961
962
963
964
#else
	y1 = _mm256_add_pd(y1, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
	tmp2 = _mm256_mul_pd(h2_imag, x2);
965
966
#ifdef __ELPA_USE_FMA__
	y2 = _mm256_add_pd(y2, _mm256_FMSUBADD_pd(h2_real, x2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
967
968
969
970
971
972
973
974
975
976
977
#else
	y2 = _mm256_add_pd(y2, _mm256_addsub_pd( _mm256_mul_pd(h2_real, x2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif

	for (i = 2; i < nb; i++)
	{
		q1 = _mm256_load_pd(&q_dbl[(2*i*ldq)+0]);
		q2 = _mm256_load_pd(&q_dbl[(2*i*ldq)+4]);

		h1_real = _mm256_broadcast_sd(&hh_dbl[(i-1)*2]);
		h1_imag = _mm256_broadcast_sd(&hh_dbl[((i-1)*2)+1]);
978
#ifndef __ELPA_USE_FMA__
979
980
981
982
983
		// conjugate
		h1_imag = _mm256_xor_pd(h1_imag, sign);
#endif

		tmp1 = _mm256_mul_pd(h1_imag, q1);
984
985
#ifdef __ELPA_USE_FMA__
		x1 = _mm256_add_pd(x1, _mm256_FMSUBADD_pd(h1_real, q1, _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
986
987
988
989
#else
		x1 = _mm256_add_pd(x1, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q1), _mm256_shuffle_pd(tmp1, tmp1, 0x5)));
#endif
		tmp2 = _mm256_mul_pd(h1_imag, q2);
990
991
#ifdef __ELPA_USE_FMA__
		x2 = _mm256_add_pd(x2, _mm256_FMSUBADD_pd(h1_real, q2, _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
992
993
994
995
996
997
#else
		x2 = _mm256_add_pd(x2, _mm256_addsub_pd( _mm256_mul_pd(h1_real, q2), _mm256_shuffle_pd(tmp2, tmp2, 0x5)));
#endif

		h2_real = _mm256_broadcast_sd(&hh_dbl[(ldh+i)*2]);
		h2_imag = _mm256_broadcast_sd(&hh_dbl[((ldh+i)*2)+1]);
998
#ifndef __ELPA_USE_FMA__
999
1000
		// conjugate
		h2_imag = _mm256_xor_pd(h2_imag, sign);