矩阵与向量乘法的 CUDA 优化

1. 矩阵与向量乘法的CUDA优化 风辰 2010年12月11日 2011年1月8日修订

2. 目的 • 对于CUDA程序开发来说，优化往往是整个开发过程的核心，不同算法，不同存储器组织的程序性能往往差几十倍，本文通过一个简单的例子来展示CUDA开发中一些重要的因素对性能的影响。

3. 假设读者拥有以下知识 • 拥有C语言编程的经验,最好拥有并行编程经验 • 懂得CUDA,最好用CUDA写过代码

4. 测试环境 Intel xeon 5405 2.0 GHz Geforce GTX 295(只使用单核) Gcc 4.3.3 CUDA toolkit 3.1 只测试计算时间,不包括数据传输

5. 符号说明 • matrix:矩阵数据指针,以行为主序或者列为主序存储 • v || vec: 向量指针 • r: 矩阵和向量乘的结果指针 • rowSize: 表示矩阵的行数,也是r的长度 • columnSize:表示矩阵的列数,也是v的长度 • 所有指向显存的指针加前缀d_

6. 编译配置 矩阵尺寸8192*8192　单精度 编译选项-O3 –funroll-loops –msse CPU计时函数采用gettimeofday, clock,GPU计时函数采用CUDA event

7. 串行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 += matrix[i*columnSize+j]*v[j]; } r[i] = re; } } 运行时间120 ms,不使用-O3运行耗时490 ms

8. 简单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(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]; } } 运行时间99ms

9. 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(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]; } 运行时间50ms

10. CUDA优化注意事项 一、选择好的并行方式 选择好的算法，以发掘更多的数据并行性 二、保持SM忙碌 尽量利用所有的SM参与计算，可以通过加大数据量或减小线程块大小达到目的 三、优化存储器利用 保证全局存储器合并访问 使用速度更快的constant或shared存储器

11. CUDA-naïve版本 算法:每个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_matrix[id*columnPitch+i]*d_vec[i]; } d_r[id] = temp; } 耗时150 ms > 串行120ms

12. CUDA-naïve 为什么比串行还慢？ columnPitch的作用是什么？ 访问d_matrix没有满足合并访问的要求 什么是合并访问？

13. 合并访问 一句话：相邻线程访问段对齐的相邻地址 为什么说访问d_matrix没有满足合并访问要求 for(int i = 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满足合并访问要求？

14. 矩阵转置 转置后访问d_matrix的模式变成了 for(int i = 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多倍,这充分说明了合并访问的重要性。

15. 更进一步 for(int i = 0; i < rowSize; i++){ temp += d_matrix[i*columnPitch+id]*d_vec[i]; } 从上面代码很明显的看到d_vec在计算的过程中不变，而且每个线程都访问相同的地址,故可以考虑将它存放在constant中

16. 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 <= id) return; float temp = 0.0f; int end = start+CONSTANTSIZE > rowSize ? rowSize : start+CONSTANTSIZE; for(int i = start; i < end; i++){ temp += d_matrix[i*columnPitch+id]*c_v[i-start]; } d_r[id] += temp;} 其中:c_v中constant存储器数组, 大小为CONSTANTSIZE。 耗时4.17 ms

17. constant优化(续) 问题：如果d_v的大小超过constant的64KB大小限制，怎么办？ 解决方法：分批，多次传输和启动内核

18. 更进一步 很明显, 对于block内线程来说,向量都是共享的,因此我们可以使用比constant更快的shared memory来存储,此时相比使用constant，我们免掉了在向量比较大时多次数据拷贝和启动kernel的开销,而且没有使用全局变量,代码的可扩展性更好. 由于可能因为shared memory大小存储不了向量,因此需要将向量分块,每次传一小块到shared中,计算完这一小块后,再传一小块接着计算.

19. 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; start += sharedSize){ __syncthreads(); #pragma unroll 4 for(int i = threadIdx.x; i < sharedSize&&i+start<rowSize; i += blockDim.x){ s_v[i] = d_v[start+i]; } __syncthreads(); if(columnSize <= id) continue; int end = start+sharedSize > rowSize ? rowSize : start+sharedSize;

20. shared优化(续) #pragma unroll 8 for(int i = start; i < end; i++){ temp += d_matrix[i*columnPitch+id]*s_v[i-start]; } } if(id < columnSize) d_r[id] = temp; } 耗时2.62ms

21. 矩阵转置的性能 前面的CUDA代码都是基于转置后的矩阵来计算的,因此矩阵转置的性能非常重要,下面的sdk中的transposeNew转置8192*8192的float在GTX 295上的数据 由于矩阵转置比较慢，因此在很多情况下，我们要使用不转置矩阵的办法

22. 关于block和warp Block，CUDA线程以block为单位分发到SM上执行，因此使用block线程为单位来处理数据是一个很nature的选择。 Warp，block中的线程会以32个为单位划分，这32个线程称为warp, warp中线程的id是连续的，由于SM调度线程的单位是warp，因此在某些情况下，显式的使用warp可获得更好的性能。

23. 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_matrix[blockIdx.x*pitchItem+i]*d_vec[i]; } s_r[tid] = temp; __syncthreads(); ……//省略归约代码 } } 耗时5.42 ms

24. Warp模式 具体的计算和block模式差不多,只是使用一个warp线程计算矩阵的一行与向量的乘积,在我的测试中发现,这个算法对于行大于列的矩阵效果很好,很多时候性能是block的两倍以上。 耗时4.10 ms

25. cublas 函数: cublasSgemv 耗时2.61 ms

26. 总结一下

27. 总结一下（续） • 矩阵转置以满足合并访问 • 使用常量存储器，共享存储器 • 使用block模式和warp模式 • 其它的一些优化方法 • 手动循环展开 • 数据预取 • 指令混合

28. 感谢itpub提供的这次机会,谢谢大家,欢迎提问!