Good afternoon, from what ive been working on (Matlab and CUDA), here is a possible solution for anyone interested in Generic Matrix Multiplication. I use the cublas library, but the following code is to demonstrate how its possible to call cublas directly from cuda. I hope it serves all who are interested. Good luck to everyone in CUDA, David Lisin.
This is the matrixMult.m file:
[codebox]
dimension1=4000;
dimension2=4000;
A = diag(rand(dimension1,dimension1))*ones(1,dimension2);
B=ones(dimension2,dimension1);
disp(‘Matlab:’)
tic
C1=A*B;
toc
disp(‘Cuda:’)
tic
C2=matrixMulCUDA(A,B);
toc
[/codebox]
This is the actual cuda file (matrixMulCUDA.cu):
[codebox]
#pragma comment (lib,“C:\CUDA\lib\cublas.lib”)
#include “cuda.h”
#include “mex.h”
#include <string.h>
#include “cublas.h”
/* Gateway function */
void mexFunction(int nlhs, mxArray *plhs,int nrhs, const mxArray *prhs){
int mMatrizA;
int nMatrizA;
int mMatrizB;
int nMatrizB;
int mSalidaFuncion;
int nSalidaFuncion;
const mxArray *paramMatrizA;
double *matrizA;
const mxArray *paramMatrizB;
double *matrizB;
double *matrizAf, *matrizA_gpu;
double *matrizBf, *matrizB_gpu;
double *salidaFuncion;
double *salidaFuncionf, *salidaFuncion_gpu;
int j;
/* ############################################################
######## */
/* Get the dimensions of the matrizA */
paramMatrizA=prhs[0];
/* We collect the number of params */
mMatrizA= mxGetM(paramMatrizA);
nMatrizA=mxGetN(paramMatrizA);
/* ############################################################
######## */
/* Get the dimensions of the matrizB */
paramMatrizB=prhs[1];
/* We collect the number of params */
mMatrizB= mxGetM(paramMatrizB);
nMatrizB= mxGetN(paramMatrizB);
/* ############################################################
######## */
/*###########################################################
############*/
mSalidaFuncion=mMatrizA;
nSalidaFuncion=nMatrizB;
/*###########################################################
############*/
/* ############################################################
######## */
/* Create an mxArray for the output data */
if(!(plhs[0]=mxCreateDoubleMatrix(mSalidaFuncion,nSalidaFunc
ion,mxREAL))){
printf("\nERROR\n");
return;
}
/* ############################################################
######## */
/* Creamos los arrays de input y arrays de salida del GPU */
cudaMalloc((void**)&matrizA_gpu,sizeof(double)*mMatrizA*nMatrizA);
cudaMalloc((void**)&matrizB_gpu,sizeof(double)*mMatrizB*nMatrizB);
cudaMalloc((void**)&salidaFuncion_gpu,sizeof(double)*mSalidaFuncion*nSalidaF
uncion);
/* ############################################################
######## */
/* Recogemos los datos del input */
matrizA=mxGetPr(prhs[0]);
matrizB=mxGetPr(prhs[1]);
/* ############################################################
######## */
matrizAf=(double*)mxMalloc(sizeof(double)*mMatrizA*nMatrizA)
;
for(j=0;j<nMatrizA*mMatrizA;j++){
matrizAf[j]=matrizA[j];
}
cudaMemcpy(matrizA_gpu,matrizAf,sizeof(double)*mMatrizA*nMat
rizA,cudaMemcpyHostToDevice);
/* ############################################################
######## */
matrizBf=(double*)mxMalloc(sizeof(double)*mMatrizB*nMatrizB)
;
for(j=0;j<nMatrizB*mMatrizB;j++){
matrizBf[j]=matrizB[j];
}
cudaMemcpy(matrizB_gpu,matrizBf,sizeof(double)*mMatrizB*nMat
rizB,cudaMemcpyHostToDevice);
/* ############################################################
######## */
/* ############################################################
####### */
/* Reservamos la memoria de la salida de doble precision */
salidaFuncionf=(double*)mxMalloc(sizeof(double)*mSalidaFunci
on*nSalidaFuncion);
/* ############################################################
*/
/* READY TO CALL SGEMM */
//(void) cublasDgemm ('n','n',mMatrizA,nMatrizB,nMatrizA,1,matrizA,mMatrizA,matrizB_g
pu,mMatrizB,0,salidaFuncion_gpu,mMatrizA);
/* READY TO CALL SGEMM */
(void) cublasDgemm ('n','n',mMatrizA,nMatrizB,nMatrizA,1,matrizA_gpu,mMatrizA,matri
zB_gpu,mMatrizB,0,salidaFuncion_gpu,mMatrizA);
// block until the device has completed
cudaThreadSynchronize();
/* Copy result back to host */
cudaMemcpy( salidaFuncionf, salidaFuncion_gpu, sizeof(double)*mSalidaFuncion*nSalidaFuncion, cudaMemcpyDeviceToHost);
/* Create a pointer to the output data */
salidaFuncion = mxGetPr(plhs[0]);
for(j=0;j<nSalidaFuncion*mSalidaFuncion;j++){
salidaFuncion[j]=(double)salidaFuncionf[j];
}
/* Clean-up memory on device and host */
mxFree(matrizAf);
mxFree(matrizBf);
mxFree(salidaFuncionf);
cudaFree(matrizA_gpu);
cudaFree(matrizB_gpu);
cudaFree(salidaFuncion_gpu);
}
[/codebox]
The output is:
Matlab:
Elapsed time is 21.544856 seconds.
Cuda:
Elapsed time is 2.545485 seconds.
So not too bad :)
Im using double precision matrix multiplication DGEMM, but if you dont need double precision, use single precision SGEMM, as it has better times.
Again, good luck to everyone working with CUDA :)