《矩阵与向量乘法的cuda优化.ppt》由会员分享，可在线阅读，更多相关《矩阵与向量乘法的cuda优化.ppt（28页珍藏版）》请在三一办公上搜索。
1、矩阵与向量乘法的CUDA优化,风辰2010年12月11日2011年1月8日修订,1,目的,对于CUDA程序开发来说，优化往往是整个开发过程的核心，不同算法，不同存储器组织的程序性能往往差几十倍，本文通过一个简单的例子来展示CUDA开发中一些重要的因素对性能的影响。,2,假设读者拥有以下知识,拥有C语言编程的经验,最好拥有并行编程经验懂得CUDA,最好用CUDA写过代码,3,测试环境,Intel xeon 5405 2.0 GHzGeforce GTX 295(只使用单核)Gcc 4.3.3 CUDA toolkit 3.1只测试计算时间,不包括数据传输,4,符号说明,matrix:矩阵数据指针
2、,以行为主序或者列为主序存储v|vec:向量指针r:矩阵和向量乘的结果指针rowSize:表示矩阵的行数,也是r的长度columnSize:表示矩阵的列数,也是v的长度所有指向显存的指针加前缀d_,5,编译配置,矩阵尺寸8192*8192单精度编译选项-O3 funroll-loops msseCPU计时函数采用gettimeofday,clock,GPU计时函数采用CUDA event,6,串行C版本,算法:遍历矩阵行,每行和向量相乘,最终结果为一向量void mxv(const int rowSize,const int columnSize,const float*matrix,cons
3、t float*v,float*r)for(int i=0;i rowSize;i+)float re=0.0f;for(int j=0;j columnSize;j+)re+=matrixi*columnSize+j*vj;ri=re;,7,运行时间120 ms,不使用-O3运行耗时490 ms,简单SSE版本,算法:利用sse指令计算矩阵每行和向量的乘积void mxvSSE(const int rowSize,const int columnSize,const float*matrix,const float*v,float*r)_m128*mv=(_m128*)v;_m128*mm=
4、(_m128*)matrix;for(int i=0;i rowSize;i+)_m128 re=_mm_set_ps(0.0f,0.0f,0.0f,0.0f);for(int j=0;j columnSize/4;j+)re=_mm_add_ps(re,_mm_mul_ps(mmi*columnSize/4+j,mvj);float _attribute(aligned(16)a4;_mm_store_ps(a,re);ri=a0+a1+a2+a3;,运行时间99ms,8,SSE+openmp,算法:使用二线程并行计算行循环void mxvSSEOpenmp(const int rowSiz
5、e,const int columnSize,float*matrix,float*vec,float*r)_m128*mv=(_m128*)v;_m128*mm=(_m128*)matrix;#pragma omp parallel for num_threads(2)for(int i=0;i rowSize;i+)_m128 re=_mm_set_ps(0.0f,0.0f,0.0f,0.0f);for(int j=0;j columnSize/4;j+)re=_mm_add_ps(re,_mm_mul_ps(mmi*columnSize/4+j,mvj);float _attribute
6、(aligned(16)a4;_mm_store_ps(a,re);ri=a0+a1+a2+a3;,运行时间50ms,9,CUDA优化注意事项,一、选择好的并行方式选择好的算法，以发掘更多的数据并行性二、保持SM忙碌尽量利用所有的SM参与计算，可以通过加大数据量或减小线程块大小达到目的三、优化存储器利用保证全局存储器合并访问使用速度更快的constant或shared存储器,10,CUDA-nave版本,算法:每个CUDA线程计算矩阵的一行与向量乘积static void _global_ mxvNaive(int rowSize,int columnSize,int columnPitch,
7、const float*d_matrix,const float*d_vec,float*d_r)uint id=blockDim.x*blockIdx.x+threadIdx.x;if(rowSize=id)return;float temp=0.0f;#pragma unroll 4 for(int i=0;i columnSize;i+)temp+=d_matrixid*columnPitch+i*d_veci;d_rid=temp;,耗时150 ms 串行120ms,11,CUDA-nave,为什么比串行还慢？columnPitch的作用是什么？访问d_matrix没有满足合并访问的要
8、求什么是合并访问？,12,合并访问,一句话：相邻线程访问段对齐的相邻地址为什么说访问d_matrix没有满足合并访问要求for(int i=0;i columnSize;i+)temp+=d_matrixid*columnPitch+i*d_veci;假设i=0,线程0访问d_matrix0,线程1访问d_matrixcolumnPitch,线程2访问d_matrix2*columnPitch,这些数据的地址并不相邻，因此没有满足合并访问的要求。columnPitch由函数cudaMallocPitch返回，保证段对齐。,怎样才能使用访问d_matrix满足合并访问要求？,13,矩阵转置,转置
9、后访问d_matrix的模式变成了for(int i=0;i rowSize;i+)temp+=d_matrixi*columnPitch+id*d_veci;假设i=0,线程0访问d_matrix0,线程1访问d_matrix,线程2访问d_matrix2,此时满足合并访问的要求。此时运行时间下降到了4.65ms,性能提高到原来的30多倍,这充分说明了合并访问的重要性。,14,更进一步,for(int i=0;i rowSize;i+)temp+=d_matrixi*columnPitch+id*d_veci;从上面代码很明显的看到d_vec在计算的过程中不变，而且每个线程都访问相同的地址,
10、故可以考虑将它存放在constant中,15,constant优化,static void _global_ mxvNaiveTransposeConstant(int rowSize,int columnSize,int columnPitch,const float*d_matrix,const int start,float*d_r)uint id=blockDim.x*blockIdx.x+threadIdx.x;if(columnSize rowSize?rowSize:start+CONSTANTSIZE;for(int i=start;i end;i+)temp+=d_matri
11、xi*columnPitch+id*c_vi-start;d_rid+=temp;,其中:c_v中constant存储器数组,大小为CONSTANTSIZE。,16,耗时4.17 ms,constant优化(续),问题：如果d_v的大小超过constant的64KB大小限制，怎么办？解决方法：分批，多次传输和启动内核,17,更进一步,很明显,对于block内线程来说,向量都是共享的,因此我们可以使用比constant更快的shared memory来存储,此时相比使用constant，我们免掉了在向量比较大时多次数据拷贝和启动kernel的开销,而且没有使用全局变量,代码的可扩展性更好.由于可
12、能因为shared memory大小存储不了向量,因此需要将向量分块,每次传一小块到shared中,计算完这一小块后,再传一小块接着计算.,18,shared优化,static void _global_ mxvNaiveTransposeShared(int rowSize,int columnSize,int columnPitch,const float*d_matrix,const float*d_v,const int sharedSize,float*d_r)uint id=blockDim.x*blockIdx.x+threadIdx.x;float temp=0.0f;exte
13、rn _shared_ float s_v;for(int start=0;start rowSize?rowSize:start+sharedSize;,19,shared优化(续),#pragma unroll 8 for(int i=start;i end;i+)temp+=d_matrixi*columnPitch+id*s_vi-start;if(id columnSize)d_rid=temp;,20,耗时2.62 ms,矩阵转置的性能,前面的CUDA代码都是基于转置后的矩阵来计算的,因此矩阵转置的性能非常重要,下面的sdk中的transposeNew转置8192*8192的flo
14、at在GTX 295上的数据,21,由于矩阵转置比较慢，因此在很多情况下，我们要使用不转置矩阵的办法,关于block和warp,Block，CUDA线程以block为单位分发到SM上执行，因此使用block线程为单位来处理数据是一个很nature的选择。Warp，block中的线程会以32个为单位划分，这32个线程称为warp,warp中线程的id是连续的，由于SM调度线程的单位是warp，因此在某些情况下，显式的使用warp可获得更好的性能。,22,Block模式,算法:一个block处理矩阵的一行和向量乘积，其中block中的每个线程处理该行中的一个与对应向量元素的乘积,然后归约。stat
15、ic void _global_ mxvBlock(int rowSize,int columnSize,int pitchItem,const float*_restrict_ d_matrix,const float*_restrict_ d_vec,float*_restrict_ d_r)unsigned int tid=threadIdx.x;extern _shared_ float s_r;float temp=0.0f;for(int i=tid;i columnSize;i+=blockDim.x)temp+=d_matrixblockIdx.x*pitchItem+i*d_
16、veci;s_rtid=temp;_syncthreads();/省略归约代码,23,耗时5.42 ms,Warp模式,耗时4.10 ms,24,具体的计算和block模式差不多,只是使用一个warp线程计算矩阵的一行与向量的乘积,在我的测试中发现,这个算法对于行大于列的矩阵效果很好,很多时候性能是block的两倍以上。,cublas,函数:cublasSgemv,25,耗时2.61 ms,总结一下,26,总结一下（续）,矩阵转置以满足合并访问使用常量存储器，共享存储器使用block模式和warp模式,27,其它的一些优化方法手动循环展开数据预取指令混合,感谢itpub提供的这次机会,谢谢大家,欢迎提问!,28,
