/*
pbcY.cu

copy vector (Y-PBC)
*/

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

// copy
__host__ __device__ __forceinline__
static void copyYpbcEy_(cuDoubleComplex *v, int k, int i, param_t *q)
{
	v[nEY(q, i,       - 1, k)] = v[nEY(q, i, q->Ny - 1, k)];
	v[nEY(q, i, q->Ny + 0, k)] = v[nEY(q, i,         0, k)];
}
__host__ __device__ __forceinline__
static void copyYpbcEz_(cuDoubleComplex *v, int k, int i, param_t *q)
{
	v[nEZ(q, i,       - 1, k)] = v[nEZ(q, i, q->Ny - 1, k)];
	v[nEZ(q, i, q->Ny + 1, k)] = v[nEZ(q, i,         1, k)];
}
__host__ __device__ __forceinline__
static void copyYpbcEx_(cuDoubleComplex *v, int k, int i, param_t *q)
{
	v[nEX(q, i,       - 1, k)] = v[nEX(q, i, q->Ny - 1, k)];
	v[nEX(q, i, q->Ny + 1, k)] = v[nEX(q, i,         1, k)];
}

// factor
__host__ __device__ __forceinline__
static void factorYpbcEy_(cuDoubleComplex *v, int j0, int j1, int k, int i, param_t *q)
{
	int64_t n;
	int j;
	const cuDoubleComplex fp = q->pbc_fy;
	const cuDoubleComplex fm = cuCinv(fp);
	j =       - 1;
	n = nEY(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fm);
	}
	j = q->Ny + 0;
	n = nEY(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fp);
	}
}
__host__ __device__ __forceinline__
static void factorYpbcEz_(cuDoubleComplex *v, int j0, int j1, int k, int i, param_t *q)
{
	int64_t n;
	int j;
	const cuDoubleComplex fp = q->pbc_fy;
	const cuDoubleComplex fm = cuCinv(fp);
	j =       - 1;
	n = nEZ(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fm);
	}
	j = q->Ny + 1;
	n = nEZ(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fp);
	}
}
__host__ __device__ __forceinline__
static void factorYpbcEx_(cuDoubleComplex *v, int j0, int j1, int k, int i, param_t *q)
{
	int64_t n;
	int j;
	const cuDoubleComplex fp = q->pbc_fy;
	const cuDoubleComplex fm = cuCinv(fp);
	j =       - 1;
	n = nEX(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fm);
	}
	j = q->Ny + 1;
	n = nEX(q, i, j, k);
	if ((j >= j0 - 1) && (j <= j1 + 1)) {  // MPI
		v[n] = cuCmul(v[n], fp);
	}
}

// (GPU)
__global__
static void copyYpbc_gpu(char c, cuDoubleComplex *v, int k0, int k1, int i0, int i1)
{
	int k = (k0 - 1) + threadIdx.x + (blockIdx.x * blockDim.x);
	int i = (i0 - 1) + threadIdx.y + (blockIdx.y * blockDim.y);
	if ((k <= k1 + 1) &&
	    (i <= i1 + 1)) {
		if (c == 'X') copyYpbcEx_(v, k, i, &d_Param);
		if (c == 'Y') copyYpbcEy_(v, k, i, &d_Param);
		if (c == 'Z') copyYpbcEz_(v, k, i, &d_Param);
	}
}
__global__
static void factorYpbc_gpu(char c, cuDoubleComplex *v, int i0, int i1, int j0, int j1, int k0, int k1)
{
	int k = (k0 - 1) + threadIdx.x + (blockIdx.x * blockDim.x);
	int i = (i0 - 1) + threadIdx.y + (blockIdx.y * blockDim.y);
	if ((k <= k1 + 1) &&
	    (i <= i1 + 1)) {
		if (c == 'X') factorYpbcEx_(v, j0, j1, k, i, &d_Param);
		if (c == 'Y') factorYpbcEy_(v, j0, j1, k, i, &d_Param);
		if (c == 'Z') factorYpbcEz_(v, j0, j1, k, i, &d_Param);
	}
}

