Matrix Multiplication with CUBLAS and MATLAB

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 :)