real_128bit_BLOCK_template.c

//    This file is part of ELPA.
//
//    The ELPA library was originally created by the ELPA consortium,
//    consisting of the following organizations:
//
//    - Max Planck Computing and Data Facility (MPCDF), formerly known as
//      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
//    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
//      Informatik,
//    - Technische Universität München, Lehrstuhl für Informatik mit
//      Schwerpunkt Wissenschaftliches Rechnen ,
//    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
//    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
//      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
//      and
//    - IBM Deutschland GmbH
//
//
//    This particular source code file contains additions, changes and
//    enhancements authored by Intel Corporation which is not part of
//    the ELPA consortium.
//
//    More information can be found here:
//    http://elpa.mpcdf.mpg.de/
//
//    ELPA is free software: you can redistribute it and/or modify
//    it under the terms of the version 3 of the license of the
//    GNU Lesser General Public License as published by the Free
//    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.
//
// Author: Andreas Marek, MPCDF, based on the double precision case of A. Heinecke
//
#include "config-f90.h"

#define CONCAT_8ARGS(a, b, c, d, e, f, g, h) CONCAT2_8ARGS(a, b, c, d, e, f, g, h)
#define CONCAT2_8ARGS(a, b, c, d, e, f, g, h) a ## b ## c ## d ## e ## f ## g ## h

#define CONCAT_7ARGS(a, b, c, d, e, f, g) CONCAT2_7ARGS(a, b, c, d, e, f, g)
#define CONCAT2_7ARGS(a, b, c, d, e, f, g) a ## b ## c ## d ## e ## f ## g

#define CONCAT_6ARGS(a, b, c, d, e, f) CONCAT2_6ARGS(a, b, c, d, e, f)
#define CONCAT2_6ARGS(a, b, c, d, e, f) a ## b ## c ## d ## e ## f

#define CONCAT_5ARGS(a, b, c, d, e) CONCAT2_5ARGS(a, b, c, d, e)
#define CONCAT2_5ARGS(a, b, c, d, e) a ## b ## c ## d ## e

#define CONCAT_4ARGS(a, b, c, d) CONCAT2_4ARGS(a, b, c, d)
#define CONCAT2_4ARGS(a, b, c, d) a ## b ## c ## d

#define CONCAT_3ARGS(a, b, c) CONCAT2_3ARGS(a, b, c)
#define CONCAT2_3ARGS(a, b, c) a ## b ## c

#if VEC_SET == 128 || VEC_SET == 256
#include <x86intrin.h>
#endif
#if VEC_SET == 1281
#include <fjmfunc.h>
#include <emmintrin.h>
#endif
#include <stdio.h>
#include <stdlib.h>

#ifdef BLOCK6
#define PREFIX hexa
#define BLOCK 6
#endif

#ifdef BLOCK4
#define PREFIX quad
#define BLOCK 4
#endif

#ifdef BLOCK2
#define PREFIX double
#define BLOCK 2
#endif

#if VEC_SET == 128
#define SIMD_SET SSE
#endif

#if VEC_SET == 1281
#define SIMD_SET SPARC64
#endif

#if VEC_SET == 256
#define SIMD_SET AVX_AVX2
#endif

#define __forceinline __attribute__((always_inline)) static

#if VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define offset 2
#define __SIMD_DATATYPE __m128d
#define _SIMD_LOAD _mm_load_pd
#define _SIMD_STORE _mm_store_pd
#define _SIMD_ADD _mm_add_pd
#define _SIMD_MUL _mm_mul_pd
#define _SIMD_SUB _mm_sub_pd
#define _SIMD_XOR _mm_xor_pd
#endif
#ifdef SINGLE_PRECISION_REAL
#define offset 4
#define __SIMD_DATATYPE __m128
#define _SIMD_LOAD _mm_load_ps
#define _SIMD_STORE _mm_store_ps
#define _SIMD_ADD _mm_add_ps
#define _SIMD_MUL _mm_mul_ps
#define _SIMD_SUB _mm_sub_ps
#define _SIMD_XOR _mm_xor_ps
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281 */

#if VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define offset 4
#define __SIMD_DATATYPE __m256d
#define _SIMD_LOAD _mm256_load_pd
#define _SIMD_STORE _mm256_store_pd
#define _SIMD_ADD _mm256_add_pd
#define _SIMD_MUL _mm256_mul_pd
#define _SIMD_SUB _mm256_sub_pd
#define _SIMD_XOR _mm256_xor_pd
#define _SIMD_BROADCAST _mm256_broadcast_sd
#ifdef HAVE_AVX2
#ifdef __FMA4__
#define __ELPA_USE_FMA__
#define _mm256_FMA_pd(a,b,c) _mm256_macc_pd(a,b,c)
#define _mm256_NFMA_pd(a,b,c) _mm256_nmacc_pd(a,b,c)
#error "This should be prop _mm256_msub_pd instead of _mm256_msub"
#define _mm256_FMSUB_pd(a,b,c) _mm256_msub(a,b,c)
#endif /* __FMA4__ */
#ifdef __AVX2__
#define __ELPA_USE_FMA__
#define _mm256_FMA_pd(a,b,c) _mm256_fmadd_pd(a,b,c)
#define _mm256_NFMA_pd(a,b,c) _mm256_fnmadd_pd(a,b,c)
#define _mm256_FMSUB_pd(a,b,c) _mm256_fmsub_pd(a,b,c)
#endif /* __AVX2__ */
#ifdef __ELPA_USE_FMA__
#define _SIMD_FMA _mm256_FMA_pd
#define _SIMD_NFMA _mm256_NFMA_pd
#define _SIMD_FMSUB _mm256_FMSUB_pd
#endif
#endif /* HAVE_AVX2 */
#endif /* DOUBLE_PRECISION_REAL */

#ifdef SINGLE_PRECISION_REAL
#define offset 8
#define __SIMD_DATATYPE __m256
#define _SIMD_LOAD _mm256_load_ps
#define _SIMD_STORE _mm256_store_ps
#define _SIMD_ADD _mm256_add_ps
#define _SIMD_MUL _mm256_mul_ps
#define _SIMD_SUB _mm256_sub_ps
#define _SIMD_XOR _mm256_xor_ps
#define _SIMD_BROADCAST _mm256_broadcast_ss
#ifdef HAVE_AVX2
#ifdef __FMA4__
#define __ELPA_USE_FMA__
#define _mm256_FMA_ps(a,b,c) _mm256_macc_ps(a,b,c)
#define _mm256_NFMA_ps(a,b,c) _mm256_nmacc_ps(a,b,c)
#error "This should be prop _mm256_msub_ps instead of _mm256_msub"
#define _mm256_FMSUB_ps(a,b,c) _mm256_msub(a,b,c)
#endif /* __FMA4__ */
#ifdef __AVX2__
#define __ELPA_USE_FMA__
#define _mm256_FMA_ps(a,b,c) _mm256_fmadd_ps(a,b,c)
#define _mm256_NFMA_ps(a,b,c) _mm256_fnmadd_ps(a,b,c)
#define _mm256_FMSUB_ps(a,b,c) _mm256_fmsub_ps(a,b,c)
#endif /* __AVX2__ */
#ifdef __ELPA_USE_FMA__
#define _SIMD_FMA _mm256_FMA_ps
#define _SIMD_NFMA _mm256_NFMA_ps
#define _SIMD_FMSUB _mm256_FMSUB_ps
#endif
#endif /* HAVE_AVX2 */
#endif /* SINGLE_PRECISION_REAL */
#endif /* VEC_SET == 256 */

#ifdef DOUBLE_PRECISION_REAL
#define WORD_LENGTH double
#define DATA_TYPE double
#define DATA_TYPE_PTR double*
#endif
#ifdef SINGLE_PRECISION_REAL
#define WORD_LENGTH single
#define DATA_TYPE float
#define DATA_TYPE_PTR float*
#endif


#ifdef DOUBLE_PRECISION_REAL


#define _SSE_SET _mm_set_pd
#define _SSE_SET1 _mm_set1_pd
#define _SSE_SET _mm_set_pd
#endif

#ifdef SINGLE_PRECISION_REAL

#define _SSE_SET _mm_set_ps
#define _SSE_SET1 _mm_set1_ps
#define _SSE_SET _mm_set_ps

#endif

#if VEC_SET == 128
#undef __AVX__
#endif


#if VEC_SET == 128 || VEC_SET == 1281
//Forward declaration
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 2
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 4
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281 */

#if VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 8
#endif
#endif /* VEC_SET == 256 */
__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh, 
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

#if VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 8
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281 */

#if VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 8
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 16
#endif
#endif /* VEC_SET == 256 */
__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh, 
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

#if VEC_SET == 128 || VEC_SET == 1281 
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 6
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 12
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281  */

#if VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 12
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 24
#endif
#endif /* VEC_SET == 256 */

__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh,
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

#if VEC_SET == 128 || VEC_SET == 1281 
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 8
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 16
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281  */

#if VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 16
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 32
#endif
#endif /* VEC_SET == 256 */

__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh, 
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 10
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 20
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 20
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 40
#endif
#endif /*  VEC_SET == 256 */


__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh, 
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 12
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 24
#endif
#endif /* VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 24
#endif
#ifdef SINGLE_PRECISION_REAL
#undef ROW_LENGTH
#define ROW_LENGTH 48
#endif
#endif /*  VEC_SET == 256 */

__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh,
#ifdef BLOCK2
	DATA_TYPE s);
#endif
#ifdef BLOCK4
	DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4);
#endif
#ifdef BLOCK6
	DATA_TYPE_PTR scalarprods);
#endif

