cpu實現a*b = c
的矩陣乘法(矩陣尺寸是n*m的,n和m大於1000)
將cpu**移植到cuda。將cpu值傳入gpu,使用cuda計算,與cpu結果對比。
優化思路1:將矩陣分塊進行計算
優化思路2:使用share memory進行優化
優化思路3:將資料繫結在texture上
廢話不多說,直接上原始碼
/* cpumatmultiply:cpu下矩陣乘法
* a:第乙個矩陣指標,表示a[m][n]
* b:第二個矩陣指標,表示b[n][s]
* result:結果矩陣,表示為result[m][s]
*/void cpumatmultiply(const int * a,const int * b, int *result,const int m,const int n,const int s)}}
}
ps:有關cuda的環境搭建,hello_world工程的建立可以移步
下面直接進入正題,矩陣乘法的移植。
從cpu上直接移植矩陣乘法到gpu上是非常簡單的,不需要for迴圈,直接通過cuda執行緒的id號,即threadidx.x
和threadidx.y
即可操作相應的資料。
gpu矩陣乘法核函式-源**:
/* gpumatmultkernel:gpu下矩陣乘法核函式
* a:第乙個矩陣指標,表示a[m][n]
* b:第二個矩陣指標,表示b[n][s]
* result:結果矩陣,表示result[m][s]
*/__global__ void gpumatmultkernel(const int *a, const int *b, int *result, const int m, const int n, const int s)
}}
其中,blockidx
、blockdim
、threadidx
、griddim
都是cuda的內建變數。
blockidx
、threadidx
:分別cuda中線程塊的id、執行緒的id。
blockdim
、griddim
:分別是cuda中線程塊的維度,執行緒網格的維度。
if
語句是為了判斷是否是對矩陣中的資料進行操作,防止在當執行緒threadid
超過了m*s時,使用result[threadid]
出現越界行為。
後續的操作則是矩陣乘法的簡單實現。
cuda c支援共享記憶體。關鍵字__shared__新增到宣告中,這將使這個變數駐留在共享記憶體中。當宣告變數為shared變數時,執行緒塊中的每乙個執行緒都共享這塊記憶體,使得乙個執行緒塊中的多個執行緒能夠在計算上進行通訊和協作。而且,共享記憶體緩衝區駐留在物理gpu上,而不是駐留在gpu之外的系統記憶體中。因此,在訪問共享記憶體時的延遲要遠遠低於訪問普通緩衝區的延遲,提高了效率。
使用共享記憶體的核函式:
/* gpumatmultwithsharedkernel:gpu下使用shared記憶體的矩陣乘法
* a:第乙個矩陣指標,表示a[height_a][width_a]
* b:第二個矩陣指標,表示b[width_a][width_b]
* result:結果矩陣,表示result[height_a][width_b]
*/template
__global__ void gpumatmultwithsharedkernel(const
int *a, const
int *b, int *result, const
int height_a, const
int width_a, const
int width_b)
const
int begin_a = block_y * blockdim.y * width_a;
const
int end_a = begin_a + width_a - 1;
const
int step_a = blockdim.x;
const
int begin_b = block_x * blockdim.x;
const
int step_b = blockdim.y * width_b;
int result_temp = 0;
for (int index_a = begin_a, int index_b = begin_b;
index_a < end_a; index_a += step_a, index_b += step_b)
__syncthreads();
}int begin_result = block_y * blockdim.y * width_b + begin_b;
result[begin_result + thread_y * width_b + thread_x] = result_temp;
}
矩陣乘法的並行運算,每次計算矩陣的一塊資料。利用共享記憶體的共享功能,每次將一塊資料儲存到共享記憶體中使得乙個執行緒塊同時呼叫資料進行計算當前塊相對應得矩陣乘法結果值。
注意:cuda架構將確保,除非執行緒塊中的每個執行緒都執行了__syncthreads()
,否則沒有任何執行緒能執行__syncthreads()
之後的指令。例如:
if (cacheindex < i)
此時,如果有執行緒中的cacheindex>i
則會出現一直等待的情況。
紋理記憶體,cuda c程式中的另一種唯讀記憶體。紋理記憶體是專門為那些在記憶體訪問模式中存在大量空間區域性性的圖形應用程式而設計的。在某個計算應用程式中,這意味著乙個執行緒讀取的位置可能與鄰近執行緒讀取的位置「非常接近」,如下圖所示。
從數學的角度看,這四個位址並非連續的,在一般的cpu快取模式中,這些位址將不會快取。但由於gpu紋理快取是專門為了加速這種訪問模式而設計的,因此如果在這種情況中使用紋理記憶體而不是全域性記憶體,那麼將會獲得效能提公升。
使用紋理記憶體的核函式-源**:
/* gpumatmultwithtexturekernel:gpu下使用texture記憶體的矩陣乘法
* result:結果矩陣,表示為result[m][s];
* m:表示為矩陣a與矩陣result的行數
* n:表示矩陣a的列數,矩陣b的行數
* s:表示矩陣b和矩陣result的列數
*/__global__ void gpumatmultwithtexturekernel(int * result, const
int m, const
int n, const
int s)
result[offset] = temp_result;}}
紋理記憶體的運用與普通記憶體運用時的演算法大致相同,只不過資料是在核函式中呼叫tex1dfetch
從紋理中提取。
在使用紋理記憶體時,主要注意的是紋理記憶體的使用。
首先,需要將輸入的資料宣告為texture型別的引用。
注意,輸入的資料是什麼型別,相應的紋理也應該與之一致。並且紋理引用必須宣告為檔案作用域內的全域性變數。
//這些變數將位於gpu上
texture
texa;
//二維紋理引用,增加了代表維數的引數2
texture
texb;
在為這兩個緩衝區分配了gpu記憶體後,需要通過cudabindtexture
將這些變數繫結到記憶體緩衝區。這相當於告訴cuda執行時兩件事:
cudabindtexture(null, texa, dev_a, desc, m * n * sizeof(int));
cudabindtexture(null, texb, dev_b, desc, n * s * sizeof(int));
在繫結紋理時,cuda執行時要求提供乙個cudachannelformatdesc。此時,需要呼叫cudacreatechanneldesc()
。
最後,通過cudaunbindtexture()
函式來取消紋理的繫結。
矩陣乘法計算中,執行時間上有:共享記憶體 < 紋理記憶體 ≈ 普通gpu移植 < cpu運算。
具體源**
cuda矩陣相乘 CUDA的矩陣乘法
2 那麼下面就是不使用shared memory的並行化演算法的思路。簡單地來說,就是將上述可並行化的部分傳遞給gpu,使用cuda來計算。如下 void matrixmulondevice float m,float n,float p,intwidth int size width width ...
《多核程式設計》學習筆記 矩陣乘法並行化
說到矩陣乘法,最先想到的就是用兩個for迴圈,迴圈矩陣a的行再迴圈矩陣b的列,從而實現矩陣a與b的相乘。1 下面是序列演算法的實現 include includetypedef struct matrix void initial matrix m,int row,int col void init...
cuda實現任意尺寸的矩陣乘法
nvidia gpu常見的塊內線程數最大為1024,當矩陣的行數和列數均小於1024時,我們可以簡單的採用行和列點到點依次相乘構建核函式,即塊內的每個執行緒負責一對元素的乘積計算,然後將所有塊內線程相乘的結果累加求和,得到結果矩陣對應行和列的元素值。code 參照cuda11指導手冊,給出核函式 如...