矩阵与向量乘法的CUDA优化.ppt

矩阵与向量乘法的CUDA优化,风辰 2010年12月11日 2011年1月8日修订,1,目的,对于CUDA程序开发来说，优化往往是整个开发过程的核心，不同算法，不同存储器组织的程序性能往往差几十倍，本文通过一个简单的例子来展示CUDA开发中一些重要的因素对性能的影响。,2,假设读者拥有以下知识,拥有C语言编程的经验,最好拥有并行编程经验 懂得CUDA,最好用CUDA写过代码,3,测试环境,Intel xeon 5405 2.0 GHz Geforce GTX 295(只使用单核) Gcc 4.3.3 CUDA toolkit 3.1 只测试计算时间,不包括数据传输,4,符号说明,matrix:矩阵数据指针,以行为主序或者列为主序存储 v | vec: 向量指针 r: 矩阵和向量乘的结果指针 rowSize: 表示矩阵的行数,也是r的长度 columnSize:表示矩阵的列数,也是v的长度 所有指向显存的指针加前缀d_,5,编译配置,矩阵尺寸8192*8192 单精度 编译选项-O3 funroll-loops msse CPU计时函数采用gettimeofday, clock,GPU计时函数采用CUDA event,6,串行C版本,算法:遍历矩阵行,每行和向量相乘,最终结果为一向量 void mxv(const int rowSize, const int columnSize, const float *matrix, const 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 = (_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 rowSize, 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(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, 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没有满足合并访问的要求 什么是合并访问？,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,矩阵转置,转置后访问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在计算的过程中不变，而且每个线程都访问相同的地址,故可以考虑将它存放在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_matrixi*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的开销,而且没有使用全局变量,代码的可扩展性更好. 由于可能因为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; extern _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的float在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中的每个线程处理该行中的一个与对应向量元素的乘积,然后归约。 static 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_veci; s_rtid = temp; _syncthreads(); /省略归约代码 ,23,耗时5.42 ms,Warp模式