不同工具下的矩阵乘法速度测试
来源:互联网 发布: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常用的那几招(引用于葵花宝典)》。
最后欢迎大家测试,如果有不对的地方欢迎大家指正。
- 不同工具下的矩阵乘法速度测试
- 测试矩阵乘法的例子
- 矩阵乘法测试
- Hadoop的安装、测试、以及为伪分布式下矩阵乘法的实现
- c++ struct下的矩阵乘法
- 从矩阵乘法的不同计算方式来看局部性原理
- 使用iperf工具测试infiniband的速度
- cuda开发矩阵乘法测试你的GPU效率
- CUDA开发矩阵乘法测试你的GPU效率
- 【矩阵乘法】:矩阵乘法的基本实现
- 测试矩阵连续运算的速度问题的代码
- 矩阵的思考-矩阵乘法
- 使用不同方法拷贝字节流文件的速度测试
- 通过矩阵乘法看内存访问对CPU运算速度的影响
- 矩阵的乘法
- 矩阵的乘法问题
- 矩阵乘法的模板
- 矩阵乘法的疑惑
- js oop写法小例子
- ECharts使用心得(最初的版本)
- CSS(3) Text(文本)
- 第三方登录与分享
- Web高并发量网站解决方案
- 不同工具下的矩阵乘法速度测试
- 循环队列和链队列
- 《剑指offer》:[40]数组中只出现一次的数字
- 【c++程序】数字颠倒
- JAVA_ListIterator
- struts2-拦截器(一)
- Android视图SurfaceView的实现原理分析
- PageRank算法
- Linux大小端判断