void CONCAT_7ARGS(PREFIX,_hh_trafo_real_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int* pnb, int* pnq, int* pldq, int* pldh);

/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine double_hh_trafo_real_SSE_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="double_hh_trafo_real_SSE_2hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine double_hh_trafo_real_SSE_2hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="double_hh_trafo_real_SSE_2hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_float)  :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine double_hh_trafo_real_SPARC64_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="double_hh_trafo_real_SPARC64_2hv_double")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_double) :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine double_hh_trafo_real_SPARC64_2hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="double_hh_trafo_real_SPARC64_2hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_float)  :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine double_hh_trafo_real_AVX_AVX2_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="double_hh_trafo_real_AVX_AVX2_2hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine double_hh_trafo_real_AVX_AVX2_2hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="double_hh_trafo_real_AVX_AVX2_2hv_single")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)       :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_float)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine quad_hh_trafo_real_SSE_4hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="quad_hh_trafo_real_SSE_4hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine quad_hh_trafo_real_SSE_4hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="quad_hh_trafo_real_SSE_4hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_float)  :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine quad_hh_trafo_real_SPARC64_4hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="quad_hh_trafo_real_SPARC64_4hv_double")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_double) :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine quad_hh_trafo_real_SPARC64_4hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="quad_hh_trafo_real_SPARC64_4hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_float)  :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine quad_hh_trafo_real_AVX_AVX2_4hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="quad_hh_trafo_real_AVX_AVX2_4hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine quad_hh_trafo_real_AVX_AVX2_4hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>              bind(C, name="quad_hh_trafo_real_AVX_AVX2_4hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int) :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value  :: q
!f>     real(kind=c_float)  :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine hexa_hh_trafo_real_sse_6hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="hexa_hh_trafo_real_SSE_6hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine hexa_hh_trafo_real_sparc64_6hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="hexa_hh_trafo_real_SPARC64_6hv_double")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_double)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#ifdef HAVE_SSE_INTRINSICS
!f> interface
!f>   subroutine hexa_hh_trafo_real_sse_6hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="hexa_hh_trafo_real_SSE_6hv_single")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_float)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#ifdef HAVE_SPARC64_SSE
!f> interface
!f>   subroutine hexa_hh_trafo_real_sparc64_6hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>                                bind(C, name="hexa_hh_trafo_real_SPARC64_6hv_single")
!f>        use, intrinsic :: iso_c_binding
!f>        integer(kind=c_int)        :: pnb, pnq, pldq, pldh
!f>        type(c_ptr), value        :: q
!f>        real(kind=c_float)        :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/

/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine hexa_hh_trafo_real_AVX_AVX2_6hv_double(q, hh, pnb, pnq, pldq, pldh) &
!f>                             bind(C, name="hexa_hh_trafo_real_AVX_AVX2_6hv_double")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int)     :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value      :: q
!f>     real(kind=c_double)     :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/
/*
!f>#if defined(HAVE_AVX) || defined(HAVE_AVX2)
!f> interface
!f>   subroutine hexa_hh_trafo_real_AVX_AVX2_6hv_single(q, hh, pnb, pnq, pldq, pldh) &
!f>                             bind(C, name="hexa_hh_trafo_real_AVX_AVX2_6hv_single")
!f>     use, intrinsic :: iso_c_binding
!f>     integer(kind=c_int)     :: pnb, pnq, pldq, pldh
!f>     type(c_ptr), value      :: q
!f>     real(kind=c_float)      :: hh(pnb,6)
!f>   end subroutine
!f> end interface
!f>#endif
*/


void CONCAT_7ARGS(PREFIX,_hh_trafo_real_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int* pnb, int* pnq, int* pldq, int* pldh)
{
  int i;
  int nb = *pnb;
  int nq = *pldq;
  int ldq = *pldq;
  int ldh = *pldh;
  int worked_on;

  worked_on = 0;

#ifdef BLOCK2
  // calculating scalar product to compute
  // 2 householder vectors simultaneously
  DATA_TYPE s = hh[(ldh)+1]*1.0;
#endif

#ifdef BLOCK4
  // calculating scalar products to compute
  // 4 householder vectors simultaneously
  DATA_TYPE s_1_2 = hh[(ldh)+1];  
  DATA_TYPE s_1_3 = hh[(ldh*2)+2];
  DATA_TYPE s_2_3 = hh[(ldh*2)+1];
  DATA_TYPE s_1_4 = hh[(ldh*3)+3];
  DATA_TYPE s_2_4 = hh[(ldh*3)+2];
  DATA_TYPE s_3_4 = hh[(ldh*3)+1];

  // calculate scalar product of first and fourth householder Vector
  // loop counter = 2
  s_1_2 += hh[2-1] * hh[(2+ldh)];          
  s_2_3 += hh[(ldh)+2-1] * hh[2+(ldh*2)];  
  s_3_4 += hh[(ldh*2)+2-1] * hh[2+(ldh*3)];

  // loop counter = 3
  s_1_2 += hh[3-1] * hh[(3+ldh)];          
  s_2_3 += hh[(ldh)+3-1] * hh[3+(ldh*2)];  
  s_3_4 += hh[(ldh*2)+3-1] * hh[3+(ldh*3)];

  s_1_3 += hh[3-2] * hh[3+(ldh*2)];        
  s_2_4 += hh[(ldh*1)+3-2] * hh[3+(ldh*3)];
#endif /* BLOCK4 */

#ifdef BLOCK6
  // calculating scalar products to compute
  // 6 householder vectors simultaneously
  DATA_TYPE scalarprods[15];

  scalarprods[0] = hh[(ldh+1)];  
  scalarprods[1] = hh[(ldh*2)+2];
  scalarprods[2] = hh[(ldh*2)+1];
  scalarprods[3] = hh[(ldh*3)+3];
  scalarprods[4] = hh[(ldh*3)+2];
  scalarprods[5] = hh[(ldh*3)+1];
  scalarprods[6] = hh[(ldh*4)+4];
  scalarprods[7] = hh[(ldh*4)+3];
  scalarprods[8] = hh[(ldh*4)+2];
  scalarprods[9] = hh[(ldh*4)+1];
  scalarprods[10] = hh[(ldh*5)+5];
  scalarprods[11] = hh[(ldh*5)+4];
  scalarprods[12] = hh[(ldh*5)+3];
  scalarprods[13] = hh[(ldh*5)+2];
  scalarprods[14] = hh[(ldh*5)+1];

  // calculate scalar product of first and fourth householder Vector
  // loop counter = 2
  scalarprods[0] += hh[1] * hh[(2+ldh)];           
  scalarprods[2] += hh[(ldh)+1] * hh[2+(ldh*2)];   
  scalarprods[5] += hh[(ldh*2)+1] * hh[2+(ldh*3)]; 
  scalarprods[9] += hh[(ldh*3)+1] * hh[2+(ldh*4)]; 
  scalarprods[14] += hh[(ldh*4)+1] * hh[2+(ldh*5)];

  // loop counter = 3
  scalarprods[0] += hh[2] * hh[(3+ldh)];          
  scalarprods[2] += hh[(ldh)+2] * hh[3+(ldh*2)];  
  scalarprods[5] += hh[(ldh*2)+2] * hh[3+(ldh*3)];
  scalarprods[9] += hh[(ldh*3)+2] * hh[3+(ldh*4)];
  scalarprods[14] += hh[(ldh*4)+2] * hh[3+(ldh*5)];

  scalarprods[1] += hh[1] * hh[3+(ldh*2)];         
  scalarprods[4] += hh[(ldh*1)+1] * hh[3+(ldh*3)]; 
  scalarprods[8] += hh[(ldh*2)+1] * hh[3+(ldh*4)]; 
  scalarprods[13] += hh[(ldh*3)+1] * hh[3+(ldh*5)];

  // loop counter = 4
  scalarprods[0] += hh[3] * hh[(4+ldh)];           
  scalarprods[2] += hh[(ldh)+3] * hh[4+(ldh*2)];   
  scalarprods[5] += hh[(ldh*2)+3] * hh[4+(ldh*3)]; 
  scalarprods[9] += hh[(ldh*3)+3] * hh[4+(ldh*4)]; 
  scalarprods[14] += hh[(ldh*4)+3] * hh[4+(ldh*5)];

  scalarprods[1] += hh[2] * hh[4+(ldh*2)];         
  scalarprods[4] += hh[(ldh*1)+2] * hh[4+(ldh*3)]; 
  scalarprods[8] += hh[(ldh*2)+2] * hh[4+(ldh*4)]; 
  scalarprods[13] += hh[(ldh*3)+2] * hh[4+(ldh*5)];

  scalarprods[3] += hh[1] * hh[4+(ldh*3)];         
  scalarprods[7] += hh[(ldh)+1] * hh[4+(ldh*4)];   
  scalarprods[12] += hh[(ldh*2)+1] * hh[4+(ldh*5)];

  // loop counter = 5
  scalarprods[0] += hh[4] * hh[(5+ldh)];           
  scalarprods[2] += hh[(ldh)+4] * hh[5+(ldh*2)];   
  scalarprods[5] += hh[(ldh*2)+4] * hh[5+(ldh*3)]; 
  scalarprods[9] += hh[(ldh*3)+4] * hh[5+(ldh*4)]; 
  scalarprods[14] += hh[(ldh*4)+4] * hh[5+(ldh*5)];

  scalarprods[1] += hh[3] * hh[5+(ldh*2)];         
  scalarprods[4] += hh[(ldh*1)+3] * hh[5+(ldh*3)]; 
  scalarprods[8] += hh[(ldh*2)+3] * hh[5+(ldh*4)]; 
  scalarprods[13] += hh[(ldh*3)+3] * hh[5+(ldh*5)];

  scalarprods[3] += hh[2] * hh[5+(ldh*3)];         
  scalarprods[7] += hh[(ldh)+2] * hh[5+(ldh*4)];   
  scalarprods[12] += hh[(ldh*2)+2] * hh[5+(ldh*5)];

  scalarprods[6] += hh[1] * hh[5+(ldh*4)];         
  scalarprods[11] += hh[(ldh)+1] * hh[5+(ldh*5)];  


#endif /* BLOCK6 */

#if VEC_SET == 128 || VEC_SET == 256
  #pragma ivdep
#endif
  for (i = BLOCK; i < nb; i++)
    {
#ifdef BLOCK2
      s += hh[i-1] * hh[(i+ldh)];
#endif
#ifdef BLOCK4
      s_1_2 += hh[i-1] * hh[(i+ldh)];           
      s_2_3 += hh[(ldh)+i-1] * hh[i+(ldh*2)];   
      s_3_4 += hh[(ldh*2)+i-1] * hh[i+(ldh*3)]; 

      s_1_3 += hh[i-2] * hh[i+(ldh*2)];         
      s_2_4 += hh[(ldh*1)+i-2] * hh[i+(ldh*3)]; 

      s_1_4 += hh[i-3] * hh[i+(ldh*3)];         
#endif /* BLOCK4 */
#ifdef BLOCK6
      scalarprods[0] += hh[i-1] * hh[(i+ldh)];           
      scalarprods[2] += hh[(ldh)+i-1] * hh[i+(ldh*2)];   
      scalarprods[5] += hh[(ldh*2)+i-1] * hh[i+(ldh*3)]; 
      scalarprods[9] += hh[(ldh*3)+i-1] * hh[i+(ldh*4)]; 
      scalarprods[14] += hh[(ldh*4)+i-1] * hh[i+(ldh*5)];

      scalarprods[1] += hh[i-2] * hh[i+(ldh*2)];         
      scalarprods[4] += hh[(ldh*1)+i-2] * hh[i+(ldh*3)]; 
      scalarprods[8] += hh[(ldh*2)+i-2] * hh[i+(ldh*4)]; 
      scalarprods[13] += hh[(ldh*3)+i-2] * hh[i+(ldh*5)];

      scalarprods[3] += hh[i-3] * hh[i+(ldh*3)];         
      scalarprods[7] += hh[(ldh)+i-3] * hh[i+(ldh*4)];   
      scalarprods[12] += hh[(ldh*2)+i-3] * hh[i+(ldh*5)];

      scalarprods[6] += hh[i-4] * hh[i+(ldh*4)];         
      scalarprods[11] += hh[(ldh)+i-4] * hh[i+(ldh*5)];  

      scalarprods[10] += hh[i-5] * hh[i+(ldh*5)];        
#endif /* BLOCK6 */

    }

  // Production level kernel calls with padding
#ifdef BLOCK2

#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define STEP_SIZE 12
#define ROW_LENGTH 12
#define UPPER_BOUND 10
#endif
#ifdef SINGLE_PRECISION_REAL
#define STEP_SIZE 24
#define ROW_LENGTH 24
#define UPPER_BOUND 20
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define STEP_SIZE 24
#define ROW_LENGTH 24
#define UPPER_BOUND 20
#endif
#ifdef SINGLE_PRECISION_REAL
#define STEP_SIZE 48
#define ROW_LENGTH 48
#define UPPER_BOUND 40
#endif
#endif /*  VEC_SET == 256 */

  for (i = 0; i < nq - UPPER_BOUND; i+= STEP_SIZE )
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

  if (nq == i)
    {
      return;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 10
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 20
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 20
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 40
#endif
#endif /* VEC_SET == 256 */

  if (nq-i == ROW_LENGTH)
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 16
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 16
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 32
#endif
#endif /* VEC_SET == 256 */

  if (nq-i == ROW_LENGTH)
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 6
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 12
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 12
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 24
#endif
#endif /* VEC_SET == 256 */


  if (nq-i == ROW_LENGTH)
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 16
#endif
#endif /* VEC_SET == 256 */


  if (nq-i == ROW_LENGTH)
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 2
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#endif /* VEC_SET == 256 */

  if (nq-i == ROW_LENGTH)
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_2hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s);
      worked_on += ROW_LENGTH;
    }

