/*
prodmvEy.cu
u = Av (Ey)
*/

#include "oth.h"
#include "oth_cuda.h"

__host__ __device__ __forceinline__
static void prodmvEy_(
	int i, int j, int k,
	const cuDoubleComplex *a, const cuDoubleComplex *v, cuDoubleComplex *u,
	param_t *q)
{
	int64_t n = nEY(q, i, j, k);
	int64_t p[NNZ];
	p[ 0] = nEY(q, i,     j,     k    );
	p[ 1] = nEY(q, i,     j,     k + 1);
	p[ 2] = nEY(q, i,     j,     k - 1);
	p[ 3] = nEY(q, i + 1, j,     k    );
	p[ 4] = nEY(q, i - 1, j,     k    );
	p[ 5] = nEZ(q, i,     j + 1, k    );
	p[ 6] = nEZ(q, i,     j,     k    );
	p[ 7] = nEZ(q, i,     j + 1, k - 1);
	p[ 8] = nEZ(q, i,     j,     k - 1);
	p[ 9] = nEX(q, i,     j + 1, k    );
	p[10] = nEX(q, i,     j,     k    );
	p[11] = nEX(q, i - 1, j + 1, k    );
	p[12] = nEX(q, i - 1, j,     k    );

	cuDoubleComplex sum = make_cuDoubleComplex(0, 0);
	for (int m = 0; m < NNZ; m++) {
		sum = cuCadd(sum, cuCmul(a[(n * NNZ) + m], v[p[m]]));
	}
	u[n] = sum;
}

__global__
static void prodmvEy_gpu(
	int i0, int i1, int j0, int j1, int k0, int k1,
	const cuDoubleComplex *a, const cuDoubleComplex *v, cuDoubleComplex *u)
{
	int i = i0 + threadIdx.z + (blockIdx.z * blockDim.z);
	int j = j0 + threadIdx.y + (blockIdx.y * blockDim.y);
	int k = k0 + threadIdx.x + (blockIdx.x * blockDim.x);
	if ((i <= i1) &&
	    (j <= j1) &&
	    (k <= k1)) {
		prodmvEy_(
			i, j, k,
			a, v, u,
			&d_Param);
	}
}

static void prodmvEy_cpu(
	int i0, int i1, int j0, int j1, int k0, int k1,
	const cuDoubleComplex *a, const cuDoubleComplex *v, cuDoubleComplex *u)
{
	for (int i = i0; i <= i1; i++) {
	for (int j = j0; j <= j1; j++) {
	for (int k = k0; k <= k1; k++) {
		prodmvEy_(
			i, j, k,
			a, v, u,
			&h_Param);
	}
	}
	}
}

void prodmvEy(const cuDoubleComplex *a, const cuDoubleComplex *v, cuDoubleComplex *u)
{
	if (GPU) {
		cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
		dim3 block(BLOCK_X, BLOCK_Y, 1);
		dim3 grid(
			CEIL(k1Ey - k0Ey + 1, block.x),
			CEIL(j1Ey - j0Ey + 1, block.y),
			CEIL(i1Ey - i0Ey + 1, block.z));
		prodmvEy_gpu<<<grid, block>>>(
			i0Ey, i1Ey, j0Ey, j1Ey, k0Ey, k1Ey,
			a, v, u);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		prodmvEy_cpu(
			i0Ey, i1Ey, j0Ey, j1Ey, k0Ey, k1Ey,
			a, v, u);
	}
}
