矩阵与向量乘法的CUDA优化Stillwatersrundeep.流静水深流静水深,人静心深人静心深Wherethereislife,thereishope。有生命必有希望。有生命必有希望目的对于CUDA程序开发来说，优化往往是整个开发过程的核心，不同算法，不同存储器组织的程序性能往往差几十倍，本文通过一个简单的例子来展示CUDA开发中一些重要的因素对性能的影响。2假设读者拥有以下知识l拥有C语言编程的经验,最好拥有并行编程经验l懂得CUDA,最好用CUDA写过代码3测试环境Intel xeon 5405 2.0 GHzGeforce 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 msseCPU计时函数采用gettimeofday, clock,GPU计时函数采用CUDA event6串行C版本算法:遍历矩阵行,每行和向量相乘,最终结果为一向量void mxv(const int rowSize, const int columnSize,void mxv(const int rowSize, const int columnSize, const float *matrix, const float *v, float *r) const float *matrix, const float *v, float *r) for(int i = 0; i rowSize; i+) for(int i = 0; i rowSize; i+) float re = 0.0f; float re = 0.0f; for(int j = 0; j columnSize; j+) for(int j = 0; j columnSize; j+) re += matrixi*columnSize+j*vj; re += matrixi*columnSize+j*vj; ri = re; ri = re; 7运行时间运行时间120 ms,120 ms,不使用不使用-O3-O3运行耗时运行耗时490 ms490 ms简单SSE版本算法算法: :利用利用ssesse指令计算矩阵每行和向量的乘积指令计算矩阵每行和向量的乘积void mxvSSE(const int rowSize, const int columnSize, const float *matrix, const void mxvSSE(const int rowSize, const int columnSize, const float *matrix, const float *v, float *r)float *v, float *r) _m128 *mv = (_m128*)v;_m128 *mv = (_m128*)v; _m128 *mm = (_m128*)matrix; _m128 *mm = (_m128*)matrix; for(int i = 0; i rowSize; i+) for(int i = 0; i rowSize; i+) _m128 re = _mm_set_ps(0.0f, 0.0f, 0.0f, 0.0f); _m128 re = _mm_set_ps(0.0f, 0.0f, 0.0f, 0.0f); for(int j = 0; j columnSize/4; j+) for(int j = 0; j columnSize/4; j+) re = _mm_add_ps(re, _mm_mul_ps(mmi*columnSize/4+j, mvj); re = _mm_add_ps(re, _mm_mul_ps(mmi*columnSize/4+j, mvj); float _attribute(aligned(16) a4;float _attribute(aligned(16) a4; _mm_store_ps(a, re); _mm_store_ps(a, re); ri = a0 + a1 + a2 + a3; ri = a0 + a1 + a2 + a3; 运行时间99ms8SSE + openmp算法算法: :使用二线程并行计算行循环使用二线程并行计算行循环void mxvSSEOpenmp(const int rowSize, const int columnSize, float *matrix, float void mxvSSEOpenmp(const int rowSize, const int columnSize, float *matrix, float *vec, float *r)*vec, float *r)_m128 *mv = (_m128*)v;_m128 *mv = (_m128*)v; _m128 *mm = (_m128*)matrix; _m128 *mm = (_m128*)matrix;#pragma omp parallel for num_threads(2)#pragma omp parallel for num_threads(2)for(int i = 0; i rowSize; i+)for(int i = 0; i rowSize; i+) _m128 re = _mm_set_ps(0.0f, 0.0f, 0.0f, 0.0f);_m128 re = _mm_set_ps(0.0f, 0.0f, 0.0f, 0.0f); for(int j = 0; j columnSize/4; j+)for(int j = 0; j columnSize/4; j+) re = _mm_add_ps(re, _mm_mul_ps(mmi*columnSize/4+j, mvj); re = _mm_add_ps(re, _mm_mul_ps(mmi*columnSize/4+j, mvj); float _attribute(aligned(16) a4; float _attribute(aligned(16) a4; _mm_store_ps(a, re); _mm_store_ps(a, re); ri = a0 + a1 + a2 + a3; ri = a0 + a1 + a2 + a3; 运行时间50ms9CUDA优化注意事项一、选择好的并行方式选择好的算法，以发掘更多的数据并行性二、保持SM忙碌尽量利用所有的SM参与计算，可以通过加大数据量或减小线程块大小达到目的三、优化存储器利用保证全局存储器合并访问使用速度更快的constant或shared存储器10CUDA-nave版本算法算法: :每个每个CUDACUDA线程计算矩阵的一行与向量乘积线程计算矩阵的一行与向量乘积static void _global_ mxvNaive(int rowSize, int columnSize, int columnPitch, const static void _global_ mxvNaive(int rowSize, int columnSize, int columnPitch, const float *d_matrix, const float *d_vec, float *d_r) float *d_matrix, const float *d_vec, float *d_r) uint id = blockDim.x*blockIdx.x + threadIdx.x;uint id = blockDim.x*blockIdx.x + threadIdx.x; if(rowSize = id) return; if(rowSize = id) return; float temp = 0.0f; float temp = 0.0f;#pragma unroll 4 #pragma unroll 4 for(int i = 0; i columnSize; i+)for(int i = 0; i 串行120ms11CUDA-nave为什么比串行还慢？为什么比串行还慢？ columnPitch columnPitch的作用是什么？的作用是什么？访问访问d_matrixd_matrix没有满足合并访问的要求没有满足合并访问的要求什么是合并访问？什么是合并访问？12合并访问一句话：相邻线程访问段对齐的相邻地址为什么说访问d_matrix没有满足合并访问要求for(int i = 0; i columnSize; i+) temp += d_matrixid*columnPitch+i*d_veci; 假设假设i=0, i=0, 线程线程0 0访问访问d_matrix0,d_matrix0,线程线程1 1访问访问d_matrixcolumnPitch,d_matrixcolumnPitch,线程线程2 2访问访问d_matrixd_matrix2*columnPitch,2*columnPitch,这些数据的地址并不相邻，因此没有满足合并访问的要求。这些数据的地址并不相邻，因此没有满足合并访问的要求。columnPitchcolumnPitch由函数由函数cudaMallocPitchcudaMallocPitch返回，保证段对齐。返回，保证段对齐。怎样才能使用访问d_matrix满足合并访问要求？13矩阵转置转置后访问d_matrix的模式变成了for(int i = 0; i rowSize; i+) temp += d_matrixi*columnPitch+id*d_veci;假设假设i=0, i=0, 线程线程0 0访问访问d_matrix0,d_matrix0,线程线程1 1访问访问d_matrixd_matrix,线程线程2 2访问访问d_md_matrix2,atrix2,此时满足合并访问的要求。此时满足合并访问的要求。此时运行时间下降到了此时运行时间下降到了4.65ms4.65ms, ,性能提高到原来的性能提高到原来的3030多倍多倍, ,这充这充分说明了合并访问的重要性。分说明了合并访问的重要性。14更进一步for(int i = 0; i rowSize; i+) temp += d_matrixi*columnPitch+id*d_veci;从上面代码很明显的看到d_vec在计算的过程中不变，而且每个线程都访问相同的地址,故可以考虑将它存放在constant中15constant优化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 msconstant优化(续)问题：如果d_v的大小超过constant的64KB大小限制，怎么办？解决方法：分批，多次传输和启动内核17更进一步很明显, 对于block内线程来说,向量都是共享的,因此我们可以使用比constant更快的shared memory来存储,此时相比使用constant，我们免掉了在向量比较大时多次数据拷贝和启动kernel的开销,而且没有使用全局变量,代码的可扩展性更好.由于可能因为shared memory大小存储不了向量,因此需要将向量分块,每次传一小块到shared中,计算完这一小块后,再传一小块接着计算.18shared优化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; start += sharedSize) _syncthreads();#pragma unroll 4 for(int i = threadIdx.x; i sharedSize&i+startrowSize; i += blockDim.x) s_vi = d_vstart+i; _syncthreads(); if(columnSize rowSize ? rowSize : start+sharedSize;19shared优化(续)#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方法说明方法说明吞吐量吞吐量Kernel运行时间运行时间transposeNew-Outer-fine-grained67.7686 GB/s7.37804 stransposeNew-Inner-fine-grained72.7973 GB/s6.86839 stransposeNew-Outer-diagonal transpose28.4115 GB/s17.59853 stransposeNew-Inner-diagonal transpose33.8458 GB/s14.77287 stransposeNew-Outer-no bank conflict trans17.2629 GB/s28.96379 stransposeNew-Inner-no bank conflict transs17.0058 GB/s29.40170 由于矩阵转置比较慢，因此在很多情况下，我们要使用不转置矩阵的办法关于block和warpBlock，CUDA线程以block为单位分发到SM上执行，因此使用block线程为单位来处理数据是一个很nature的选择。Warp，block中的线程会以32个为单位划分，这32个线程称为warp, warp中线程的id是连续的，由于SM调度线程的单位是warp，因此在某些情况下，显式的使用warp可获得更好的性能。22Block模式算法:一个block处理矩阵的一行和向量乘积，其中block中的每个线程处理该行中的一个与对应向量元素的乘积,然后归约。static void _global_ mxvBlock(int rowSize, int columnSize, int pitchItem, const static void _global_ mxvBlock(int rowSize, int columnSize, int pitchItem, const float* _restrict_ d_matrix,const float* _restrict_ d_vec, float* float* _restrict_ d_matrix,const float* _restrict_ d_vec, float* _restrict_ d_r)_restrict_ d_r)unsigned int tid = threadIdx.x;unsigned int tid = threadIdx.x;extern _shared_ float s_r;extern _shared_ float s_r;float temp = 0.0f;float temp = 0.0f;for(int i = tid; i columnSize; i += blockDim.x)for(int i = tid; i columnSize; i += blockDim.x)temp += d_matrixblockIdx.x*pitchItem+i*d_veci;temp += d_matrixblockIdx.x*pitchItem+i*d_veci; s_rtid = temp; _syncthreads();s_rtid = temp; _syncthreads();/省略归约代码省略归约代码 23耗时5.42 msWarp模式耗时4.10 ms24具体的计算和具体的计算和blockblock模式差不多模式差不多, ,只是使用一个只是使用一个warpwarp线程计算矩阵的一行与向量的乘积线程计算矩阵的一行与向量的乘积, ,在我的测试中在我的测试中发现发现, ,这个算法对于行大于列的矩阵效果很好这个算法对于行大于列的矩阵效果很好, ,很多时很多时候性能是候性能是blockblock的两倍以上。的两倍以上。cublas函数: cublasSgemv25耗时2.61 ms总结一下函数名函数名说明说明时间时间/ms加速比加速比mxv串行C120mxvSSE串行C+SSE991.2mxvSSEOpenmp串行C+SSE+openmp502.4mxvNaive1500.8mxvNaiveTranspose矩阵转置4.626.1mxvNaiveTransposeConstant矩阵转置+constant memory4.228.6mxvNaiveTransposeShared矩阵转置+shared memory2.646.2mxvBlockblock模式5.422.2mxvWarpwarp模式4.129.3cublas调用cublasSgemv函数2.646.226总结一下（续）矩阵转置以满足合并访问使用常量存储器，共享存储器使用block模式和warp模式27其它的一些优化方法其它的一些优化方法l手动循环展开l数据预取l指令混合感谢itpub提供的这次机会,谢谢大家,欢迎提问!28