#endif /* BLOCK2 */

#ifdef BLOCK4


#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 6
#define STEP_SIZE 6
#define UPPER_BOUND 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 12
#define STEP_SIZE 12
#define UPPER_BOUND 8
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 12
#define STEP_SIZE 12
#define UPPER_BOUND 8
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 24
#define STEP_SIZE 24
#define UPPER_BOUND 16
#endif
#endif /* VEC_SET == 256 */

  for (i = 0; i < nq - UPPER_BOUND; i+= STEP_SIZE )
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_4hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
      worked_on += ROW_LENGTH;
    }

  if (nq == i)
    {
      return;
    }


#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 16
#endif
#endif /* VEC_SET == 256 */

  if (nq-i == ROW_LENGTH )
    {
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_4hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
      worked_on += ROW_LENGTH;
    }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 2
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#endif /* VEC_SET == 256 */

   if (nq-i == ROW_LENGTH )
     {
       CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_4hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
       worked_on += ROW_LENGTH;
     }

#endif /* BLOCK4 */

#ifdef BLOCK6

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#define STEP_SIZE 4
#define UPPER_BOUND 2
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#define STEP_SIZE 8
#define UPPER_BOUND 4
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 8
#define STEP_SIZE 8
#define UPPER_BOUND 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 16
#define STEP_SIZE 16
#define UPPER_BOUND 8
#endif
#endif /* VEC_SET == 256 */

  for (i = 0; i < nq - UPPER_BOUND; i+= STEP_SIZE)
    { 
      CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_6hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, scalarprods);
      worked_on += ROW_LENGTH;
    }
    if (nq == i)
      {
        return;
      }

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 2
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 4
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 8
#endif
#endif /* VEC_SET == 256 */

    if (nq -i == ROW_LENGTH )
      {
        CONCAT_6ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_6hv_,WORD_LENGTH) (&q[i], hh, nb, ldq, ldh, scalarprods);
        worked_on += ROW_LENGTH;
      }
  
#endif /* BLOCK6 */

#ifdef WITH_DEBUG
  if (worked_on != nq)
    {
      printf("Error in real SIMD_SET BLOCK BLOCK kernel %d %d\n", worked_on, nq);
      abort();
    }
#endif
}

#undef ROW_LENGTH
#if  VEC_SET == 128 || VEC_SET == 1281
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 12
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 24
#endif
#endif /*  VEC_SET == 128 || VEC_SET == 1281 */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
#define ROW_LENGTH 24
#endif
#ifdef SINGLE_PRECISION_REAL
#define ROW_LENGTH 48
#endif
#endif /* VEC_SET == 256 */
/*
 * Unrolled kernel that computes
 * ROW_LENGTH rows of Q simultaneously, a
 * matrix Vector product with two householder
 */
#ifdef BLOCK2
/*
 * vectors + a rank 2 update is performed
 */
#endif
#ifdef BLOCK4
/*
 * vectors + a rank 1 update is performed
 */
#endif
__forceinline void CONCAT_8ARGS(hh_trafo_kernel_,ROW_LENGTH,_,SIMD_SET,_,BLOCK,hv_,WORD_LENGTH) (DATA_TYPE_PTR q, DATA_TYPE_PTR hh, int nb, int ldq, int ldh,
#ifdef BLOCK2
               DATA_TYPE s)
#endif
#ifdef BLOCK4
               DATA_TYPE s_1_2, DATA_TYPE s_1_3, DATA_TYPE s_2_3, DATA_TYPE s_1_4, DATA_TYPE s_2_4, DATA_TYPE s_3_4)
#endif 
#ifdef BLOCK6
               DATA_TYPE_PTR scalarprods)
