/*
* Ce fichier MEX + CUDA doit être compilé en deux étapes:
*
*	nvcc -c AddMatCUDA.cu -Xcompiler -fPIC -arch sm_13 -I /usr/local/math/matlab/matlab-2010a/extern/include/
*	/usr/local/math/matlab/matlab-2010a/bin/mex AddMatCUDA.o -L /usr/local/cuda/cuda31/lib/ -lcudart -cxx
*
* Le fichier généré, "AddMatCUDA.mexglx", est le programme à exécuter dans MATLAB.
*
*
* Le programme retourne la somme des deux vecteurs (ou matrices) de précision double qui lui sont fournis à l'entrée; le 
* calcul est effectué par le GPU.
*
* Pour importer des objets de précision simple, remplacez les "double" par "float" ainsi que "mxDOUBLE_CLASS" (dans la 
* fonction "mxCreateNumericMatrix") par "mxSINGLE_CLASS". Pour les autres types de données (entiers, signés ou non, etc.), 
* remplacez "double" par le type désiré et "mxDOUBLE_CLASS" par le type correspondant (voir
* http://www.mathworks.com/help/techdoc/apiref/mxcreatenumericmatrix.html pour la liste).
*/
#include "cuda.h"
#include "mex.h"
/*
* La fonction qui suit additionne les vecteurs "A" et "B" et place le résultat dans "C". Cette fonction, de type 
* "__global__" est démarrée par le CPU dans le code principal mais s'exécute sur le GPU. Ce dernier démarrera autant de copies 
* de cette fonction qu'il y aura de "threads" à exécuter. Les éléments à additionner sont indicés par "i" et celui-ci est 
* obtenu en connaissant la position du "thread" dans le bloc (threadIdx), la position du bloc dans la grille (blockIdx) ainsi 
* qu'en connaissant le format de la grille (gridDim) et des blocs (blockDim).
* La condition "if" doit être ajoutée pour être certain de ne pas lire et écrire en dehors de la mémoire prévue pour les 
* vecteurs "A", "B" et "C". En effet, le dernier bloc de la grille pourrait contenir un nombre inférieur de "threads".
*/
__global__ void Addition(unsigned int N, double *A, double *B, double *C)
{
	unsigned int i = (blockIdx.y*gridDim.x + blockIdx.x)*blockDim.x + threadIdx.x; 
		
	if (i < N)
		C[i] = A[i] + B[i];
}
/*
* Ci-dessous, la fonction "mexFunction" fait le lien entre le code C et MATLAB; c'est l'équivalent de la fonction "main", pour 
* un programme s'exécutant sous MATLAB.
*
* Description courte des arguments de "mexFunction":
*	nlhs: Nombre d'arguments de sortie fournis à la fonction "AddMat".
*	*plhs[]: Tableau de pointeurs vers les vecteurs (ou matrices) de sortie passés à la fonction "AddMat".
*	nrhs: Nombre d'arguments d'entrée fournis à la fonction "AddMat".
*	*prhs[]: Tableau de pointeurs vers les vecteurs (ou matrices) d'entrée passés à la fonction "AddMat".
*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	unsigned int m, n;
	double *dataIN1cpu, *dataIN2cpu, *dataOUTcpu;
	double *dataIN1gpu, *dataIN2gpu, *dataOUTgpu;
	
	if(nrhs != 2)
		mexErrMsgTxt("Deux arguments (vecteurs ou matrices) doivent être fournis à la fonction.\n");
	
	m = mxGetM(prhs[0]);			/* Nombre de rangées de la première matrice d'entrée.  			*/
	n = mxGetN(prhs[0]);			/* Nombre de colonnes de la première matrice d'entrée. 			*/
		
	if( (m != mxGetM(prhs[1])) || (n != mxGetN(prhs[1])) )
		mexErrMsgTxt("Les deux vecteurs (ou matrices) doivent être de même format.\n");
	
	dataIN1cpu = (double *) mxGetPr(prhs[0]);		/* Les données en entrée sont des tableaux unidimensionnels, peu importe */
	dataIN2cpu = (double *) mxGetPr(prhs[1]);		/* qu'il s'agissent de matrices ou de vecteurs.				 */
