/*
solve.cu (CUDA)
*/

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

void solve(FILE *fp, int iout)
{
	const size_t size = NE * sizeof(d_complex_t);

	// alloc host memory
	memalloc();
	Bcg = (d_complex_t *)malloc(size);
	Ecg = (d_complex_t *)malloc(size);
	Pcg = (d_complex_t *)malloc(size);

	// alloc device memory
	memalloc_gpu();
	memalloc_cg_gpu();
	memcopy_gpu();

	// alloc work array
	size_t work_size = BLAS_GRID_SIZE * sizeof(cuDoubleComplex);
	h_Wblas =     (cuDoubleComplex *)malloc(work_size);
	cuda_malloc(GPU, UM, (void **)&d_Wblas, work_size);

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		// monitor
		if (iout) {
			monitor5(fp, ifreq);
		}

		// setup
		setup2(ifreq);
		memcopy2_gpu();
		setup_parameter_gpu();

		// setup RHS vector
		memset(Bcg, 0, size);
		rhs(Bcg, ifreq);
		cuda_memcpy(GPU, d_Bcg, Bcg, size, cudaMemcpyHostToDevice);

		// (setup matrix and) scaling RHS vector
		if (MATRIX) cuda_memset(GPU, d_Acg, 0, NNZ * size);
		prodmv(1, ifreq, d_Acg, NULL, NULL, d_Bcg);

		// clear
		cuda_memset(GPU, d_Rcg, 0, size);
		cuda_memset(GPU, d_Pcg, 0, size);
		cuda_memset(GPU, d_Qcg, 0, size);
		cuda_memset(GPU, d_Tcg, 0, size);
		cuda_memset(GPU, d_Ucg, 0, size);
		cuda_memset(GPU, d_R0cg, 0, size);

		// initial value
		if (!Iseed) {
			cuda_memcpy(GPU, d_Ecg, d_Bcg, size, cudaMemcpyDeviceToDevice);
		}
		else {
			cuda_memset(GPU, d_Ecg, 0, size);

			// r0* : random vector
			//srand((unsigned int) time(NULL));
			srand(Iseed);
			for (int64_t n = 0; n < NE; n++) {
				Pcg[n].r = rand() / (double)RAND_MAX;
				Pcg[n].i = rand() / (double)RAND_MAX;
			}
			cuda_memcpy(GPU, d_Pcg, Pcg, size, cudaMemcpyHostToDevice);
		}

		// solve
		bicgstab(ifreq, NE, d_Acg, d_Ecg, d_Bcg, d_Rcg, d_Pcg, d_Qcg, d_Tcg, d_Ucg, d_R0cg, NIter, Residual[ifreq], Eaverage[ifreq], fp, iout);

		// copy device to host
		cuda_memcpy(GPU, Ecg, d_Ecg, size, cudaMemcpyDeviceToHost);

		// copy E to Ex+Ey+Ez
		copyE(ifreq);
	}

	// free
	free(Bcg);
	free(Ecg);
	memfree_cg_gpu();
}
