#include <stdlib.h>
#include <stdio.h>
#include <math.h>

///Fonction qui fait la réduction
/**	@param in : tableau à réduire, sa taille est égale au double du nombre de threads dans la grille
 * 	@param sizeIn : taille du tableau in
 * 	@param out : tableau de valeurs réduites, sa taille est égale au nombre de block
 * 	@param sizeOut : taille du tableau out
*/
__global__ void kernelReduction(float *in, size_t sizeIn, float *out, size_t sizeOut){
	//On sait que la taille de in est 2*blockDim.x
	int i = 2*(blockIdx.x*blockDim.x + threadIdx.x);
	int indexCalcul = 1;
	while(indexCalcul < 2*blockDim.x){
		if(i + indexCalcul < sizeIn){           //on test si le prochain élément est toujours dans le tableau in
			in[i] += in[i + indexCalcul];
		}
		indexCalcul *= 2;
		__syncthreads();                       //on synchronise tout les threads d'un même block, avant de passé à la suite
	}
	if(threadIdx.x == 0) out[blockIdx.x] = in[i];
}

///fonction qui permet d'initialiser les tableaux de valeurs sur l'hote
/**	@param h_in : tableaux dont on veut initialiser les valeurs
 * 	@param ntot : nombre d'éléments du tableau
*/
__host__ void init_vec(float *h_in, int ntot){
	for(int i = 0 ; i < ntot ; i++){
		h_in[i] = sinf(float(i));
		//h_in[i] = 1.0;   //pour faire un test simple
	}
}

///fonction qui permet d'initialiser les tableaux de valeurs sur l'hote
/**	@param h_in : tableaux dont l'on veut afficher les valeurs
 * 	@param ntot : nombre d'éléments du tableau
*/
__host__ void printf_vec(float *h_in, int ntot){
	for(int i = 0 ; i < ntot ; i++){
		printf("%f,\t", h_in[i]);
		if(i%8 == 0 && i != 0) printf("\n");  //on affiche 8 colonnes
	}
	printf("\n");
}

///fonction qui permet de vérifier les résultats contenus dans les tableaux de valeurs
/**	@param sun : somme que l'on a trouvé avec le GPU
 * 	@param h_in : tableaux de valeurs
 * 	@param ntot : nombre d'éléments du tableau
*/
__host__ void verifCalcul(float sumGPU, float sumCPU){
	float err = fabsf((sumGPU - sumCPU)/sumGPU);
	printf("GPU sum : %.4e\n", sumGPU);
	printf("CPU sum : %.4e\n", sumCPU);
	if(err < 1.e-4){
		printf("TEST PASSED (err %.4e < 1.e-4).\n", err);
	}else{
		printf("TEST FAILED (err %.4e > 1.e-4).\n", err);
	}
}

float reductionCPU(float *h_in, int ntot){
	float sum = 0.0;
	for(int i = 0 ; i < ntot ; ++i){
		sum += h_in[i];
	}
	return sum;
}

void checkNoErr(int linenr, float * d_sum, float * d_bl, float * d_in, float * h_in){
	cudaError_t err = cudaGetLastError();
	if (err != cudaSuccess) {
		printf("Erreur ligne %d : %s \n", linenr, cudaGetErrorString(err));
		cudaFree(d_sum);
		cudaFree(d_bl);
		cudaFree(d_in);
		cudaFree(h_in);
		abort();
	}
}

int main(int argc, char **argv){
	float sum;
	size_t nthreads, nblocks = 16, ntot = 1024;         //le nombre de valeurs dans un tableau doit être une puissance de 2
	
	nthreads = ntot/(2*nblocks);                   //le nombre de threads doit être ntot/nblocks avec ntot et nblocks puissances de 2
	
	printf("Ntot     : %lu\n", ntot);
	printf("nthreads : %lu\n", nthreads);
	printf("nblocks  : %lu\n", nblocks);
	
	float *d_sum, *d_bl, *h_bl, *d_in, *h_in;
	
	cudaMallocHost((void**)&h_in, ntot*sizeof(float));
	cudaMallocHost((void**)&h_bl, nblocks*sizeof(float));
	if(h_in == NULL){printf("Impossible d'allouer h_in == NULL\n");return 1;}
	
	cudaMalloc((void**)&d_sum, sizeof(float));
	cudaMalloc((void**)&d_bl, nblocks*sizeof(float));
	cudaMalloc((void**)&d_in, ntot*sizeof(float));
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
	init_vec(h_in, ntot);
	float sommeCPU = reductionCPU(h_in, ntot);
// 	printf_vec(h_in, ntot);
	cudaMemcpy(d_in, h_in, ntot*sizeof(float), cudaMemcpyHostToDevice);
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
	dim3 dimGrid(nblocks, 1, 1);
	dim3 dimBlock(nthreads, 1, 1);
	
	dim3 dimGridHalf(nblocks/2, 1, 1);
	dim3 unique(1,1,1);
	
	//réduction à l'intérieur de chaque block
	kernelReduction<<<dimGrid, dimBlock>>>(d_in, ntot, d_bl, nblocks);             //on récupère le résultat dans d_bl
	//rédution entre les blocks, un seul block fait la réduction
	kernelReduction<<<unique, dimGridHalf>>>(d_bl, nblocks, d_sum, 1);             //on récupère le résultat dans d_sum
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
	cudaThreadSynchronize();
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
	cudaMemcpy(h_in, d_in, ntot*sizeof(float), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_bl, d_bl, nblocks*sizeof(float), cudaMemcpyDeviceToHost);
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
	cudaMemcpy(&sum, d_sum, sizeof(float), cudaMemcpyDeviceToHost);
	checkNoErr(__LINE__, d_sum, d_bl, d_in, h_in);
// 	printf("\nAprès le calcul h_in :\n");
// 	printf_vec(h_in, ntot);
// 	printf("\nAprès le calcul h_dl :\n");
// 	printf_vec(h_bl, nblocks);
	verifCalcul(sum, sommeCPU);

	cudaFree(d_sum);
	cudaFree(d_bl);
	cudaFreeHost(h_bl);
	cudaFree(d_in);
	cudaFreeHost(h_in);
	return 0;
}