#endif
  {
#ifdef BLOCK2
    /////////////////////////////////////////////////////
    // Matrix Vector Multiplication, Q [10 x nb+1] * hh
    // hh contains two householder vectors, with offset 1
    /////////////////////////////////////////////////////
#endif
#if defined(BLOCK4) || defined(BLOCK6)
    /////////////////////////////////////////////////////
    // Matrix Vector Multiplication, Q [10 x nb+3] * hh
    // hh contains four householder vectors
    /////////////////////////////////////////////////////
#endif

    int i;

#ifdef BLOCK2
#if VEC_SET == 128
    // Needed bit mask for floating point sign flip
#ifdef DOUBLE_PRECISION_REAL
    __SIMD_DATATYPE sign = (__SIMD_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
#endif
#ifdef SINGLE_PRECISION_REAL
    __SIMD_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
#endif
#endif /* HAVE_SSE_INTRINSICS */

#if  VEC_SET == 256
#ifdef DOUBLE_PRECISION_REAL
        __SIMD_DATATYPE sign = (__SIMD_DATATYPE)_mm256_set1_epi64x(0x8000000000000000);
#endif
#ifdef SINGLE_PRECISION_REAL
        __SIMD_DATATYPE sign = (__SIMD_DATATYPE)_mm256_set1_epi32(0x80000000);
#endif
#endif /* VEC_SET == 256 */

    __SIMD_DATATYPE x1 = _SIMD_LOAD(&q[ldq]);
    __SIMD_DATATYPE x2 = _SIMD_LOAD(&q[ldq+offset]);
    __SIMD_DATATYPE x3 = _SIMD_LOAD(&q[ldq+2*offset]);
    __SIMD_DATATYPE x4 = _SIMD_LOAD(&q[ldq+3*offset]);
    __SIMD_DATATYPE x5 = _SIMD_LOAD(&q[ldq+4*offset]);
    __SIMD_DATATYPE x6 = _SIMD_LOAD(&q[ldq+5*offset]);

#if VEC_SET == 128
    __SIMD_DATATYPE h1 = _SSE_SET1(hh[ldh+1]);
#endif
#if VEC_SET == 1281
    __SIMD_DATATYPE h1 = _SSE_SET(hh[ldh+1], hh[ldh+1]);
#endif
#if VEC_SET == 256
    __SIMD_DATATYPE h1 = _SIMD_BROADCAST(&hh[ldh+1]);
#endif
 
    __SIMD_DATATYPE h2;
#ifdef __ELPA_USE_FMA__
    __SIMD_DATATYPE q1 = _SIMD_LOAD(q);
    __SIMD_DATATYPE y1 = _SIMD_FMA(x1, h1, q1);
    __SIMD_DATATYPE q2 = _SIMD_LOAD(&q[offset]);
    __SIMD_DATATYPE y2 = _SIMD_FMA(x2, h1, q2);
    __SIMD_DATATYPE q3 = _SIMD_LOAD(&q[2*offset]);
    __SIMD_DATATYPE y3 = _SIMD_FMA(x3, h1, q3);
    __SIMD_DATATYPE q4 = _SIMD_LOAD(&q[3*offset]);
    __SIMD_DATATYPE y4 = _SIMD_FMA(x4, h1, q4);
    __SIMD_DATATYPE q5 = _SIMD_LOAD(&q[4*offset]);
    __SIMD_DATATYPE y5 = _SIMD_FMA(x5, h1, q5);
    __SIMD_DATATYPE q6 = _SIMD_LOAD(&q[5*offset]);
    __SIMD_DATATYPE y6 = _SIMD_FMA(x6, h1, q6);
#else
    __SIMD_DATATYPE q1 = _SIMD_LOAD(q);
    __SIMD_DATATYPE y1 = _SIMD_ADD(q1, _SIMD_MUL(x1, h1));
    __SIMD_DATATYPE q2 = _SIMD_LOAD(&q[offset]);
    __SIMD_DATATYPE y2 = _SIMD_ADD(q2, _SIMD_MUL(x2, h1));
    __SIMD_DATATYPE q3 = _SIMD_LOAD(&q[2*offset]);
    __SIMD_DATATYPE y3 = _SIMD_ADD(q3, _SIMD_MUL(x3, h1));
    __SIMD_DATATYPE q4 = _SIMD_LOAD(&q[3*offset]);
    __SIMD_DATATYPE y4 = _SIMD_ADD(q4, _SIMD_MUL(x4, h1));
    __SIMD_DATATYPE q5 = _SIMD_LOAD(&q[4*offset]);
    __SIMD_DATATYPE y5 = _SIMD_ADD(q5, _SIMD_MUL(x5, h1));
    __SIMD_DATATYPE q6 = _SIMD_LOAD(&q[5*offset]);
    __SIMD_DATATYPE y6 = _SIMD_ADD(q6, _SIMD_MUL(x6, h1));
#endif
#endif /* BLOCK2 */

#ifdef BLOCK4
    __SIMD_DATATYPE a1_1 = _SIMD_LOAD(&q[ldq*3]);
    __SIMD_DATATYPE a2_1 = _SIMD_LOAD(&q[ldq*2]);
    __SIMD_DATATYPE a3_1 = _SIMD_LOAD(&q[ldq]);  
    __SIMD_DATATYPE a4_1 = _SIMD_LOAD(&q[0]);    

#if VEC_SET == 128
    __SIMD_DATATYPE h_2_1 = _SSE_SET1(hh[ldh+1]);    
    __SIMD_DATATYPE h_3_2 = _SSE_SET1(hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SSE_SET1(hh[(ldh*2)+2]);
    __SIMD_DATATYPE h_4_3 = _SSE_SET1(hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SSE_SET1(hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SSE_SET1(hh[(ldh*3)+3]);
#endif

#if VEC_SET == 1281
    __SIMD_DATATYPE h_2_1 = _SSE_SET(hh[ldh+1], hh[ldh+1]);
    __SIMD_DATATYPE h_3_2 = _SSE_SET(hh[(ldh*2)+1], hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SSE_SET(hh[(ldh*2)+2], hh[(ldh*2)+2]);
    __SIMD_DATATYPE h_4_3 = _SSE_SET(hh[(ldh*3)+1], hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SSE_SET(hh[(ldh*3)+2], hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SSE_SET(hh[(ldh*3)+3], hh[(ldh*3)+3]);
#endif

#if VEC_SET == 256
    __SIMD_DATATYPE h_2_1 = _SIMD_BROADCAST(&hh[ldh+1]);
    __SIMD_DATATYPE h_3_2 = _SIMD_BROADCAST(&hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SIMD_BROADCAST(&hh[(ldh*2)+2]);
    __SIMD_DATATYPE h_4_3 = _SIMD_BROADCAST(&hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SIMD_BROADCAST(&hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SIMD_BROADCAST(&hh[(ldh*3)+3]);
#endif

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w1 = _SIMD_FMA(a3_1, h_4_3, a4_1);
    w1 = _SIMD_FMA(a2_1, h_4_2, w1);
    w1 = _SIMD_FMA(a1_1, h_4_1, w1);
    register __SIMD_DATATYPE z1 = _SIMD_FMA(a2_1, h_3_2, a3_1);
    z1 = _SIMD_FMA(a1_1, h_3_1, z1);
    register __SIMD_DATATYPE y1 = _SIMD_FMA(a1_1, h_2_1, a2_1);
    register __SIMD_DATATYPE x1 = a1_1;
#else
    register __SIMD_DATATYPE w1 = _SIMD_ADD(a4_1, _SIMD_MUL(a3_1, h_4_3));
    w1 = _SIMD_ADD(w1, _SIMD_MUL(a2_1, h_4_2));                          
    w1 = _SIMD_ADD(w1, _SIMD_MUL(a1_1, h_4_1));                          
    register __SIMD_DATATYPE z1 = _SIMD_ADD(a3_1, _SIMD_MUL(a2_1, h_3_2));
    z1 = _SIMD_ADD(z1, _SIMD_MUL(a1_1, h_3_1));                          
    register __SIMD_DATATYPE y1 = _SIMD_ADD(a2_1, _SIMD_MUL(a1_1, h_2_1));
    register __SIMD_DATATYPE x1 = a1_1;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE a1_2 = _SIMD_LOAD(&q[(ldq*3)+offset]);                  
    __SIMD_DATATYPE a2_2 = _SIMD_LOAD(&q[(ldq*2)+offset]);
    __SIMD_DATATYPE a3_2 = _SIMD_LOAD(&q[ldq+offset]);
    __SIMD_DATATYPE a4_2 = _SIMD_LOAD(&q[0+offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w2 = _SIMD_FMA(a3_2, h_4_3, a4_2);
    w2 = _SIMD_FMA(a2_2, h_4_2, w2);
    w2 = _SIMD_FMA(a1_2, h_4_1, w2);
    register __SIMD_DATATYPE z2 = _SIMD_FMA(a2_2, h_3_2, a3_2);
    z2 = _SIMD_FMA(a1_2, h_3_1, z2);
    register __SIMD_DATATYPE y2 = _SIMD_FMA(a1_2, h_2_1, a2_2);
    register __SIMD_DATATYPE x2 = a1_2;
#else
    register __SIMD_DATATYPE w2 = _SIMD_ADD(a4_2, _SIMD_MUL(a3_2, h_4_3));
    w2 = _SIMD_ADD(w2, _SIMD_MUL(a2_2, h_4_2));
    w2 = _SIMD_ADD(w2, _SIMD_MUL(a1_2, h_4_1));
    register __SIMD_DATATYPE z2 = _SIMD_ADD(a3_2, _SIMD_MUL(a2_2, h_3_2));
    z2 = _SIMD_ADD(z2, _SIMD_MUL(a1_2, h_3_1));
    register __SIMD_DATATYPE y2 = _SIMD_ADD(a2_2, _SIMD_MUL(a1_2, h_2_1));
    register __SIMD_DATATYPE x2 = a1_2;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE a1_3 = _SIMD_LOAD(&q[(ldq*3)+2*offset]);
    __SIMD_DATATYPE a2_3 = _SIMD_LOAD(&q[(ldq*2)+2*offset]);
    __SIMD_DATATYPE a3_3 = _SIMD_LOAD(&q[ldq+2*offset]);
    __SIMD_DATATYPE a4_3 = _SIMD_LOAD(&q[0+2*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w3 = _SIMD_FMA(a3_3, h_4_3, a4_3);
    w3 = _SIMD_FMA(a2_3, h_4_2, w3);
    w3 = _SIMD_FMA(a1_3, h_4_1, w3);
    register __SIMD_DATATYPE z3 = _SIMD_FMA(a2_3, h_3_2, a3_3);
    z3 = _SIMD_FMA(a1_3, h_3_1, z3);
    register __SIMD_DATATYPE y3 = _SIMD_FMA(a1_3, h_2_1, a2_3);
    register __SIMD_DATATYPE x3 = a1_3;
#else
    register __SIMD_DATATYPE w3 = _SIMD_ADD(a4_3, _SIMD_MUL(a3_3, h_4_3));
    w3 = _SIMD_ADD(w3, _SIMD_MUL(a2_3, h_4_2));
    w3 = _SIMD_ADD(w3, _SIMD_MUL(a1_3, h_4_1));
    register __SIMD_DATATYPE z3 = _SIMD_ADD(a3_3, _SIMD_MUL(a2_3, h_3_2));
    z3 = _SIMD_ADD(z3, _SIMD_MUL(a1_3, h_3_1));
    register __SIMD_DATATYPE y3 = _SIMD_ADD(a2_3, _SIMD_MUL(a1_3, h_2_1));
    register __SIMD_DATATYPE x3 = a1_3;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE a1_4 = _SIMD_LOAD(&q[(ldq*3)+3*offset]);
    __SIMD_DATATYPE a2_4 = _SIMD_LOAD(&q[(ldq*2)+3*offset]);
    __SIMD_DATATYPE a3_4 = _SIMD_LOAD(&q[ldq+3*offset]);
    __SIMD_DATATYPE a4_4 = _SIMD_LOAD(&q[0+3*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w4 = _SIMD_FMA(a3_4, h_4_3, a4_4);
    w4 = _SIMD_FMA(a2_4, h_4_2, w4);
    w4 = _SIMD_FMA(a1_4, h_4_1, w4);
    register __SIMD_DATATYPE z4 = _SIMD_FMA(a2_4, h_3_2, a3_4);
    z4 = _SIMD_FMA(a1_4, h_3_1, z4);
    register __SIMD_DATATYPE y4 = _SIMD_FMA(a1_4, h_2_1, a2_4);
    register __SIMD_DATATYPE x4 = a1_4;
#else
    register __SIMD_DATATYPE w4 = _SIMD_ADD(a4_4, _SIMD_MUL(a3_4, h_4_3));
    w4 = _SIMD_ADD(w4, _SIMD_MUL(a2_4, h_4_2));
    w4 = _SIMD_ADD(w4, _SIMD_MUL(a1_4, h_4_1));
    register __SIMD_DATATYPE z4 = _SIMD_ADD(a3_4, _SIMD_MUL(a2_4, h_3_2));
    z4 = _SIMD_ADD(z4, _SIMD_MUL(a1_4, h_3_1));
    register __SIMD_DATATYPE y4 = _SIMD_ADD(a2_4, _SIMD_MUL(a1_4, h_2_1));
    register __SIMD_DATATYPE x4 = a1_4;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE a1_5 = _SIMD_LOAD(&q[(ldq*3)+4*offset]);
    __SIMD_DATATYPE a2_5 = _SIMD_LOAD(&q[(ldq*2)+4*offset]);
    __SIMD_DATATYPE a3_5 = _SIMD_LOAD(&q[ldq+4*offset]);
    __SIMD_DATATYPE a4_5 = _SIMD_LOAD(&q[0+4*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w5 = _SIMD_FMA(a3_5, h_4_3, a4_5);
    w5 = _SIMD_FMA(a2_5, h_4_2, w5);
    w5 = _SIMD_FMA(a1_5, h_4_1, w5);
    register __SIMD_DATATYPE z5 = _SIMD_FMA(a2_5, h_3_2, a3_5);
    z5 = _SIMD_FMA(a1_5, h_3_1, z5);
    register __SIMD_DATATYPE y5 = _SIMD_FMA(a1_5, h_2_1, a2_5);
    register __SIMD_DATATYPE x5 = a1_5;
#else
    register __SIMD_DATATYPE w5 = _SIMD_ADD(a4_5, _SIMD_MUL(a3_5, h_4_3));
    w5 = _SIMD_ADD(w5, _SIMD_MUL(a2_5, h_4_2));
    w5 = _SIMD_ADD(w5, _SIMD_MUL(a1_5, h_4_1));
    register __SIMD_DATATYPE z5 = _SIMD_ADD(a3_5, _SIMD_MUL(a2_5, h_3_2));
    z5 = _SIMD_ADD(z5, _SIMD_MUL(a1_5, h_3_1));
    register __SIMD_DATATYPE y5 = _SIMD_ADD(a2_5, _SIMD_MUL(a1_5, h_2_1));
    register __SIMD_DATATYPE x5 = a1_5;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE a1_6 = _SIMD_LOAD(&q[(ldq*3)+5*offset]);
    __SIMD_DATATYPE a2_6 = _SIMD_LOAD(&q[(ldq*2)+5*offset]);
    __SIMD_DATATYPE a3_6 = _SIMD_LOAD(&q[ldq+5*offset]);
    __SIMD_DATATYPE a4_6 = _SIMD_LOAD(&q[0+5*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w6 = _SIMD_FMA(a3_6, h_4_3, a4_6);
    w6 = _SIMD_FMA(a2_6, h_4_2, w6);
    w6 = _SIMD_FMA(a1_6, h_4_1, w6);
    register __SIMD_DATATYPE z6 = _SIMD_FMA(a2_6, h_3_2, a3_6);
    z6 = _SIMD_FMA(a1_6, h_3_1, z6);
    register __SIMD_DATATYPE y6 = _SIMD_FMA(a1_6, h_2_1, a2_6);
    register __SIMD_DATATYPE x6 = a1_6;
#else
    register __SIMD_DATATYPE w6 = _SIMD_ADD(a4_6, _SIMD_MUL(a3_6, h_4_3));
    w6 = _SIMD_ADD(w6, _SIMD_MUL(a2_6, h_4_2));
    w6 = _SIMD_ADD(w6, _SIMD_MUL(a1_6, h_4_1));
    register __SIMD_DATATYPE z6 = _SIMD_ADD(a3_6, _SIMD_MUL(a2_6, h_3_2));
    z6 = _SIMD_ADD(z6, _SIMD_MUL(a1_6, h_3_1));
    register __SIMD_DATATYPE y6 = _SIMD_ADD(a2_6, _SIMD_MUL(a1_6, h_2_1));
    register __SIMD_DATATYPE x6 = a1_6;
#endif /* __ELPA_USE_FMA__ */

    __SIMD_DATATYPE q1;
    __SIMD_DATATYPE q2;
    __SIMD_DATATYPE q3;
    __SIMD_DATATYPE q4;
    __SIMD_DATATYPE q5;
    __SIMD_DATATYPE q6;

    __SIMD_DATATYPE h1;
    __SIMD_DATATYPE h2;
    __SIMD_DATATYPE h3;
    __SIMD_DATATYPE h4;
#endif /* BLOCK4 */

#ifdef BLOCK6
    
    __SIMD_DATATYPE a1_1 = _SIMD_LOAD(&q[ldq*5]);
    __SIMD_DATATYPE a2_1 = _SIMD_LOAD(&q[ldq*4]);
    __SIMD_DATATYPE a3_1 = _SIMD_LOAD(&q[ldq*3]);
    __SIMD_DATATYPE a4_1 = _SIMD_LOAD(&q[ldq*2]);
    __SIMD_DATATYPE a5_1 = _SIMD_LOAD(&q[ldq]);  
    __SIMD_DATATYPE a6_1 = _SIMD_LOAD(&q[0]);    

#if VEC_SET == 128
    __SIMD_DATATYPE h_6_5 = _SSE_SET1(hh[(ldh*5)+1]);
    __SIMD_DATATYPE h_6_4 = _SSE_SET1(hh[(ldh*5)+2]);
    __SIMD_DATATYPE h_6_3 = _SSE_SET1(hh[(ldh*5)+3]);
    __SIMD_DATATYPE h_6_2 = _SSE_SET1(hh[(ldh*5)+4]);
    __SIMD_DATATYPE h_6_1 = _SSE_SET1(hh[(ldh*5)+5]);
#endif

#if VEC_SET == 1281
    __SIMD_DATATYPE h_6_5 = _SSE_SET(hh[(ldh*5)+1], hh[(ldh*5)+1]);
    __SIMD_DATATYPE h_6_4 = _SSE_SET(hh[(ldh*5)+2], hh[(ldh*5)+2]);
    __SIMD_DATATYPE h_6_3 = _SSE_SET(hh[(ldh*5)+3], hh[(ldh*5)+3]);
    __SIMD_DATATYPE h_6_2 = _SSE_SET(hh[(ldh*5)+4], hh[(ldh*5)+4]);
    __SIMD_DATATYPE h_6_1 = _SSE_SET(hh[(ldh*5)+5], hh[(ldh*5)+5]);
#endif

#if VEC_SET == 256
    __SIMD_DATATYPE h_6_5 = _SIMD_BROADCAST(&hh[(ldh*5)+1]);
    __SIMD_DATATYPE h_6_4 = _SIMD_BROADCAST(&hh[(ldh*5)+2]);
    __SIMD_DATATYPE h_6_3 = _SIMD_BROADCAST(&hh[(ldh*5)+3]);
    __SIMD_DATATYPE h_6_2 = _SIMD_BROADCAST(&hh[(ldh*5)+4]);
    __SIMD_DATATYPE h_6_1 = _SIMD_BROADCAST(&hh[(ldh*5)+5]);
#endif

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t1 = _SIMD_FMA(a5_1, h_6_5, a6_1);
    t1 = _SIMD_FMA(a4_1, h_6_4, t1);
    t1 = _SIMD_FMA(a3_1, h_6_3, t1);
    t1 = _SIMD_FMA(a2_1, h_6_2, t1);
    t1 = _SIMD_FMA(a1_1, h_6_1, t1);
#else
    register __SIMD_DATATYPE t1 = _SIMD_ADD(a6_1, _SIMD_MUL(a5_1, h_6_5)); 
    t1 = _SIMD_ADD(t1, _SIMD_MUL(a4_1, h_6_4));
    t1 = _SIMD_ADD(t1, _SIMD_MUL(a3_1, h_6_3));
    t1 = _SIMD_ADD(t1, _SIMD_MUL(a2_1, h_6_2));
    t1 = _SIMD_ADD(t1, _SIMD_MUL(a1_1, h_6_1));
#endif /* __ELPA_USE_FMA__ */

#if VEC_SET == 128
    __SIMD_DATATYPE h_5_4 = _SSE_SET1(hh[(ldh*4)+1]);
    __SIMD_DATATYPE h_5_3 = _SSE_SET1(hh[(ldh*4)+2]);
    __SIMD_DATATYPE h_5_2 = _SSE_SET1(hh[(ldh*4)+3]);
    __SIMD_DATATYPE h_5_1 = _SSE_SET1(hh[(ldh*4)+4]);
#endif

#if VEC_SET == 1281
    __SIMD_DATATYPE h_5_4 = _SSE_SET(hh[(ldh*4)+1], hh[(ldh*4)+1]);
    __SIMD_DATATYPE h_5_3 = _SSE_SET(hh[(ldh*4)+2], hh[(ldh*4)+2]);
    __SIMD_DATATYPE h_5_2 = _SSE_SET(hh[(ldh*4)+3], hh[(ldh*4)+3]);
    __SIMD_DATATYPE h_5_1 = _SSE_SET(hh[(ldh*4)+4], hh[(ldh*4)+4]);
#endif

#if VEC_SET == 256
    __SIMD_DATATYPE h_5_4 = _SIMD_BROADCAST(&hh[(ldh*4)+1]);
    __SIMD_DATATYPE h_5_3 = _SIMD_BROADCAST(&hh[(ldh*4)+2]);
    __SIMD_DATATYPE h_5_2 = _SIMD_BROADCAST(&hh[(ldh*4)+3]);
    __SIMD_DATATYPE h_5_1 = _SIMD_BROADCAST(&hh[(ldh*4)+4]);
#endif

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE v1 = _SIMD_FMA(a4_1, h_5_4, a5_1);
    v1 = _SIMD_FMA(a3_1, h_5_3, v1);
    v1 = _SIMD_FMA(a2_1, h_5_2, v1);
    v1 = _SIMD_FMA(a1_1, h_5_1, v1);
#else
    register __SIMD_DATATYPE v1 = _SIMD_ADD(a5_1, _SIMD_MUL(a4_1, h_5_4)); 
    v1 = _SIMD_ADD(v1, _SIMD_MUL(a3_1, h_5_3));
    v1 = _SIMD_ADD(v1, _SIMD_MUL(a2_1, h_5_2));
    v1 = _SIMD_ADD(v1, _SIMD_MUL(a1_1, h_5_1));
#endif /* __ELPA_USE_FMA__ */

#if VEC_SET == 128
    __SIMD_DATATYPE h_4_3 = _SSE_SET1(hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SSE_SET1(hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SSE_SET1(hh[(ldh*3)+3]);
#endif

#if VEC_SET == 1281
    __SIMD_DATATYPE h_4_3 = _SSE_SET(hh[(ldh*3)+1], hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SSE_SET(hh[(ldh*3)+2], hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SSE_SET(hh[(ldh*3)+3], hh[(ldh*3)+3]);
#endif

#if VEC_SET == 256
    __SIMD_DATATYPE h_4_3 = _SIMD_BROADCAST(&hh[(ldh*3)+1]);
    __SIMD_DATATYPE h_4_2 = _SIMD_BROADCAST(&hh[(ldh*3)+2]);
    __SIMD_DATATYPE h_4_1 = _SIMD_BROADCAST(&hh[(ldh*3)+3]);
#endif

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE w1 = _SIMD_FMA(a3_1, h_4_3, a4_1);
    w1 = _SIMD_FMA(a2_1, h_4_2, w1);
    w1 = _SIMD_FMA(a1_1, h_4_1, w1);
#else
    register __SIMD_DATATYPE w1 = _SIMD_ADD(a4_1, _SIMD_MUL(a3_1, h_4_3)); 
    w1 = _SIMD_ADD(w1, _SIMD_MUL(a2_1, h_4_2));
    w1 = _SIMD_ADD(w1, _SIMD_MUL(a1_1, h_4_1));
#endif /* __ELPA_USE_FMA__ */

#if VEC_SET == 128
    __SIMD_DATATYPE h_2_1 = _SSE_SET1(hh[ldh+1]);    
    __SIMD_DATATYPE h_3_2 = _SSE_SET1(hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SSE_SET1(hh[(ldh*2)+2]);
#endif

#if VEC_SET == 1281
    __SIMD_DATATYPE h_2_1 = _SSE_SET(hh[ldh+1], hh[ldh+1]);
    __SIMD_DATATYPE h_3_2 = _SSE_SET(hh[(ldh*2)+1], hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SSE_SET(hh[(ldh*2)+2], hh[(ldh*2)+2]);
#endif

#if VEC_SET == 256
    __SIMD_DATATYPE h_2_1 = _SIMD_BROADCAST(&hh[ldh+1]);
    __SIMD_DATATYPE h_3_2 = _SIMD_BROADCAST(&hh[(ldh*2)+1]);
    __SIMD_DATATYPE h_3_1 = _SIMD_BROADCAST(&hh[(ldh*2)+2]);
#endif

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE z1 = _SIMD_FMA(a2_1, h_3_2, a3_1);
    z1 = _SIMD_FMA(a1_1, h_3_1, z1);
    register __SIMD_DATATYPE y1 = _SIMD_FMA(a1_1, h_2_1, a2_1);
#else
    register __SIMD_DATATYPE z1 = _SIMD_ADD(a3_1, _SIMD_MUL(a2_1, h_3_2));
    z1 = _SIMD_ADD(z1, _SIMD_MUL(a1_1, h_3_1));
    register __SIMD_DATATYPE y1 = _SIMD_ADD(a2_1, _SIMD_MUL(a1_1, h_2_1)); 
#endif /* __ELPA_USE_FMA__ */

    register __SIMD_DATATYPE x1 = a1_1;

    __SIMD_DATATYPE a1_2 = _SIMD_LOAD(&q[(ldq*5)+offset]);
    __SIMD_DATATYPE a2_2 = _SIMD_LOAD(&q[(ldq*4)+offset]);
    __SIMD_DATATYPE a3_2 = _SIMD_LOAD(&q[(ldq*3)+offset]);
    __SIMD_DATATYPE a4_2 = _SIMD_LOAD(&q[(ldq*2)+offset]);
    __SIMD_DATATYPE a5_2 = _SIMD_LOAD(&q[(ldq)+offset]);
    __SIMD_DATATYPE a6_2 = _SIMD_LOAD(&q[offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t2 = _SIMD_FMA(a5_2, h_6_5, a6_2);
    t2 = _SIMD_FMA(a4_2, h_6_4, t2);
    t2 = _SIMD_FMA(a3_2, h_6_3, t2);
    t2 = _SIMD_FMA(a2_2, h_6_2, t2);
    t2 = _SIMD_FMA(a1_2, h_6_1, t2);
    register __SIMD_DATATYPE v2 = _SIMD_FMA(a4_2, h_5_4, a5_2);
    v2 = _SIMD_FMA(a3_2, h_5_3, v2);
    v2 = _SIMD_FMA(a2_2, h_5_2, v2);
    v2 = _SIMD_FMA(a1_2, h_5_1, v2);
    register __SIMD_DATATYPE w2 = _SIMD_FMA(a3_2, h_4_3, a4_2);
    w2 = _SIMD_FMA(a2_2, h_4_2, w2);
    w2 = _SIMD_FMA(a1_2, h_4_1, w2);
    register __SIMD_DATATYPE z2 = _SIMD_FMA(a2_2, h_3_2, a3_2);
    z2 = _SIMD_FMA(a1_2, h_3_1, z2);
    register __SIMD_DATATYPE y2 = _SIMD_FMA(a1_2, h_2_1, a2_2);
#else
    register __SIMD_DATATYPE t2 = _SIMD_ADD(a6_2, _SIMD_MUL(a5_2, h_6_5));
    t2 = _SIMD_ADD(t2, _SIMD_MUL(a4_2, h_6_4));
    t2 = _SIMD_ADD(t2, _SIMD_MUL(a3_2, h_6_3));
    t2 = _SIMD_ADD(t2, _SIMD_MUL(a2_2, h_6_2));
    t2 = _SIMD_ADD(t2, _SIMD_MUL(a1_2, h_6_1));
    register __SIMD_DATATYPE v2 = _SIMD_ADD(a5_2, _SIMD_MUL(a4_2, h_5_4));
    v2 = _SIMD_ADD(v2, _SIMD_MUL(a3_2, h_5_3));
    v2 = _SIMD_ADD(v2, _SIMD_MUL(a2_2, h_5_2));
    v2 = _SIMD_ADD(v2, _SIMD_MUL(a1_2, h_5_1));
    register __SIMD_DATATYPE w2 = _SIMD_ADD(a4_2, _SIMD_MUL(a3_2, h_4_3));
    w2 = _SIMD_ADD(w2, _SIMD_MUL(a2_2, h_4_2));
    w2 = _SIMD_ADD(w2, _SIMD_MUL(a1_2, h_4_1));
    register __SIMD_DATATYPE z2 = _SIMD_ADD(a3_2, _SIMD_MUL(a2_2, h_3_2));
    z2 = _SIMD_ADD(z2, _SIMD_MUL(a1_2, h_3_1));
    register __SIMD_DATATYPE y2 = _SIMD_ADD(a2_2, _SIMD_MUL(a1_2, h_2_1));
#endif /* __ELPA_USE_FMA__ */

    register __SIMD_DATATYPE x2 = a1_2;

    __SIMD_DATATYPE a1_3 = _SIMD_LOAD(&q[(ldq*5)+2*offset]);
    __SIMD_DATATYPE a2_3 = _SIMD_LOAD(&q[(ldq*4)+2*offset]);
    __SIMD_DATATYPE a3_3 = _SIMD_LOAD(&q[(ldq*3)+2*offset]);
    __SIMD_DATATYPE a4_3 = _SIMD_LOAD(&q[(ldq*2)+2*offset]);
    __SIMD_DATATYPE a5_3 = _SIMD_LOAD(&q[(ldq)+2*offset]);
    __SIMD_DATATYPE a6_3 = _SIMD_LOAD(&q[2*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t3 = _SIMD_FMA(a5_3, h_6_5, a6_3);
    t3 = _SIMD_FMA(a4_3, h_6_4, t3);
    t3 = _SIMD_FMA(a3_3, h_6_3, t3);
    t3 = _SIMD_FMA(a2_3, h_6_2, t3);
    t3 = _SIMD_FMA(a1_3, h_6_1, t3);
    register __SIMD_DATATYPE v3 = _SIMD_FMA(a4_3, h_5_4, a5_3);
    v3 = _SIMD_FMA(a3_3, h_5_3, v3);
    v3 = _SIMD_FMA(a2_3, h_5_2, v3);
    v3 = _SIMD_FMA(a1_3, h_5_1, v3);
    register __SIMD_DATATYPE w3 = _SIMD_FMA(a3_3, h_4_3, a4_3);
    w3 = _SIMD_FMA(a2_3, h_4_2, w3);
    w3 = _SIMD_FMA(a1_3, h_4_1, w3);
    register __SIMD_DATATYPE z3 = _SIMD_FMA(a2_3, h_3_2, a3_3);
    z3 = _SIMD_FMA(a1_3, h_3_1, z3);
    register __SIMD_DATATYPE y3 = _SIMD_FMA(a1_3, h_2_1, a2_3);
#else
    register __SIMD_DATATYPE t3 = _SIMD_ADD(a6_3, _SIMD_MUL(a5_3, h_6_5));
    t3 = _SIMD_ADD(t3, _SIMD_MUL(a4_3, h_6_4));
    t3 = _SIMD_ADD(t3, _SIMD_MUL(a3_3, h_6_3));
    t3 = _SIMD_ADD(t3, _SIMD_MUL(a2_3, h_6_2));
    t3 = _SIMD_ADD(t3, _SIMD_MUL(a1_3, h_6_1));
    register __SIMD_DATATYPE v3 = _SIMD_ADD(a5_3, _SIMD_MUL(a4_3, h_5_4));
    v3 = _SIMD_ADD(v3, _SIMD_MUL(a3_3, h_5_3));
    v3 = _SIMD_ADD(v3, _SIMD_MUL(a2_3, h_5_2));
    v3 = _SIMD_ADD(v3, _SIMD_MUL(a1_3, h_5_1));
    register __SIMD_DATATYPE w3 = _SIMD_ADD(a4_3, _SIMD_MUL(a3_3, h_4_3));
    w3 = _SIMD_ADD(w3, _SIMD_MUL(a2_3, h_4_2));
    w3 = _SIMD_ADD(w3, _SIMD_MUL(a1_3, h_4_1));
    register __SIMD_DATATYPE z3 = _SIMD_ADD(a3_3, _SIMD_MUL(a2_3, h_3_2));
    z3 = _SIMD_ADD(z3, _SIMD_MUL(a1_3, h_3_1));
    register __SIMD_DATATYPE y3 = _SIMD_ADD(a2_3, _SIMD_MUL(a1_3, h_2_1));
#endif /* __ELPA_USE_FMA__ */

    register __SIMD_DATATYPE x3 = a1_3;

    __SIMD_DATATYPE a1_4 = _SIMD_LOAD(&q[(ldq*5)+3*offset]);
    __SIMD_DATATYPE a2_4 = _SIMD_LOAD(&q[(ldq*4)+3*offset]);
    __SIMD_DATATYPE a3_4 = _SIMD_LOAD(&q[(ldq*3)+3*offset]);
    __SIMD_DATATYPE a4_4 = _SIMD_LOAD(&q[(ldq*2)+3*offset]);
    __SIMD_DATATYPE a5_4 = _SIMD_LOAD(&q[(ldq)+3*offset]);
    __SIMD_DATATYPE a6_4 = _SIMD_LOAD(&q[3*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t4 = _SIMD_FMA(a5_4, h_6_5, a6_4);
    t4 = _SIMD_FMA(a4_4, h_6_4, t4);
    t4 = _SIMD_FMA(a3_4, h_6_3, t4);
    t4 = _SIMD_FMA(a2_4, h_6_2, t4);
    t4 = _SIMD_FMA(a1_4, h_6_1, t4);
    register __SIMD_DATATYPE v4 = _SIMD_FMA(a4_4, h_5_4, a5_4);
    v4 = _SIMD_FMA(a3_4, h_5_3, v4);
    v4 = _SIMD_FMA(a2_4, h_5_2, v4);
    v4 = _SIMD_FMA(a1_4, h_5_1, v4);
    register __SIMD_DATATYPE w4 = _SIMD_FMA(a3_4, h_4_3, a4_4);
    w4 = _SIMD_FMA(a2_4, h_4_2, w4);
    w4 = _SIMD_FMA(a1_4, h_4_1, w4);
    register __SIMD_DATATYPE z4 = _SIMD_FMA(a2_4, h_3_2, a3_4);
    z4 = _SIMD_FMA(a1_4, h_3_1, z4);
    register __SIMD_DATATYPE y4 = _SIMD_FMA(a1_4, h_2_1, a2_4);
#else
    register __SIMD_DATATYPE t4 = _SIMD_ADD(a6_4, _SIMD_MUL(a5_4, h_6_5));
    t4 = _SIMD_ADD(t4, _SIMD_MUL(a4_4, h_6_4));
    t4 = _SIMD_ADD(t4, _SIMD_MUL(a3_4, h_6_3));
    t4 = _SIMD_ADD(t4, _SIMD_MUL(a2_4, h_6_2));
    t4 = _SIMD_ADD(t4, _SIMD_MUL(a1_4, h_6_1));
    register __SIMD_DATATYPE v4 = _SIMD_ADD(a5_4, _SIMD_MUL(a4_4, h_5_4));
    v4 = _SIMD_ADD(v4, _SIMD_MUL(a3_4, h_5_3));
    v4 = _SIMD_ADD(v4, _SIMD_MUL(a2_4, h_5_2));
    v4 = _SIMD_ADD(v4, _SIMD_MUL(a1_4, h_5_1));
    register __SIMD_DATATYPE w4 = _SIMD_ADD(a4_4, _SIMD_MUL(a3_4, h_4_3));
    w4 = _SIMD_ADD(w4, _SIMD_MUL(a2_4, h_4_2));
    w4 = _SIMD_ADD(w4, _SIMD_MUL(a1_4, h_4_1));
    register __SIMD_DATATYPE z4 = _SIMD_ADD(a3_4, _SIMD_MUL(a2_4, h_3_2));
    z4 = _SIMD_ADD(z4, _SIMD_MUL(a1_4, h_3_1));
    register __SIMD_DATATYPE y4 = _SIMD_ADD(a2_4, _SIMD_MUL(a1_4, h_2_1));
#endif /* __ELPA_USE_FMA__ */

    register __SIMD_DATATYPE x4 = a1_4;

    __SIMD_DATATYPE a1_5 = _SIMD_LOAD(&q[(ldq*5)+4*offset]);
    __SIMD_DATATYPE a2_5 = _SIMD_LOAD(&q[(ldq*4)+4*offset]);
    __SIMD_DATATYPE a3_5 = _SIMD_LOAD(&q[(ldq*3)+4*offset]);
    __SIMD_DATATYPE a4_5 = _SIMD_LOAD(&q[(ldq*2)+4*offset]);
    __SIMD_DATATYPE a5_5 = _SIMD_LOAD(&q[(ldq)+4*offset]);
    __SIMD_DATATYPE a6_5 = _SIMD_LOAD(&q[4*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t5 = _SIMD_FMA(a5_5, h_6_5, a6_5);
    t5 = _SIMD_FMA(a4_5, h_6_4, t5);
    t5 = _SIMD_FMA(a3_5, h_6_3, t5);
    t5 = _SIMD_FMA(a2_5, h_6_2, t5);
    t5 = _SIMD_FMA(a1_5, h_6_1, t5);
    register __SIMD_DATATYPE v5 = _SIMD_FMA(a4_5, h_5_4, a5_5);
    v5 = _SIMD_FMA(a3_5, h_5_3, v5);
    v5 = _SIMD_FMA(a2_5, h_5_2, v5);
    v5 = _SIMD_FMA(a1_5, h_5_1, v5);
    register __SIMD_DATATYPE w5 = _SIMD_FMA(a3_5, h_4_3, a4_5);
    w5 = _SIMD_FMA(a2_5, h_4_2, w5);
    w5 = _SIMD_FMA(a1_5, h_4_1, w5);
    register __SIMD_DATATYPE z5 = _SIMD_FMA(a2_5, h_3_2, a3_5);
    z5 = _SIMD_FMA(a1_5, h_3_1, z5);
    register __SIMD_DATATYPE y5 = _SIMD_FMA(a1_5, h_2_1, a2_5);
#else
    register __SIMD_DATATYPE t5 = _SIMD_ADD(a6_5, _SIMD_MUL(a5_5, h_6_5));
    t5 = _SIMD_ADD(t5, _SIMD_MUL(a4_5, h_6_4));
    t5 = _SIMD_ADD(t5, _SIMD_MUL(a3_5, h_6_3));
    t5 = _SIMD_ADD(t5, _SIMD_MUL(a2_5, h_6_2));
    t5 = _SIMD_ADD(t5, _SIMD_MUL(a1_5, h_6_1));
    register __SIMD_DATATYPE v5 = _SIMD_ADD(a5_5, _SIMD_MUL(a4_5, h_5_4));
    v5 = _SIMD_ADD(v5, _SIMD_MUL(a3_5, h_5_3));
    v5 = _SIMD_ADD(v5, _SIMD_MUL(a2_5, h_5_2));
    v5 = _SIMD_ADD(v5, _SIMD_MUL(a1_5, h_5_1));
    register __SIMD_DATATYPE w5 = _SIMD_ADD(a4_5, _SIMD_MUL(a3_5, h_4_3));
    w5 = _SIMD_ADD(w5, _SIMD_MUL(a2_5, h_4_2));
    w5 = _SIMD_ADD(w5, _SIMD_MUL(a1_5, h_4_1));
    register __SIMD_DATATYPE z5 = _SIMD_ADD(a3_5, _SIMD_MUL(a2_5, h_3_2));
    z5 = _SIMD_ADD(z5, _SIMD_MUL(a1_5, h_3_1));
    register __SIMD_DATATYPE y5 = _SIMD_ADD(a2_5, _SIMD_MUL(a1_5, h_2_1));
#endif /* __ELPA_USE_FMA__ */

    register __SIMD_DATATYPE x5 = a1_5;

    __SIMD_DATATYPE a1_6 = _SIMD_LOAD(&q[(ldq*5)+5*offset]);
    __SIMD_DATATYPE a2_6 = _SIMD_LOAD(&q[(ldq*4)+5*offset]);
    __SIMD_DATATYPE a3_6 = _SIMD_LOAD(&q[(ldq*3)+5*offset]);
    __SIMD_DATATYPE a4_6 = _SIMD_LOAD(&q[(ldq*2)+5*offset]);
    __SIMD_DATATYPE a5_6 = _SIMD_LOAD(&q[(ldq)+5*offset]);
    __SIMD_DATATYPE a6_6 = _SIMD_LOAD(&q[5*offset]);

#ifdef __ELPA_USE_FMA__
    register __SIMD_DATATYPE t6 = _SIMD_FMA(a5_6, h_6_5, a6_6);
    t6 = _SIMD_FMA(a4_6, h_6_4, t6);
    t6 = _SIMD_FMA(a3_6, h_6_3, t6);
    t6 = _SIMD_FMA(a2_6, h_6_2, t6);
    t6 = _SIMD_FMA(a1_6, h_6_1, t6);
    register __SIMD_DATATYPE v6 = _SIMD_FMA(a4_6, h_5_4, a5_6);
    v6 = _SIMD_FMA(a3_6, h_5_3, v6);
    v6 = _SIMD_FMA(a2_6, h_5_2, v6);
    v6 = _SIMD_FMA(a1_6, h_5_1, v6);
    register __SIMD_DATATYPE w6 = _SIMD_FMA(a3_6, h_4_3, a4_6);
    w6 = _SIMD_FMA(a2_6, h_4_2, w6);
    w6 = _SIMD_FMA(a1_6, h_4_1, w6);
    register __SIMD_DATATYPE z6 = _SIMD_FMA(a2_6, h_3_2, a3_6);
    z6 = _SIMD_FMA(a1_6, h_3_1, z6);
    register __SIMD_DATATYPE y6 = _SIMD_FMA(a1_6, h_2_1, a2_6);
#else
    register __SIMD_DATATYPE t6 = _SIMD_ADD(a6_6, _SIMD_MUL(a5_6, h_6_5));
    t6 = _SIMD_ADD(t6, _SIMD_MUL(a4_6, h_6_4));
    t6 = _SIMD_ADD(t6, _SIMD_MUL(a3_6, h_6_3));
    t6 = _SIMD_ADD(t6, _SIMD_MUL(a2_6, h_6_2));
    t6 = _SIMD_ADD(t6, _SIMD_MUL(a1_6, h_6_1));
    register __SIMD_DATATYPE v6 = _SIMD_ADD(a5_6, _SIMD_MUL(a4_6, h_5_4));
    v6 = _SIMD_ADD(v6, _SIMD_MUL(a3_6, h_5_3));
    v6 = _SIMD_ADD(v6, _SIMD_MUL(a2_6, h_5_2));
    v6 = _SIMD_ADD(v6, _SIMD_MUL(a1_6, h_5_1));
    register __SIMD_DATATYPE w6 = _SIMD_ADD(a4_6, _SIMD_MUL(a3_6, h_4_3));
    w6 = _SIMD_ADD(w6, _SIMD_MUL(a2_6, h_4_2));
    w6 = _SIMD_ADD(w6, _SIMD_MUL(a1_6, h_4_1));
    register __SIMD_DATATYPE z6 = _SIMD_ADD(a3_6, _SIMD_MUL(a2_6, h_3_2));
    z6 = _SIMD_ADD(z6, _SIMD_MUL(a1_6, h_3_1));
    register __SIMD_DATATYPE y6 = _SIMD_ADD(a2_6, _SIMD_MUL(a1_6, h_2_1));
#endif /* __ELPA_USE_FMA__ */
 
    register __SIMD_DATATYPE x6 = a1_6;

    __SIMD_DATATYPE q1;
    __SIMD_DATATYPE q2;
    __SIMD_DATATYPE q3;
    __SIMD_DATATYPE q4;
    __SIMD_DATATYPE q5;
    __SIMD_DATATYPE q6;

    __SIMD_DATATYPE h1;
    __SIMD_DATATYPE h2;
    __SIMD_DATATYPE h3;
    __SIMD_DATATYPE h4;
    __SIMD_DATATYPE h5;
    __SIMD_DATATYPE h6;

#endif /* BLOCK6 */


    for(i = BLOCK; i < nb; i++)
      {

#if VEC_SET == 128
        h1 = _SSE_SET1(hh[i-(BLOCK-1)]);
        h2 = _SSE_SET1(hh[ldh+i-(BLOCK-2)]);
#endif
#if VEC_SET == 1281
        h1 = _SSE_SET(hh[i-(BLOCK-1)], hh[i-(BLOCK-1)]);
        h2 = _SSE_SET(hh[ldh+i-(BLOCK-2)], hh[ldh+i-(BLOCK-2)]);
#endif
#if  VEC_SET == 256
        h1 = _SIMD_BROADCAST(&hh[i-(BLOCK-1)]);
        h2 = _SIMD_BROADCAST(&hh[ldh+i-(BLOCK-2)]);
#endif /*   VEC_SET == 256 */

#ifdef __ELPA_USE_FMA__
        q1 = _SIMD_LOAD(&q[i*ldq]);
        x1 = _SIMD_FMA(q1, h1, x1);
        y1 = _SIMD_FMA(q1, h2, y1);
        q2 = _SIMD_LOAD(&q[(i*ldq)+offset]);
        x2 = _SIMD_FMA(q2, h1, x2);
        y2 = _SIMD_FMA(q2, h2, y2);
        q3 = _SIMD_LOAD(&q[(i*ldq)+2*offset]);
        x3 = _SIMD_FMA(q3, h1, x3);
        y3 = _SIMD_FMA(q3, h2, y3);
        q4 = _SIMD_LOAD(&q[(i*ldq)+3*offset]);
        x4 = _SIMD_FMA(q4, h1, x4);
        y4 = _SIMD_FMA(q4, h2, y4);
        q5 = _SIMD_LOAD(&q[(i*ldq)+4*offset]);
        x5 = _SIMD_FMA(q5, h1, x5);
        y5 = _SIMD_FMA(q5, h2, y5);
        q6 = _SIMD_LOAD(&q[(i*ldq)+5*offset]);
        x6 = _SIMD_FMA(q6, h1, x6);
        y6 = _SIMD_FMA(q6, h2, y6);
#else
        q1 = _SIMD_LOAD(&q[i*ldq]);
        x1 = _SIMD_ADD(x1, _SIMD_MUL(q1,h1));
        y1 = _SIMD_ADD(y1, _SIMD_MUL(q1,h2));
        q2 = _SIMD_LOAD(&q[(i*ldq)+offset]);
        x2 = _SIMD_ADD(x2, _SIMD_MUL(q2,h1));
        y2 = _SIMD_ADD(y2, _SIMD_MUL(q2,h2));
        q3 = _SIMD_LOAD(&q[(i*ldq)+2*offset]);
        x3 = _SIMD_ADD(x3, _SIMD_MUL(q3,h1));
        y3 = _SIMD_ADD(y3, _SIMD_MUL(q3,h2));
        q4 = _SIMD_LOAD(&q[(i*ldq)+3*offset]);
        x4 = _SIMD_ADD(x4, _SIMD_MUL(q4,h1));
        y4 = _SIMD_ADD(y4, _SIMD_MUL(q4,h2));
        q5 = _SIMD_LOAD(&q[(i*ldq)+4*offset]);
        x5 = _SIMD_ADD(x5, _SIMD_MUL(q5,h1));
        y5 = _SIMD_ADD(y5, _SIMD_MUL(q5,h2));
        q6 = _SIMD_LOAD(&q[(i*ldq)+5*offset]);
        x6 = _SIMD_ADD(x6, _SIMD_MUL(q6,h1));
        y6 = _SIMD_ADD(y6, _SIMD_MUL(q6,h2));
#endif /* __ELPA_USE_FMA__ */

#if defined(BLOCK4) || defined(BLOCK6)
#if VEC_SET == 128
        h3 = _SSE_SET1(hh[(ldh*2)+i-(BLOCK-3)]);
#endif

#if VEC_SET == 1281