// (CPU)
static void copyYpbc_cpu(char c, cuDoubleComplex *v, int k0, int k1, int i0, int i1)
{
	for (int k = k0 - 1; k <= k1 + 1; k++) {
	for (int i = i0 - 1; i <= i1 + 1; i++) {
		if (c == 'X') copyYpbcEx_(v, k, i, &h_Param);
		if (c == 'Y') copyYpbcEy_(v, k, i, &h_Param);
		if (c == 'Z') copyYpbcEz_(v, k, i, &h_Param);
	}
	}
}
static void factorYpbc_cpu(char c, cuDoubleComplex *v, int i0, int i1, int j0, int j1, int k0, int k1)
{
	for (int k = k0 - 1; k <= k1 + 1; k++) {
	for (int i = i0 - 1; i <= i1 + 1; i++) {
		if (c == 'X') factorYpbcEx_(v, j0, j1, k, i, &h_Param);
		if (c == 'Y') factorYpbcEy_(v, j0, j1, k, i, &h_Param);
		if (c == 'Z') factorYpbcEz_(v, j0, j1, k, i, &h_Param);
	}
	}
}

void copyYpbc(cuDoubleComplex *v)
{
	assert(Pbc.ypbc);

	if (GPU) {
		dim3 block(PBC_BLOCK_X, PBC_BLOCK_Y);
		dim3 gridx(CEIL(k1Ex - k0Ex + 3, block.x),
		           CEIL(i1Ex - i0Ex + 3, block.y));
		dim3 gridy(CEIL(k1Ey - k0Ey + 3, block.x),
		           CEIL(i1Ey - i0Ey + 3, block.y));
		dim3 gridz(CEIL(k1Ez - k0Ez + 3, block.x),
		           CEIL(i1Ez - i0Ez + 3, block.y));
		cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
		copyYpbc_gpu<<<gridx, block>>>('X', v, k0Ex, k1Ex, i0Ex, i1Ex);
		copyYpbc_gpu<<<gridy, block>>>('Y', v, k0Ey, k1Ey, i0Ey, i1Ey);
		copyYpbc_gpu<<<gridz, block>>>('Z', v, k0Ez, k1Ez, i0Ez, i1Ez);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		copyYpbc_cpu('X', v, k0Ex, k1Ex, i0Ex, i1Ex);
		copyYpbc_cpu('Y', v, k0Ey, k1Ey, i0Ey, i1Ey);
		copyYpbc_cpu('Z', v, k0Ez, k1Ez, i0Ez, i1Ez);
	}
}

void factorYpbc(cuDoubleComplex *v)
{
	assert(Pbc.ypbc);

	if (GPU) {
		dim3 block(PBC_BLOCK_X, PBC_BLOCK_Y);
		dim3 gridx(CEIL(k1Ex - k0Ex + 3, block.x),
		           CEIL(i1Ex - i0Ex + 3, block.y));
		dim3 gridy(CEIL(k1Ey - k0Ey + 3, block.x),
		           CEIL(i1Ey - i0Ey + 3, block.y));
		dim3 gridz(CEIL(k1Ez - k0Ez + 3, block.x),
		           CEIL(i1Ez - i0Ez + 3, block.y));
		cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
		factorYpbc_gpu<<<gridx, block>>>('X', v, i0Ex, i1Ex, j0Ex, j1Ex, k0Ex, k1Ex);
		factorYpbc_gpu<<<gridy, block>>>('Y', v, i0Ey, i1Ey, j0Ey, j1Ey, k0Ey, k1Ey);
		factorYpbc_gpu<<<gridz, block>>>('Z', v, i0Ez, i1Ez, j0Ez, j1Ez, k0Ez, k1Ez);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		factorYpbc_cpu('X', v, i0Ex, i1Ex, j0Ex, j1Ex, k0Ex, k1Ex);
		factorYpbc_cpu('Y', v, i0Ey, i1Ey, j0Ey, j1Ey, k0Ey, k1Ey);
		factorYpbc_cpu('Z', v, i0Ez, i1Ez, j0Ez, j1Ez, k0Ez, k1Ez);
	}
}
