/*
prodmvEy.c
matrix mode
*/

#include "oth.h"

void prodmvEy(d_complex_t a[], d_complex_t v[], d_complex_t u[])
{
	const int64_t i0 = i0Ey;
	const int64_t i1 = i1Ey;
	const int64_t j0 = j0Ey;
	const int64_t j1 = j1Ey;
	const int64_t k0 = k0Ey;
	const int64_t k1 = k1Ey;

	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (        i = i0; i <= i1; i++) {
	for (int64_t j = j0; j <= j1; j++) {
		int64_t n = NEY(i, j, k0);  // NEYK=1
		d_complex_t *aptr = &a[n * NNZ];
		d_complex_t *v00p = &v[NEY(i,     j,     k0    )];
		d_complex_t *v01p = &v[NEY(i,     j,     k0 + 1)];
		d_complex_t *v02p = &v[NEY(i,     j,     k0 - 1)];
		d_complex_t *v03p = &v[NEY(i + 1, j,     k0    )];
		d_complex_t *v04p = &v[NEY(i - 1, j,     k0    )];
		d_complex_t *v05p = &v[NEZ(i,     j + 1, k0    )];
		d_complex_t *v06p = &v[NEZ(i,     j,     k0    )];
		d_complex_t *v07p = &v[NEZ(i,     j + 1, k0 - 1)];
		d_complex_t *v08p = &v[NEZ(i,     j,     k0 - 1)];
		d_complex_t *v09p = &v[NEX(i,     j + 1, k0    )];
		d_complex_t *v10p = &v[NEX(i,     j,     k0    )];
		d_complex_t *v11p = &v[NEX(i - 1, j + 1, k0    )];
		d_complex_t *v12p = &v[NEX(i - 1, j,     k0    )];
#if defined(__FUJITSU)
#pragma loop novrec
#elif defined(__CLANG_FUJITSU)
#pragma clang loop vectorize(assume_safety)
#endif
	for (int64_t k = k0; k <= k1; k++, n++) {
		const d_complex_t a00 = *aptr++;
		const d_complex_t a01 = *aptr++;
		const d_complex_t a02 = *aptr++;
		const d_complex_t a03 = *aptr++;
		const d_complex_t a04 = *aptr++;
		const d_complex_t a05 = *aptr++;
		const d_complex_t a06 = *aptr++;
		const d_complex_t a07 = *aptr++;
		const d_complex_t a08 = *aptr++;
		const d_complex_t a09 = *aptr++;
		const d_complex_t a10 = *aptr++;
		const d_complex_t a11 = *aptr++;
		const d_complex_t a12 = *aptr++;

		const d_complex_t v00 = *v00p++;
		const d_complex_t v01 = *v01p++;
		const d_complex_t v02 = *v02p++;
		const d_complex_t v03 = *v03p++;
		const d_complex_t v04 = *v04p++;
		const d_complex_t v05 = *v05p++;
		const d_complex_t v06 = *v06p++;
		const d_complex_t v07 = *v07p++;
		const d_complex_t v08 = *v08p++;
		const d_complex_t v09 = *v09p++;
		const d_complex_t v10 = *v10p++;
		const d_complex_t v11 = *v11p++;
		const d_complex_t v12 = *v12p++;

		u[n].r =
			(a00.r * v00.r) - (a00.i * v00.i) +
			(a01.r * v01.r) - (a01.i * v01.i) +
			(a02.r * v02.r) - (a02.i * v02.i) +
			(a03.r * v03.r) - (a03.i * v03.i) +
			(a04.r * v04.r) - (a04.i * v04.i) +
			(a05.r * v05.r) - (a05.i * v05.i) +
			(a06.r * v06.r) - (a06.i * v06.i) +
			(a07.r * v07.r) - (a07.i * v07.i) +
			(a08.r * v08.r) - (a08.i * v08.i) +
			(a09.r * v09.r) - (a09.i * v09.i) +
			(a10.r * v10.r) - (a10.i * v10.i) +
			(a11.r * v11.r) - (a11.i * v11.i) +
			(a12.r * v12.r) - (a12.i * v12.i);

		u[n].i =
			(a00.r * v00.i) + (a00.i * v00.r) +
			(a01.r * v01.i) + (a01.i * v01.r) +
			(a02.r * v02.i) + (a02.i * v02.r) +
			(a03.r * v03.i) + (a03.i * v03.r) +
			(a04.r * v04.i) + (a04.i * v04.r) +
			(a05.r * v05.i) + (a05.i * v05.r) +
			(a06.r * v06.i) + (a06.i * v06.r) +
			(a07.r * v07.i) + (a07.i * v07.r) +
			(a08.r * v08.i) + (a08.i * v08.r) +
			(a09.r * v09.i) + (a09.i * v09.r) +
			(a10.r * v10.i) + (a10.i * v10.r) +
			(a11.r * v11.i) + (a11.i * v11.r) +
			(a12.r * v12.i) + (a12.i * v12.r);
	}
	}
	}
}
