#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>

#define blockSize 32
#define f0 0.0
#define f1 0.0

__global__ void kernel(double *u, double *uOld, double a, double b, double h, double tau, double r, int N)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;

	if (i == 0)
	{
		u[0] = f0;
	}
	else if (i == N - 1)
	{
		u[N - 1] = f1;
	}
	else
	{
		u[i] = r * (i * h) * (i * h) * ((u[i + 1] - u[i]) - (u[i] - u[i - 1])) - (uOld[i] - 2 * u[i]);
		uOld[i] = u[i];
	}
}

double g(const double x)
{
	return 0.5*x;
}

double uol(const double x)
{
	return 0.9*x;
}


int main()
{
	double a = 1.0, L = 10.0, x = 0.0, time = 0.0, tmax = 10.0, tau = 0.001, h = 0.1, H = 1.0, b = 1.0;
	int N = (int)(L / h) + 1;
	double r = a * tau * tau / (h * h);

	double *uDev = NULL, *uOldDev = NULL;
	int size = N * sizeof(double);

	double *u = new double[N];
	double *uOld = new double[N];
	for (int i = 1; i < N - 1; i++)
	{
		u[i] = g(i*h);
		uOld[i] = u[i] - tau * uol(i*h);
	}
	u[0] = uOld[0] = f0;
	u[N - 1] = uOld[N - 1] = f1;

	float calculationTime = 0.0f;
	cudaEvent_t tn, tk;
	cudaEventCreate(&tn);
	cudaEventCreate(&tk);
	cudaEventRecord(tn, 0);

	cudaMalloc((void**)&uDev, size);
	cudaMalloc((void**)&uOldDev, size);

	cudaMemcpy(uDev, u, size, cudaMemcpyHostToDevice);
	cudaMemcpy(uOldDev, uOld, size, cudaMemcpyHostToDevice);

	do {
		kernel << < dim3(N / blockSize + 1), dim3(blockSize) >> > (uDev, uOldDev, a, b, h, tau, r, N);
		cudaDeviceSynchronize();

		cudaMemcpy(uOldDev, uDev, size, cudaMemcpyDeviceToDevice);

		time += tau;
	} while (time <= tmax);

	cudaMemcpy(u, uDev, size, cudaMemcpyDeviceToHost);

	cudaEventRecord(tk, 0);
	cudaEventSynchronize(tk);
	cudaEventElapsedTime(&calculationTime, tn, tk);

	cudaError_t err = cudaGetLastError();
	const char* error = cudaGetErrorString(err);
	printf("Cuda error = %s\n", error);

	if (err == cudaSuccess)
	{
		FILE* f = fopen("result.txt", "w");

		for (int i = 0; i < N; i++)
		{
			fprintf(f, "u[%d] = %f\n", i, u[i]);
		}

		fprintf(f, "\nTime = %f", calculationTime / 1000.0);
		fclose(f);
	}

	return 0;
}