/*
blas1c.cu (CUDA)

BLAS functions, level-1, double complex
grid : 1D
block : 1D
*/

#include "oth.h"
#include "oth_cuda.h"
#include "reduction_sum.cu"

extern void cuda_memcpy(int, void *, const void *, size_t, cudaMemcpyKind);

// [y] = [x]
__global__
static void Zcopy_kernel(int64_t n, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	int64_t tid = threadIdx.x + (blockIdx.x * blockDim.x);

	if (tid < n) {
		y[tid].x = x[tid].x;
		y[tid].y = x[tid].y;
	}
}

void Zcopy(int64_t n, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	if (GPU) {
		unsigned int grid = (unsigned int)CEIL(n, BLAS_BLOCK_SIZE);
		Zcopy_kernel<<<grid, BLAS_BLOCK_SIZE>>>(n, x, y);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		for (int64_t i = 0; i < n; i++) {
			y[i].x = x[i].x;
			y[i].y = x[i].y;
		}
	}
}
/*
// [x] * [y]
__global__
static void Zdotu_kernel(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y, cuDoubleComplex *work)
{
	extern __shared__ double s[];

	int64_t tid = (blockIdx.x * blockDim.x) + threadIdx.x;

	double sum_r = 0;
	double sum_i = 0;
	while (tid < n) {
		sum_r += (x[tid].x * y[tid].x) - (x[tid].y * y[tid].y);
		sum_i += (x[tid].x * y[tid].y) + (x[tid].y * y[tid].x);
		tid += blockDim.x * gridDim.x;
	}
	s[0          + threadIdx.x] = sum_r;
	s[blockDim.x + threadIdx.x] = sum_i;

	reduction_sum(threadIdx.x, blockDim.x, &s[0         ], &work[blockIdx.x].x);
	reduction_sum(threadIdx.x, blockDim.x, &s[blockDim.x], &work[blockIdx.x].y);
}

cuDoubleComplex Zdotu(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y)
{
	cuDoubleComplex ret = make_cuDoubleComplex(0, 0);

	if (GPU) {
		size_t sm_size = 2 * BLAS_BLOCK_SIZE * sizeof(double);
		Zdotu_kernel<<<BLAS_GRID_SIZE, BLAS_BLOCK_SIZE, sm_size>>>(n, x, y, d_Wblas);
		if (UM) cudaDeviceSynchronize();
		if (!UM) {
			cudaMemcpy(h_Wblas, d_Wblas, BLAS_GRID_SIZE * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost);
			for (int i = 0; i < BLAS_GRID_SIZE; i++) {
				ret = cuCadd(ret, h_Wblas[i]);
			}
		}
		else {
			for (int i = 0; i < BLAS_GRID_SIZE; i++) {
				ret = cuCadd(ret, d_Wblas[i]);
			}
		}
	}
	else {
		for (int64_t i = 0; i < n; i++) {
			ret.x += (x[i].x * y[i].x) - (x[i].y * y[i].y);
			ret.y += (x[i].x * y[i].y) + (x[i].y * y[i].x);
		}
	}

	return ret;
}
*/
// conj[x] * [y]
__global__
static void Zdotc_kernel(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y, cuDoubleComplex *work)
{
	extern __shared__ double s[];

	int64_t tid = (blockIdx.x * blockDim.x) + threadIdx.x;

	double sum_r = 0;
	double sum_i = 0;
	while (tid < n) {
		sum_r += (x[tid].x * y[tid].x) + (x[tid].y * y[tid].y);
		sum_i += (x[tid].x * y[tid].y) - (x[tid].y * y[tid].x);
		tid += blockDim.x * gridDim.x;
	}
	s[0          + threadIdx.x] = sum_r;
	s[blockDim.x + threadIdx.x] = sum_i;

	reduction_sum(threadIdx.x, blockDim.x, &s[0         ], &work[blockIdx.x].x);
	reduction_sum(threadIdx.x, blockDim.x, &s[blockDim.x], &work[blockIdx.x].y);
}

