不同工具下的矩阵乘法速度测试

来源:互联网 发布:java两种多态机制 编辑:程序博客网 时间:2024/03/29 19:53

      今天花了一些时间将基于以上几种工具(cuda&arrayfire&matlab)的矩阵乘法的速度进行了测试比较,验证了一些想法吧。

首先是c(CPU)的乘法测试:写的有点繁琐,后面在cuda程序中进行了综合

C.matrix.cpp

void Matrix_print(double **a,long nl, long nh)/*矩阵的输出*/{int i,j;for(i=1;i<=nl;i++){for(j=1;j<=nh;j++){cout<<a[i][j]<<",";       }cout<<endl;}}double **Matrix(long nrl,long nrh,long ncl,long nch)/*分配一个矩阵 m[nrl..nrh][ncl..nch]*/{long i,nrow=nrh-nrl+1,ncol=nch-ncl+1;double **m;m=(double **)malloc((size_t)((nrow)*sizeof(double *))); //分配row个行指计if(!m) nrerror("allocation matrix row failure!");m-=nrl;m[nrl]=(double *)malloc((size_t)((nrow*ncol)*sizeof(double))); //分配row*col个内存空间    if(!m[nrl]) nrerror("allocation matrix col failure!");m[nrl]-=ncl;for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;return m;}void CMatrix_multiply(complex **a,int row1,int col1,complex **b,int row2,int col2,complex **c)/*复数矩阵乘法 c=a x b*/{<span style="white-space:pre"></span>int i,j,k;<span style="white-space:pre"></span>if(col1!=row2) nrerror("第一个矩阵的列数不等于第二个矩阵的行数,不能相乘!");<span style="white-space:pre"></span>for(i=1;i<=row1;i++)<span style="white-space:pre"></span>for(j=1;j<=col2;j++)<span style="white-space:pre"></span>{             c[i][j].real=0;<span style="white-space:pre"></span> c[i][j].imag=0;             for(k=1; k<=col1;k++) <span style="white-space:pre"></span> {<span style="white-space:pre"></span> c[i][j]=Cadd(c[i][j],Cmultiply(a[i][k],b[k][j]));<span style="white-space:pre"></span> }<span style="white-space:pre"></span>}}
</pre><pre code_snippet_id="1721460" snippet_file_name="blog_20160619_3_9040794" name="code" class="cpp">#include "stdio.h"#include "stdlib.h"#include<vector>#include<time.h>void main(){int const NN=8;double **data1,**data2,**result;data1=Matrix(1,NN,1,NN);data2=Matrix(1,NN,1,NN);result=Matrix(1,NN,1,NN);for(int i=1;i<NN+1;i++){for(int j=1;j<NN+1;j++){data1[i][j]=rand()/(double)(RAND_MAX);data2[i][j]=rand()/(double)(RAND_MAX);}}//Matrix_print(data1,NN,NN);//Matrix_print(data2,NN,NN);clock_t start,end;      start = clock(); Matrix_multiply(data1,NN,NN,data2,NN,NN,result);//Matrix_print(result,NN,NN);//system("pause");end = clock();      printf("time=%f\n",((double)end-start)/CLK_TCK); free_Matrix(data1,1,NN,1,NN);free_Matrix(data2,1,NN,1,NN);free_Matrix(result,1,NN,1,NN);system("pause");}
CUDA测试代码

matrixCuda.cu

#include <cuda_runtime.h> #include<cuda.h>#include "matrixmul.h"#include <device_launch_parameters.h>#define BLOCK_SIZE 8using namespace std;// 核函数__global__ void MatrixMulKernel( double* Md, double* Nd,double* Pd,int WidthM, int WidthN){//计算Pd中的元素索引int Row = blockIdx.y * BLOCK_SIZE + threadIdx.y; //行int Col = blockIdx.x * BLOCK_SIZE + threadIdx.x; //列float Pvalue = 0.0;for (int k = 0; k < WidthM; k++){Pvalue += Md[Row * WidthM + k] * Nd[k * WidthN + Col];}//每个线程负责计算P中的一个元素Pd[Row * WidthN + Col] = Pvalue;}// 矩阵乘法void matrixMul( double* M,double* N,double* P,int WidthM, int WidthN, int HeightM, int HeightN){double *Md, *Nd, *Pd;    //字节长度int sizeM = WidthM * HeightM * sizeof(double);    int sizeN = WidthN * HeightN * sizeof(double);    int sizeP = WidthN * HeightM * sizeof(double);    cudaMalloc((void**)&Md, sizeM);cudaMalloc((void**)&Nd, sizeN);cudaMalloc((void**)&Pd, sizeP);//Copies a matrix from the memory* area pointed to by src to the memory area pointed to by dstcudaMemcpy(Md, M, sizeM, cudaMemcpyHostToDevice);cudaMemcpy(Nd, N, sizeN, cudaMemcpyHostToDevice);        dim3 dimGrid( HeightM / BLOCK_SIZE , WidthN / BLOCK_SIZE);dim3 dimBlock( BLOCK_SIZE,  BLOCK_SIZE);MatrixMulKernel<<< dimGrid, dimBlock >>>(Md, Nd, Pd, WidthM, WidthN);<pre name="code" class="cpp">
cudaMemcpy(P, Pd, sizeP, cudaMemcpyDeviceToHost);//释放设备上的矩阵cudaFree(Md);cudaFree(Nd);cudaFree(Pd); cudaDeviceReset();}

matrixmul.h

#ifndef __MATRIXMUL_H__#define __MATRIXMUL_H__//extern void matrixMul(double* M, double* N, double* P,int Width);extern void matrixMul(double* M,double* N,double* P,int WidthM, int WidthN, int HeightM, int HeightN);#endif


cuda_array_ctest.cpp

#include "stdio.h"#include "stdlib.h"#include <iostream>#include "time.h"#include "matrixmul.h"#include "arrayfire.h"#include <cstdio>#include <cstdlib>using namespace std;using namespace af;/*产生随机矩阵M行N列,矩阵元素0~1*/void matrixRan (double *DATA,int M,int N){for(int i=0;i<M;i++)for(int j=0;j<N;j++){            DATA[i*N+j]=rand()/(double)(RAND_MAX);}}/*实数矩阵乘法*/void Matrix_multiply(double *DATAa,int M,int N,double *DATAb,int H,int K,double *Result){if(N!=H) printf("第一个矩阵的列数不等于第二个矩阵的行数,不能相乘!");for(int i=0;i<M;i++)for(int j=0;j<K;j++){             double sum=0;             for(int k=0; k<N;k++)   { sum+=DATAa[i*M+k]*DATAb[k*H+j];         // cout<<c[i][j];     } Result[i*M+j]=sum; }}/*打印矩阵*/void matrixPrint (double *DATA,int M,int N){for(int i=0;i<M;i++){for(int j=0;j<N;j++){                   cout<<DATA[i*N+j]<<",";        }cout<<endl;   }}void main(){double *DATA1,*DATA2,*ResultGPU,*ResultArr,*ResultCPU;int const MM=2048;int const NN=2048;int const HH=2048;//DATA1的列要与DATA2的行相等int const KK=2048;//矩阵维度,可以不同维度进行测试;DATA1 = (double*)malloc(sizeof(double) * MM * NN);      DATA2 = (double*)malloc(sizeof(double) * HH * KK);      ResultGPU= (double*)malloc(sizeof(double) * MM * KK);//GPU的计算结果//array ResultArr;//arrayfire的计算结果ResultCPU= (double*)malloc(sizeof(double) * MM * KK);//CPU的计算结果/*产生矩阵DATA1和DATA2*/matrixRan (DATA1,MM,NN);matrixRan (DATA2,MM,NN);/*GPU计算时间*/clock_t start,end;         start = clock(); matrixMul(DATA1,DATA2,ResultGPU,MM,HH,NN,KK);end = clock();          printf("GPUtime=%f\n",((double)end-start)/CLK_TCK); //free(DATA1);//free(DATA2);free(ResultGPU);//matrixPrint(DATA1,MM,NN);//matrixPrint(DATA2,HH,KK);//matrixPrint(ResultGPU,MM,KK);/*Arrayfire计算时间(这段一直编译不起来重新弄一下了)    start = clock(); array DATA11(MM,NN,DATA1);//将数组转化成arrayarray DATA22(MM,NN,DATA2);af_print(DATA11);af_print(DATA22);ResultArr=matmul(DATA11,DATA22);end = clock();    printf("time=%f\n",((double)end-start)/CLK_TCK); *//*CPU计算时间*/           start = clock(); Matrix_multiply(DATA1,MM,NN,DATA2,HH,KK,ResultCPU);//matrixPrint(ResultCPU,MM,KK);end = clock();     printf("CPUtime=%f\n",((double)end-start)/CLK_TCK);system("pause");}


Arrayfire模式:在cuda中试图写了,可是一起编译不起来,有兴趣可以尝试一下子。

#include <arrayfire.h>#include <cstdio>#include <cstdlib>#include<vector>using namespace af;std:: vector<float> vecl;int main(int argc, char *argv[]){    try {        // Select a device and display arrayfire info        int device = argc > 1 ? atoi(argv[1]) : 0;        af::setDevice(device);        af::info();timer::start();int const NN=1024;        array DATA1 = randu(NN, NN);double k=0.0;k=det<double>(DATA1);printf("%f",k);af_print(DATA1);array DATA2 =randu(NN,NN);af_print(DATA2);array ResultArr=matmul(DATA1,DATA2);af_print(ResultArr);printf("elapsed seconds: %g\n", timer::stop());    } catch (af::exception& e) {        fprintf(stderr, "%s\n", e.what());        throw;    }system("pause");    return 0;}

C-mex模式(matlab&GPU)

主要包括matrixCuda.cu,matrixMul.h,matrixMulCuda.cpp,matrixMul.m,其中matrixCuda.cu,matrixMul.上面已经写出来了,直接copy即可,另外还有两个文件如下所示:

matrixMulCuda.cpp

#include "mex.h"#include "matrixMul.h"void mexFunction(int nlhs, mxArray *plhs[], int nrhs, mxArray *prhs[]){     if (nrhs != 2)             mexErrMsgTxt("Invaid number of input arguments");                if (nlhs != 1)                     mexErrMsgTxt("Invalid number of outputs");                          if (!mxIsDouble(prhs[0]) && !mxIsDouble(prhs[1]))                                    mexErrMsgTxt("Input matrix data type must be Double");  /****************************************int numRowsM = (int)mxGetM(prhs[0]);    int numColsM = (int)mxGetN(prhs[0]);    int numRowsN = (int)mxGetM(prhs[1]);    int numColsN = (int)mxGetN(prhs[1]);   ******************************************/        int  HeightM = (int)mxGetM(prhs[0]);            int   WidthM = (int)mxGetN(prhs[0]);           int  HeightN = (int)mxGetM(prhs[1]);            int   WidthN = (int)mxGetN(prhs[1]);         if ((WidthM!=HeightN)||(WidthN!=HeightM))         mexErrMsgTxt("Invalid size. The sizes of two matrix must be limit");       double *M=(double*)mxGetData(prhs[0]);       double *N=(double*)mxGetData(prhs[1]);            plhs[0] = mxCreateDoubleMatrix(HeightM, WidthN, mxREAL);       double *P=(double*)mxGetData(plhs[0]);       matrixMul(M,N,P,WidthM,WidthN,HeightM,HeightN);        }

matrixMul.m

% RunmatrixMulclearclcdisp('1. nvcc MulMat.cu compiling ...');%system('nvcc -c matrixMul.cu -ccbin "C:\Program Files\Microsoft Visual Studio 10.0\VC\bin"'); system('nvcc -c matrixMul.cu  -gencode arch=compute_50,code=sm_50 -ccbin "D:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\bin"')mex matrixMulCuda.cpp matrixMul.obj -lcudart -L"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v7.5\lib\x64";disp('2.Input two matrix:') A=rand(2048); B=rand(2048);%A=ones(16,16);%B=ones(16,16);disp('3.Compare with the cpu and gpu result:')disp('gpu Result:')ticC_gpu =matrixMulCuda(B,A);tg=toc; disp(['GPU并行计算时间:',num2str(tg),'秒'])  disp('cpu Result:')ticC_cpu = A*B;tc=toc; disp(['CPU串行计算时间:',num2str(toc),'秒'])

至此4种方式的代码全部在上面了,可以帮助刚开始的童鞋熟悉各个不同的GPU编程方式,接下来就是测试结果了。

测试时间如下图表所示:




从上图和上表我们不难看出以下几点

1、随着数据量的增加,GPU加速效果越来越明显,arryfire趋于稳定;

2:随着数据量的增加,除了arrayfire各个工具处理的时间的都会有不同程度的增加,尤其是CPU的最为明显,所以在不同的方式下GPU都会有不同程度的加速;

3:这里的cuda程序没有进行内存优化,可能用上共享内存的方式或许会更快,改变blocksize的大小也会有影响,这里偷懒将blocksize设置为8,没有测其他值了

4:根据《CBF中for循环变矩阵乘法的思想(arrayfire)》这篇博客测试来看,如果不解决编程中for循环的问题,运行时间会维持在一个相对较高的状态,所以关键问题还是在不同的工具下都要用并行的思想去写程序。

5:顺便推广一下arrayfire,运行速度稳定,包含了大量的矩阵运算工具,可参见博客《Arrayfire常用的那几招(引用于葵花宝典)》。


最后欢迎大家测试,如果有不对的地方欢迎大家指正。


0 0