/*
* Ce fichier MEX + (CUDA + GPUmat) doit être compilé en deux étapes:
*
*	nvcc -c AddMatGPUmat.cu -Xcompiler -fPIC -arch sm_13 -I /usr/local/math/matlab/matlab-2010a/extern/include/
*	/usr/local/math/matlab/matlab-2010a/bin/mex AddMatGPUmat.o -L /usr/local/cuda/cuda31/lib/ -lcudart -cxx
*
* Le fichier généré, "AddMatGPUmat.mexglx", est le programme à exécuter dans MATLAB, une fois le toolbox "GPUmat" démarré. 
* Pour activer le toolbox, entrez "run /usr/local/math/GPUmat/GPUmat.260/GPUstart" depuis l'invite de MATLAB.
*
*
* Le programme retourne la somme des deux vecteurs (ou matrices) de type GPUdouble qui lui sont fournis à l'entrée; le calcul 
* est effectué par le GPU. 
* La syntaxe pour appeler la fonction est la suivante:
* 	AddMatGPUmat(getPtr(A), getPtr(B), getPtr(C), N), 
* où A et B sont les matrices de formats m x n et de type GPUdouble à additionner, C est la matrice réponse et où N = m x n.
*
* Pour importer des objets de précision simple, remplacez les "double" par "float". Pour les autres types de données (entiers,
* signés ou non, etc.), remplacez "double" par le type désiré.
*/
#include <stdint.h>	/* Cette librairie doit être incluse pour pouvoir utiliser le type "uintptr_t" employé dans le code. */
#include "mex.h"
#include "cuda.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 moins de "threads" que les autres 
* blocs.
*/
__global__ void Addition(double *A, double *B, double *C, unsigned int N)
{
	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[])
{
	if(nrhs != 4)
		mexErrMsgTxt("\nSyntaxe: AddMatGPUmat(getPtr(A), getPtr(B), getPtr(C), N), où A et B sont les matrices de formats m x n et de type GPUdouble à additionner, C est la matrice réponse et où N = m x n.\n");
/*
* Lecture des vecteurs (ou matrices) de type "GPUdouble" et autres données fournies dans MATLAB. Ici, les variables "Agpu", 
* "Bgpu", et "Cgpu" se trouvent déjà dans la mémoire du GPU. 
*/		
	double *Agpu = (double *) ( (uintptr_t) mxGetScalar(prhs[0]) );
	double *Bgpu = (double *) ( (uintptr_t) mxGetScalar(prhs[1]) );
	double *Cgpu = (double *) ( (uintptr_t) mxGetScalar(prhs[2]) );
	unsigned int N = (unsigned int) mxGetScalar(prhs[3]);
/*	
* 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 = (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>>>(Agpu, Bgpu, Cgpu, N);	
}