cuDoubleComplex Zdotc(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y)
{
	cuDoubleComplex ret = make_cuDoubleComplex(0, 0);

	if (GPU) {
		size_t sm_size = 2 * BLAS_BLOCK_SIZE * sizeof(double);
		Zdotc_kernel<<<BLAS_GRID_SIZE, BLAS_BLOCK_SIZE, sm_size>>>(n, x, y, d_Wblas);
		if (UM) cudaDeviceSynchronize();
		if (!UM) {
			cudaMemcpy(h_Wblas, d_Wblas, BLAS_GRID_SIZE * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost);
			for (int i = 0; i < BLAS_GRID_SIZE; i++) {
				ret = cuCadd(ret, h_Wblas[i]);
			}
		}
		else {
			for (int i = 0; i < BLAS_GRID_SIZE; i++) {
				ret = cuCadd(ret, d_Wblas[i]);
			}
		}
	}
	else {
		for (int64_t i = 0; i < n; i++) {
			ret.x += (x[i].x * y[i].x) + (x[i].y * y[i].y);
			ret.y += (x[i].x * y[i].y) - (x[i].y * y[i].x);
		}
	}

	return ret;
}

// y[i] = a * x[i] + y[i]
__global__
static void Zaxpy_kernel(int64_t n, cuDoubleComplex a, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	int64_t tid = threadIdx.x + (blockIdx.x * blockDim.x);

	if (tid < n) {
		y[tid].x += (a.x * x[tid].x) - (a.y * x[tid].y);
		y[tid].y += (a.x * x[tid].y) + (a.y * x[tid].x);
	}
}

void Zaxpy(int64_t n, cuDoubleComplex a, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	if (GPU) {
		unsigned int grid = (unsigned int)CEIL(n, BLAS_BLOCK_SIZE);
		Zaxpy_kernel<<<grid, BLAS_BLOCK_SIZE>>>(n, a, x, y);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		for (int64_t i = 0; i < n; i++) {
			y[i].x += (a.x * x[i].x) - (a.y * x[i].y);
			y[i].y += (a.x * x[i].y) + (a.y * x[i].x);
		}
	}
}

// |x|
double Dznrm2(int64_t n, const cuDoubleComplex *x)
{
	return sqrt(Zdotc(n, x, x).x);
}

// y[i] = a * x[i] + y[i]
__global__
static void Zaxbycz_kernel(int64_t n, cuDoubleComplex a, cuDoubleComplex b, cuDoubleComplex c, cuDoubleComplex *x, const cuDoubleComplex *y, const cuDoubleComplex *z)
{
	int64_t tid = threadIdx.x + (blockIdx.x * blockDim.x);

	if (tid < n) {
		double x_x = (a.x * x[tid].x) - (a.y * x[tid].y) + (b.x * y[tid].x) - (b.y * y[tid].y) + (c.x * z[tid].x) - (c.y * z[tid].y);
		double x_y = (a.x * x[tid].y) + (a.y * x[tid].x) + (b.x * y[tid].y) + (b.y * y[tid].x) + (c.x * z[tid].y) + (c.y * z[tid].x);
		x[tid].x = x_x;
		x[tid].y = x_y;
	}
}

// x[i] = a * x[i] + b * y[i] + c * z[i] (custom)
void Zaxbycz(int64_t n, cuDoubleComplex a, cuDoubleComplex b, cuDoubleComplex c, cuDoubleComplex *x, const cuDoubleComplex *y, const cuDoubleComplex *z)
{
	if (GPU) {
		unsigned int grid = (unsigned int)CEIL(n, BLAS_BLOCK_SIZE);
		Zaxbycz_kernel<<<grid, BLAS_BLOCK_SIZE>>>(n, a, b, c, x, y, z);
		if (UM) cudaDeviceSynchronize();
	}
	else {
		for (int64_t i = 0; i < n; i++) {
			double x_x = (a.x * x[i].x) - (a.y * x[i].y) + (b.x * y[i].x) - (b.y * y[i].y) + (c.x * z[i].x) - (c.y * z[i].y);
			double x_y = (a.x * x[i].y) + (a.y * x[i].x) + (b.x * y[i].y) + (b.y * y[i].x) + (c.x * z[i].y) + (c.y * z[i].x);
			x[i].x = x_x;
			x[i].y = x_y;
		}
	}
}
