矩阵与向量乘法的优化

矩阵与向量乘法的优化第1页，课件共28页，创作于2023年2月目的对于CUDA程序开发来说，优化往往是整个开发过程的核心，不同算法，不同存储器组织的程序性能往往差几十倍，本文通过一个简单的例子来展示CUDA开发中一些重要的因素对性能的影响。2第2页，课件共28页，创作于2023年2月假设读者拥有以下知识拥有C语言编程的经验,最好拥有并行编程经验懂得CUDA,最好用CUDA写过代码3第3页，课件共28页，创作于2023年2月测试环境Intelxeon54052.0GHzGeforceGTX295(只使用单核)Gcc4.3.3CUDAtoolkit3.1只测试计算时间,不包括数据传输4第4页，课件共28页，创作于2023年2月符号说明matrix:矩阵数据指针,以行为主序或者列为主序存储v||vec:向量指针r:矩阵和向量乘的结果指针rowSize:表示矩阵的行数,也是r的长度columnSize:表示矩阵的列数,也是v的长度所有指向显存的指针加前缀d_5第5页，课件共28页，创作于2023年2月编译配置矩阵尺寸8192*8192单精度编译选项-O3–funroll-loops–msseCPU计时函数采用gettimeofday,clock,GPU计时函数采用CUDAevent6第6页，课件共28页，创作于2023年2月串行C版本算法:遍历矩阵行,每行和向量相乘,最终结果为一向量voidmxv(constintrowSize,constintcolumnSize,constfloat*matrix,constfloat*v,float*r){for(inti=0;i<rowSize;i++){floatre=0.0f;for(intj=0;j<columnSize;j++){re+=matrix[i*columnSize+j]*v[j];}r[i]=re;}}7运行时间120ms,不使用-O3运行耗时490ms第7页，课件共28页，创作于2023年2月简单SSE版本算法:利用sse指令计算矩阵每行和向量的乘积voidmxvSSE(constintrowSize,constintcolumnSize,constfloat*matrix,constfloat*v,float*r){ __m128*mv=(__m128*)v;__m128*mm=(__m128*)matrix;for(inti=0;i<rowSize;i++){__m128re=_mm_set_ps(0.0f,0.0f,0.0f,0.0f);for(intj=0;j<columnSize/4;j++){re=_mm_add_ps(re,_mm_mul_ps(mm[i*columnSize/4+j],mv[j]));} float__attribute((aligned(16)))a[4]; _mm_store_ps(a,re);r[i]=a[0]+a[1]+a[2]+a[3];}}运行时间99ms8第8页，课件共28页，创作于2023年2月SSE+openmp算法:使用二线程并行计算行循环voidmxvSSEOpenmp(constintrowSize,constintcolumnSize,float*matrix,float*vec,float*r){ __m128*mv=(__m128*)v;__m128*mm=(__m128*)matrix;#pragmaompparallelfornum_threads(2) for(inti=0;i<rowSize;i++){ __m128re=_mm_set_ps(0.0f,0.0f,0.0f,0.0f); for(intj=0;j<columnSize/4;j++){ re=_mm_add_ps(re,_mm_mul_ps(mm[i*columnSize/4+j],mv[j]));}float__attribute((aligned(16)))a[4];_mm_store_ps(a,re);r[i]=a[0]+a[1]+a[2]+a[3];}运行时间50ms9第9页，课件共28页，创作于2023年2月CUDA优化注意事项一、选择好的并行方式

选择好的算法，以发掘更多的数据并行性二、保持SM忙碌

尽量利用所有的SM参与计算，可以通过加大数据量或减小线程块大小达到目的三、优化存储器利用

保证全局存储器合并访问

使用速度更快的constant或shared存储器10第10页，课件共28页，创作于2023年2月CUDA-naïve版本算法:每个CUDA线程计算矩阵的一行与向量乘积staticvoid__global__mxvNaive(introwSize,intcolumnSize,intcolumnPitch,constfloat*d_matrix,constfloat*d_vec,float*d_r){ uintid=blockDim.x*blockIdx.x+threadIdx.x;if(rowSize<=id)return;floattemp=0.0f; #pragmaunroll4 for(inti=0;i<columnSize;i++){temp+=d_matrix[id*columnPitch+i]*d_vec[i];}d_r[id]=temp;}耗时150ms>串行120ms11第11页，课件共28页，创作于2023年2月CUDA-naïve为什么比串行还慢？columnPitch的作用是什么？访问d_matrix没有满足合并访问的要求什么是合并访问？12第12页，课件共28页，创作于2023年2月合并访问一句话：相邻线程访问段对齐的相邻地址为什么说访问d_matrix没有满足合并访问要求 for(inti=0;i<columnSize;i++){temp+=d_matrix[id*columnPitch+i]*d_vec[i];}假设i=0,线程0访问d_matrix[0],线程1访问d_matrix[columnPitch],线程2访问d_matrix[2*columnPitch],这些数据的地址并不相邻，因此没有满足合并访问的要求。columnPitch由函数cudaMallocPitch返回，保证段对齐。怎样才能使用访问d_matrix满足合并访问要求？13第13页，课件共28页，创作于2023年2月矩阵转置转置后访问d_matrix的模式变成了for(inti=0;i<rowSize;i++){temp+=d_matrix[i*columnPitch+id]*d_vec[i];}

假设i=0,线程0访问d_matrix[0],线程1访问d_matrix[１],线程2访问d_matrix[2],此时满足合并访问的要求。此时运行时间下降到了4.65ms,性能提高到原来的30多倍,这充分说明了合并访问的重要性。14第14页，课件共28页，创作于2023年2月更进一步for(inti=0;i<rowSize;i++){temp+=d_matrix[i*columnPitch+id]*d_vec[i];}从上面代码很明显的看到d_vec在计算的过程中不变，而且每个线程都访问相同的地址,故可以考虑将它存放在constant中15第15页，课件共28页，创作于2023年2月constant优化staticvoid__global__mxvNaiveTransposeConstant(introwSize,intcolumnSize,intcolumnPitch,constfloat*d_matrix,constintstart,float*d_r){uintid=blockDim.x*blockIdx.x+threadIdx.x;if(columnSize<=id)return;floattemp=0.0f;intend=start+CONSTANTSIZE>rowSize?rowSize:start+CONSTANTSIZE;for(inti=start;i<end;i++){temp+=d_matrix[i*columnPitch+id]*c_v[i-start];}d_r[id]+=temp;}其中:

c_v中constant存储器数组,大小为CONSTANTSIZE。16耗时4.17ms第16页，课件共28页，创作于2023年2月constant优化(续)问题：如果d_v的大小超过constant的64KB大小限制，怎么办？解决方法：分批，多次传输和启动内核17第17页，课件共28页，创作于2023年2月更进一步很明显,对于block内线程来说,向量都是共享的,因此我们可以使用比constant更快的sharedmemory来存储,此时相比使用constant，我们免掉了在向量比较大时多次数据拷贝和启动kernel的开销,而且没有使用全局变量,代码的可扩展性更好.由于可能因为sharedmemory大小存储不了向量,因此需要将向量分块,每次传一小块到shared中,计算完这一小块后,再传一小块接着计算.18第18页，课件共28页，创作于2023年2月shared优化staticvoid__global__mxvNaiveTransposeShared(introwSize,intcolumnSize,intcolumnPitch,constfloat*d_matrix,constfloat*d_v,constintsharedSize,float*d_r){ uintid=blockDim.x*blockIdx.x+threadIdx.x; floattemp=0.0f; extern__shared__floats_v[]; for(intstart=0;start<rowSize;start+=sharedSize){ __syncthreads(); #pragmaunroll4 for(inti=threadIdx.x;i<sharedSize&&i+start<rowSize;i+=blockDim.x){ s_v[i]=d_v[start+i]; }__syncthreads(); if(columnSize<=id)continue; intend=start+sharedSize>rowSize?rowSize:start+sharedSize;19第19页，课件共28页，创作于2023年2月shared优化(续) #pragmaunroll8 for(inti=start;i<end;i++){ temp+=d_matrix[i*columnPitch+id]*s_v[i-start]; } } if(id<columnSize) d_r[id]=temp;}20耗时2.62

ms第20页，课件共28页，创作于2023年2月矩阵转置的性能前面的CUDA代码都是基于转置后的矩阵来计算的,因此矩阵转置的性能非常重要,下面的sdk中的transposeNew转置8192*8192的float在GTX295上的数据21方法说明吞吐量Kernel运行时间transposeNew-Outer-fine-grained67.7686GB/s7.37804stransposeNew-Inner-fine-grained72.7973GB/s6.86839stransposeNew-Outer-diagonaltranspose28.4115GB/s17.59853stransposeNew-Inner-diagonaltranspose33.8458GB/s14.77287stransposeNew-Outer-nobankconflicttrans17.2629GB/s28.96379stransposeNew-Inner-nobankconflicttranss17.0058GB/s29.40170由于矩阵转置比较慢，因此在很多情况下，我们要使用不转置矩阵的办法第21页，课件共28页，创作于2023年2月关于block和warpBlock，CUDA线程以block为单位分发到SM上执行，因此使用block线程为单位来处理数据是一个很nature的选择。Warp，block中的线程会以32个为单位划分，这32个线程称为warp,warp中线程的id是连续的，由于SM调度线程的单位是warp，因此在某些情况下，显式的使用warp可获得更好的性能。22第22页，课件共28页，创作于2023年2月Block模式算法:一个block处理矩阵的一行和向量乘积，其中block中的每个线程处理该行中的一个与对应向量元素的乘积,然后归约。staticvoid__global__mxvBlock(introwSize,intcolumnSize,intpitchItem,constfloat*__restrict__d_matrix,constfloat*__restrict__d_vec,float*__restrict__d_r){ unsignedinttid=threadIdx.x; extern__shared__floats_r[]; floattemp=0.0f; for(inti=tid;i<columnSize;i+=blockDim.x){ temp+=d_matrix[blockIdx.x*pitchItem+i]*d_vec[i]; } s_r[tid]=temp;__syncthreads(); ……//省略归约代码 }}23耗时5.4