/*
* Création des tableaux d'entrée sur la mémoire GPU.
*/
	cudaMalloc( (void **) &dataIN1gpu, sizeof(double)*m*n );
	cudaMalloc( (void **) &dataIN2gpu, sizeof(double)*m*n );
	cudaMalloc( (void **) &dataOUTgpu, sizeof(double)*m*n );
/*
* Transfert des données d'entrée sur le CPU vers les tableaux d'entrée tout juste créés sur le GPU. Le dernier paramètre
* accepté par "cudaMemcpy", réglé à "cudaMemcpyHostToDevice" ci-dessous, informe la fonction de l'orientation désirée pour le 
* transfert des données. Dans le cas présent, on souhaite que le transfert s'effectue depuis la mémoire "host" (CPU) vers la 
* mémoire "device" (GPU).
*/	
	cudaMemcpy(dataIN1gpu, dataIN1cpu, sizeof(double)*m*n, cudaMemcpyHostToDevice);
	cudaMemcpy(dataIN2gpu, dataIN2cpu, sizeof(double)*m*n, cudaMemcpyHostToDevice);	
/*	
* Avant d'appeler la fonction "__global__", il faut choisir le format de chaque bloc de "threads" ainsi que le nombre de 
* "threads" les composant. Ensuite, il faut déterminer le format de la grille ainsi que le nombre de blocs qui la 
* constitueront. Ici, les blocs ont été choisi en format unidimensionnel puisque les données d'entrée (MATLAB) sont des 
* vecteurs (MATLAB indexe tous ses tableaux par un seul indice). Le format de la grille est déterminé comme suit: si le nombre 
* de blocs nécessaires est supérieur à 65535, alors la grille est choisie bidimensionnelle. Sinon, le programme sélectionne 
* une grille unidimensionnelle. Ceci est nécessaire puisque CUDA ne permet pas qu'une des dimensions de la grille 
* (gridDim.x, gridDim.y et gridDim.z) dépasse 65535. 
*/	
	dim3 DimenBlocs(256, 1, 1);						/* Blocs unidimensionnels de 256 threads.*/
	dim3 DimenGrille;
	unsigned int NbBlocsNec = (m*n + DimenBlocs.x - 1)/DimenBlocs.x;
	
	if(NbBlocsNec > 65535)							/* Grille bidimensionnelle nécessaire.   */
	{
		DimenGrille.x = 65535; 
		DimenGrille.y = (NbBlocsNec + 65535 - 1)/ 65535;
		DimenGrille.z = 1;
	} 
	else									/* Grille unidimensionnelle suffisante.  */
	{
		DimenGrille.x = NbBlocsNec;
		DimenGrille.y = 1;
		DimenGrille.z = 1;
	}		
/*
* Finalement, on appelle la fonction complétant le calcul sur le GPU avec la configuration bloc/grille choisie.
*/
	Addition<<<DimenGrille, DimenBlocs>>>(m*n, dataIN1gpu, dataIN2gpu, dataOUTgpu);
/*
* On transfère à présent le résultat de l'addition vers le tableau de sortie sur le CPU, après avoir créé ce tableau. Le 
* transfert des données s'effectue maintenant du "device" (GPU) vers le "host" (CPU), donc on emploiera 
* "cudaMemcpyDeviceToHost" comme dernier paramètre de "cudaMemcpy".
*/
	plhs[0] = mxCreateNumericMatrix(m, n, mxDOUBLE_CLASS, mxREAL);		
	dataOUTcpu = (double *) mxGetPr(plhs[0]);
	cudaMemcpy(dataOUTcpu, dataOUTgpu, sizeof(double)*m*n, cudaMemcpyDeviceToHost);
/*
* Avant de quitter, on libère la mémoire occupée par les variables allouées avec "cudaMalloc".
*/	
	cudaFree(dataIN1gpu), cudaFree(dataIN2gpu), cudaFree(dataOUTgpu);
}