/*
matrixEx.c
u = Av (Ex)
*/

#include "oth.h"
#include "complex13.h"

void matrixEx(int pre, int ifreq, d_complex_t a[], const d_complex_t v[], d_complex_t u[], d_complex_t b[])
{
	const int m0 = ifreq * NMaterial;

	int64_t i0 = i0Ex;
	int64_t i1 = i1Ex;
	int64_t j0 = j0Ex;
	int64_t j1 = j1Ex;
	int64_t k0 = k0Ex;
	int64_t k1 = k1Ex;

	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (        i = i0; i <= i1; i++) {
	for (int64_t j = j0; j <= j1; j++) {
#if defined(__NEC__)
#pragma _NEC ivdep
#elif defined(__FUJITSU)
#pragma loop novrec
#elif defined(__CLANG_FUJITSU)
#pragma clang loop vectorize(assume_safety)
#endif
	for (int64_t k = k0; k <= k1; k++) {
		d_complex_t c[NNZ];
		d_set13(c, d_complex(0, 0));

		if (MEX(i, j, k) == PEC) {
			c[0] = d_complex(1, 0);
		}
		else {
			if ((Abc == MUR) && (
				(!Pbc.ypbc && ((j == 0) || (j == Ny))) ||
				(!Pbc.zpbc && ((k == 0) || (k == Nz))))) {
				// Mur
				c[0] = d_complex(1, 0);
				d_complex_t cv = Material[MEX(i, j, k) + m0].cv;
				d_complex_t cel = d_complex(0, 0);
				if      (!Pbc.ypbc && (j ==  0)) {
					cel = d_rmul(KDYn[ 0], cv);
					c[1] = d_neg(d_rmul(exp(+cel.i), d_exp(-cel.r)));
				}
				else if (!Pbc.ypbc && (j == Ny)) {
					cel = d_rmul(KDYn[Ny], cv);
					c[2] = d_neg(d_rmul(exp(+cel.i), d_exp(-cel.r)));
				}
				else if (!Pbc.zpbc && (k ==  0)) {
					cel = d_rmul(KDZn[ 0], cv);
					c[3] = d_neg(d_rmul(exp(+cel.i), d_exp(-cel.r)));
				}
				else if (!Pbc.zpbc && (k == Nz)) {
					cel = d_rmul(KDZn[Nz], cv);
					c[4] = d_neg(d_rmul(exp(+cel.i), d_exp(-cel.r)));
				}
			}
			else if ((Abc == PML) &&
				((i < 0) || (i >= Nx) || (j <= 0) || (j >= Ny) || (k <= 0) || (k >= Nz))) {
				// PML
				int64_t ib = MINMAX(0, Nx - 1, i);
				int64_t jb = MINMAX(0, Ny,     j);
				int64_t kb = MINMAX(0, Nz,     k);
				if (MEX(ib, jb, kb) != PEC) {
					d_complex_t vc = d_inv(Material[MEX(ib, jb, kb) + m0].cv);

					d_complex_t yc = d_div(vc, d_complex(Fypml[(2 * j + 0) + (2 * Lpml)], KDYn[jb]));
					d_complex_t yp = d_div(vc, d_complex(Fypml[(2 * j + 1) + (2 * Lpml)], KDYn[jb]));
					d_complex_t ym = d_div(vc, d_complex(Fypml[(2 * j - 1) + (2 * Lpml)], KDYn[jb]));
					d_complex_t zc = d_div(vc, d_complex(Fzpml[(2 * k + 0) + (2 * Lpml)], KDZn[kb]));
					d_complex_t zp = d_div(vc, d_complex(Fzpml[(2 * k + 1) + (2 * Lpml)], KDZn[kb]));
					d_complex_t zm = d_div(vc, d_complex(Fzpml[(2 * k - 1) + (2 * Lpml)], KDZn[kb]));
					d_complex_t xp = d_div(vc, d_complex(Fxpml[(2 * i + 1) + (2 * Lpml)], KDXc[ib]));

					c[ 0] = d_complex(1, 0);
					c[ 1] = d_neg(d_mul(yc, yp));
					c[ 2] = d_neg(d_mul(yc, ym));
					c[ 3] = d_neg(d_mul(zc, zp));
					c[ 4] = d_neg(d_mul(zc, zm));
					c[ 5] =       d_mul(yc, xp);
					c[ 6] = d_neg(c[ 5]);
					c[ 7] =       c[ 6];
					c[ 8] =       c[ 5];
					c[ 9] =       d_mul(zc, xp);
					c[10] = d_neg(c[ 9]);
					c[11] =       c[10];
					c[12] =       c[ 9];
					c[ 0] = d_sub(c[ 0], d_add4(c[ 1], c[ 2], c[ 3], c[ 4]));
				}
				else {
					c[0] = d_complex(1, 0);
				}
			}
			else {
				d_complex_t c0 = Material[MEX(i,     j,     k    ) + m0].cepsr;
				d_complex_t c1 = Material[MHZ(i,     j,     k    ) + m0].cnur;
				d_complex_t c2 = Material[MHZ(i,     j - 1, k    ) + m0].cnur;
				d_complex_t c3 = Material[MHY(i,     j,     k    ) + m0].cnur;
				d_complex_t c4 = Material[MHY(i,     j,     k - 1) + m0].cnur;

				c[ 0] = d_neg(c0);
				c[ 1] = d_rmul(- KRYp[j],           c1);
				c[ 2] = d_rmul(- KRYm[j],           c2);
				c[ 3] = d_rmul(- KRZp[k],           c3);
				c[ 4] = d_rmul(- KRZm[k],           c4);
				c[ 5] = d_rmul(+ KRXc[i] * KRYn[j], c1);
				c[ 6] = d_neg(c[ 5]);
				c[ 7] = d_rmul(- KRXc[i] * KRYn[j], c2);
				c[ 8] = d_neg(c[ 7]);
				c[ 9] = d_rmul(+ KRXc[i] * KRZn[k], c3);
				c[10] = d_neg(c[ 9]);
				c[11] = d_rmul(- KRXc[i] * KRZn[k], c4);
				c[12] = d_neg(c[11]);
				c[ 0] = d_sub(c[ 0], d_add4(c[ 1], c[ 2], c[ 3], c[ 4]));
			}
		}

		int64_t n = NEX(i, j, k);
		d_complex_t sf = SCALING ? d_inv(c[0]) : d_complex(1, 0);
		if (pre) {
			// setup matrix
			if (MATRIX) {
				d_scal13(&a[n * NNZ], c, sf);
			}
			// scaling RHS
			b[n] = d_mul(sf, b[n]);
		}
		else {
			// matrix * vector
			d_complex_t sum = d_add13(
				d_mul(c[ 0], v[NEX(i,     j,     k    )]),
				d_mul(c[ 1], v[NEX(i,     j + 1, k    )]),
				d_mul(c[ 2], v[NEX(i,     j - 1, k    )]),
				d_mul(c[ 3], v[NEX(i,     j,     k + 1)]),
				d_mul(c[ 4], v[NEX(i,     j,     k - 1)]),
				d_mul(c[ 5], v[NEY(i + 1, j,     k    )]),
				d_mul(c[ 6], v[NEY(i,     j,     k    )]),
				d_mul(c[ 7], v[NEY(i + 1, j - 1, k    )]),
				d_mul(c[ 8], v[NEY(i,     j - 1, k    )]),
				d_mul(c[ 9], v[NEZ(i + 1, j,     k    )]),
				d_mul(c[10], v[NEZ(i,     j,     k    )]),
				d_mul(c[11], v[NEZ(i + 1, j,     k - 1)]),
				d_mul(c[12], v[NEZ(i,     j,     k - 1)]));
			u[n] = d_mul(sf, sum);
		}
	}
	}
	}
}
