矩陣與向量乘法的CUDA優(yōu)化ppt課件

1、矩陣與向量乘法的CUDA優(yōu)化風(fēng)辰2021年12月11日2021年1月8日修訂目的 對于CUDA程序開發(fā)來說，優(yōu)化往往是整個(gè)開發(fā)過程的中心，不同算法，不同存儲(chǔ)器組織的程序性能往往差幾十倍，本文經(jīng)過一個(gè)簡單的例子來展現(xiàn)CUDA開發(fā)中一些重要的要素對性能的影響。假設(shè)讀者擁有以下知識(shí)l擁有C言語編程的閱歷,最好擁有并行編程閱歷l懂得CUDA,最好用CUDA寫過代碼測試環(huán)境Intel xeon 5405 2.0 GHzGeforce GTX 295(只運(yùn)用單核)Gcc 4.3.3 CUDA toolkit 3.1只測試計(jì)算時(shí)間,不包括數(shù)據(jù)傳輸符號闡明 matrix:矩陣數(shù)據(jù)指針,以行為主序或者列為主序存

2、儲(chǔ) v | vec: 向量指針 r: 矩陣和向量乘的結(jié)果指針 rowSize: 表示矩陣的行數(shù),也是r的長度 columnSize:表示矩陣的列數(shù),也是v的長度 一切指向顯存的指針加前綴d_編譯配置矩陣尺寸8192*8192單精度編譯選項(xiàng)-O3 funroll-loops msseCPU計(jì)時(shí)函數(shù)采用gettimeofday, clock,GPU計(jì)時(shí)函數(shù)采用CUDA event串行C版本算法:遍歷矩陣行,每行和向量相乘,最終結(jié)果為一向量void mxv(const int rowSize, const int columnSize, const float *matrix, const floa

3、t *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; 運(yùn)轉(zhuǎn)時(shí)間運(yùn)轉(zhuǎn)時(shí)間120 ms,120 ms,不運(yùn)用不運(yùn)用-O3-O3運(yùn)轉(zhuǎn)耗時(shí)運(yùn)轉(zhuǎn)耗時(shí)490 ms490 ms簡單SSE版本算法算法: :利用利用ssesse指令計(jì)算矩陣每行和向量的乘積指令計(jì)算矩陣每行和向量的乘積void mxvSSE(const int rowSize, const int void mxvSSE(const

4、 int rowSize, const int columnSize, const float columnSize, const float * *matrix, const matrix, const 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 rowS

5、ize; i+) _m128 re = _mm_set_ps(0.0f, 0.0f, _m128 re = _mm_set_ps(0.0f, 0.0f, 0.0f, 0.0f);0.0f, 0.0f); for(int j = 0; j columnSize/4; for(int j = 0; j columnSize/4; j+)j+) re = _mm_add_ps(re, re = _mm_add_ps(re, _mm_mul_ps(mmi_mm_mul_ps(mmi* *columnSize/4+j, mvj);columnSize/4+j, mvj); float _attribut

6、e(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; 運(yùn)轉(zhuǎn)時(shí)間99msSSE + openmp算法算法: :運(yùn)用二線程并行計(jì)算行循環(huán)運(yùn)用二線程并行計(jì)算行循環(huán)void mxvSSEOpenmp(const int rowSize, const int void mxvSSEOpenmp(const int rowSize, const int columnSize, flo

7、at columnSize, float * *matrix, float matrix, float * *vec, float vec, float * *r)r)_m128 _m128 * *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(in

8、t i = 0; i rowSize; i+) _m128 re = _mm_set_ps(0.0f, 0.0f, _m128 re = _mm_set_ps(0.0f, 0.0f, 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, re = _mm_add_ps(re, _mm_mul_ps(mmi_mm_mul_ps(mmi* *columnSize/4+j, mvj);columnSize/4+j, mvj);

9、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; 運(yùn)轉(zhuǎn)時(shí)間50msCUDA優(yōu)化本卷須知一、選擇好的并行方式選擇好的算法，以開掘更多的數(shù)據(jù)并行性二、堅(jiān)持SM忙碌盡量利用一切的SM參與計(jì)算，可以經(jīng)過加大數(shù)據(jù)量或減小線程塊大小到達(dá)目的三、優(yōu)化存儲(chǔ)器利用保證全局存儲(chǔ)器合并訪問運(yùn)用速度更快的constant或shared存儲(chǔ)器CUDA-nave版

10、本算法算法: :每個(gè)每個(gè)CUDACUDA線程計(jì)算矩陣的一行與向量乘積線程計(jì)算矩陣的一行與向量乘積static void _global_ mxvNaive(int rowSize, static void _global_ mxvNaive(int rowSize, int columnSize, int columnPitch, const float int columnSize, int columnPitch, const float * *d_matrix, const float d_matrix, const float * *d_vec, float d_vec, float

11、* *d_r) d_r) uint id = blockDim.xuint id = blockDim.x* *blockIdx.x + blockIdx.x + threadIdx.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 串行120msCUDA-nave為什么比串行還

12、慢？為什么比串行還慢？ columnPitch columnPitch的作用是什么？的作用是什么？訪問訪問d_matrixd_matrix沒有滿足合并訪問的要求沒有滿足合并訪問的要求什么是合并訪問？什么是合并訪問？合并訪問一句話：相鄰線程訪問段對齊的相鄰地址為什么說訪問d_matrix沒有滿足合并訪問要求for(int i = 0; i columnSize; i+) temp += d_matrixid*columnPitch+i*d_veci; 假設(shè)i=0, 線程0訪問d_matrix0,線程1訪問d_matrixcolumnPitch,線程2訪問d_matrix2*columnPitch

13、,這些數(shù)據(jù)的地址并不相鄰，因此沒有滿足合并訪問的要求。columnPitch由函數(shù)cudaMallocPitch前往，保證段對齊。怎樣才干運(yùn)用訪問d_matrix滿足合并訪問要求？矩陣轉(zhuǎn)置轉(zhuǎn)置后訪問d_matrix的方式變成了for(int i = 0; i rowSize; i+) temp += d_matrixi*columnPitch+id*d_veci;假設(shè)i=0, 線程0訪問d_matrix0,線程1訪問d_matrix,線程2訪問d_matrix2,此時(shí)滿足合并訪問的要求。此時(shí)運(yùn)轉(zhuǎn)時(shí)間下降到了4.65ms,性能提高到原來的30多倍,這充分闡明了合并訪問的重要性。更進(jìn)一步for(i

14、nt i = 0; i rowSize; i+) temp += d_matrixi*columnPitch+id*d_veci;從上面代碼很明顯的看到d_vec在計(jì)算的過程中不變，而且每個(gè)線程都訪問一樣的地址,故可以思索將它存放在constant中constant優(yōu)化static void _global_ mxvNaiveTransposeConstant(int rowSize, int columnSize, int columnPitch, const float *d_matrix, const int start, float *d_r) uint id = blockDim.x

15、*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存儲(chǔ)器數(shù)組, 大小為CONSTANTSIZE。耗時(shí)4.17 msconstant優(yōu)化(續(xù))問題：假設(shè)d_v的大小超越constant的64KB大小限制，怎樣辦？處理方法：分批，多次傳輸和啟動(dòng)內(nèi)核更進(jìn)一步很明顯, 對于block內(nèi)

16、線程來說,向量都是共享的,因此我們可以運(yùn)用比constant更快的shared memory來存儲(chǔ),此時(shí)相比運(yùn)用constant，我們免掉了在向量比較大時(shí)多次數(shù)據(jù)拷貝和啟動(dòng)kernel的開銷,而且沒有運(yùn)用全局變量,代碼的可擴(kuò)展性更好.由于能夠由于shared memory大小存儲(chǔ)不了向量,因此需求將向量分塊,每次傳一小塊到shared中,計(jì)算完這一小塊后,再傳一小塊接著計(jì)算.shared優(yōu)化static void _global_ mxvNaiveTransposeShared(int rowSize, int columnSize, int columnPitch, const float

17、*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 +=

18、blockDim.x) s_vi = d_vstart+i; _syncthreads(); if(columnSize rowSize ? rowSize : start+sharedSize;shared優(yōu)化(續(xù))#pragma unroll 8 for(int i = start; i end; i+) temp += d_matrixi*columnPitch+id*s_vi-start; if(id columnSize) d_rid = temp;耗時(shí)2.62 ms矩陣轉(zhuǎn)置的性能前面的CUDA代碼都是基于轉(zhuǎn)置后的矩陣來計(jì)算的,因此矩陣轉(zhuǎn)置的性能非常重要,下面的sdk中的transp

19、oseNew轉(zhuǎn)置8192*8192的float在GTX 295上的數(shù)據(jù)方法說明方法說明吞吐量吞吐量Kernel運(yùn)行時(shí)間運(yùn)行時(shí)間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-Ou

20、ter-no bank conflict trans17.2629 GB/s28.96379 stransposeNew-Inner-no bank conflict transs17.0058 GB/s29.40170 由于矩陣轉(zhuǎn)置比較慢，因此在很多情況下，我們要運(yùn)用不轉(zhuǎn)置矩陣的方法關(guān)于block和warpBlock，CUDA線程以block為單位分發(fā)到SM上執(zhí)行，因此運(yùn)用block線程為單位來處置數(shù)據(jù)是一個(gè)很nature的選擇。Warp，block中的線程會(huì)以32個(gè)為單位劃分，這32個(gè)線程稱為warp, warp中線程的id是延續(xù)的，由于SM調(diào)度線程的單位是warp，因此在某些情況下，顯式

21、的運(yùn)用warp可獲得更好的性能。Block方式算法:一個(gè)block處置矩陣的一行和向量乘積，其中block中的每個(gè)線程處置該行中的一個(gè)與對應(yīng)向量元素的乘積,然后歸約。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