矩阵和向量乘法CUDA优化

1、矩阵和向量乘法CUDA优化目的 对于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: 矩阵和

2、向量乘的结果指针 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, co

3、nst float const float * *matrix, const float matrix, const float * *v, float v, float * *r)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 re += matrixi* *columnSize+jcolu

4、mnSize+j* *vj;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 void mxvSSE(const int rowSize, const int columnSize, const float * *matrix, const matrix, co

5、nst float float * *v, float v, float * *r)r) _m128 _m128 * *mv = (_m128mv = (_m128* *)v;)v; _m128 _m128 * *mm = (_m128mm = (_m128* *)matrix; )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.0

6、f); 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 re = _mm_add_ps(re, _mm_mul_ps(mmi* *columnSize/4+j, mvj);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 =

7、 a0 + a1 + a2 + a3; ri = a0 + a1 + a2 + a3; 运行时间99ms8SSE + openmp算法算法: :使用二线程并行计算行循环使用二线程并行计算行循环void mxvSSEOpenmp(const int rowSize, const int columnSize, float void mxvSSEOpenmp(const int rowSize, const int columnSize, float * *matrix, float matrix, float * *vec, float vec, float * *r)r)_m128 _m128

8、 * *mv = (_m128mv = (_m128* *)v;)v; _m128 _m128 * *mm = (_m128mm = (_m128* *)matrix;)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.0

9、f, 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 re = _mm_add_ps(re, _mm_mul_ps(mmi* *columnSize/4+j, mvj);columnSize/4+j, mvj); float _attribute(aligned(16) a4; float _attribute(aligned(16) a4; _mm_store_ps(a, re); _mm_sto

10、re_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,

11、int columnSize, int columnPitch, static void _global_ mxvNaive(int rowSize, int columnSize, int columnPitch, const float const float * *d_matrix, const float d_matrix, const float * *d_vec, float d_vec, float * *d_r) d_r) uint id = blockDim.xuint id = blockDim.x* *blockIdx.x + threadIdx.x;blockIdx.x

12、 + 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没有满足合并访问的要求没有满足合并访问的

13、要求什么是合并访问？什么是合并访问？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_matrix2d_matrix2* *columnPitch,columnPitch,这些数据的地址并

14、不相邻，因此没有满足合并这些数据的地址并不相邻，因此没有满足合并访问的要求。访问的要求。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_ma

15、trixd_matrix,线程线程2 2访问访问d_matrix2,d_matrix2,此时满足合并访问的要求。此时满足合并访问的要求。此时运行时间下降到了此时运行时间下降到了4.65ms,4.65ms,性能提高到原来的性能提高到原来的3030多倍多倍, ,这充分说明了合并访问的重要性。这充分说明了合并访问的重要性。14更进一步for(int i = 0; i rowSize; i+) temp += d_matrixi*columnPitch+id*d_veci;从上面代码很明显的看到d_vec在计算的过程中不变，而且每个线程都访问相同的地址,故可以考虑将它存放在constant中15con

16、stant优化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

17、*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的开销,而且没有使用全局变量,代码的可扩展性更好.由于可能因为s

18、hared 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 =

19、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 unr

20、oll 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 st

21、ransposeNew-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 由于矩阵转

22、置比较慢，因此在很多情况下，我们要使用不转置矩阵的办法关于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 floatfloat* * _restrict_ d_matrix,const float _restrict_ d_matrix,const float* * _restrict_ d_vec, float _restrict_ d_vec, float* * _restrict_ _restrict_ d_r)d_r)un