CUDA 版本矩陣乘

URNOTJANET發表於2017-12-27

矩陣

說明：

1.轉載請聯絡本人

2.程式碼在最後

問題描述

矩陣乘法 C = aAB + bC
其中a,b為常數,A,B,C為矩陣

實驗要求

根據記憶體大小測不同規模矩陣的處理速度（GFLOPS/s），並給出計算公式。
請計算系統的理論峰值，如果沒有達到理論峰值，嘗試給出原因。

方法

CUDA矩陣的優化有多個思路，在本次試驗中我使用了shared memory進行訪問速度的提升，嘗試減少if-else語句的出現，避免序列化，同時做了精度優化以降低錯誤率（結果不怎麼好）。
同時，參考Nvidia給的Samples中0_simple裡的matrixMulCUBLAS相關程式碼，思考提升空間。

實驗

結果及分析

1.假設矩陣維度為n

處理速度公式=2*n/1000000000/time;
頻寬計算公式:
= ( sizeof(int)*dim + sizeof(int)*n + sizeof(float)*n
+ sizeof(float)*dim*2)/1000000000/time;

系統理論峰值(即浮點數理論峰值)
叢集理論浮點峰值
= CPU主頻（GHz）× CPU每時鐘週期執行浮點運算次數 × 節點數 × 8（每節點雙路四核）
=4.2*4*8=134.4GFLOPS/s

峰值頻寬： B=F×D/8=2133MHz*64bit/8=17.064GHz

沒有達到理論峰值的原因是:
程式並不只是在做浮點數運算或只是在訪問記憶體；
sgemm中還存在著if-else語句，使得執行緒存在著divergence；
由於大小分配的問題存在著Occupancy；
存在著空閒的執行緒；
以及作業系統的執行緒排程，和伺服器本身的不穩定性等等。

2.優化過程

2.1嘗試shared memory

Shared memory的作用在於降低對於全域性資料的訪問，充分利用Cuda中執行緒可以有獨立的記憶體空間及暫存器，以及block中執行緒之間可以通訊的特點
在shared memory大小定義中，Width要保證不能大於XY對應dim的最小值，另外在測試的時候發現，如果width_size大於32，那麼得到的結果是全錯（無論XY的dim有多大）暫時不清楚為什麼。

2.2嘗試減少if-else語句

在Sgemm函式中，if-else語句主要用於進行邊界判斷。
這是因為在分配block大小的時候，矩陣的維度可能不能剛好被32整除。例如dim=500時，不進行邊界判斷會引起很多問題。
一個有效的解決方案是，利用ceil的取整函式，在for迴圈中有效限制i的上界。使得對矩陣維度的限制沒有那麼大。

在程式碼中對grid, block 定義如下
dim3 block(DIM_THREAD_BLOCK_Y, DIM_THREAD_BLOCK_Y);
  dim3 grid((size_t)ceil( ((float)N) / ((float)block.x) ), 
            (size_t)ceil( ((float)N) / ((float)block.y)) );
  //取整函式ceil
複製程式碼

當然，經過反覆測試表明，矩陣的維度若能被32整除，其效能表現要比不能整除的要好。
另外在搜尋查詢的時候看到有一個方式是利用了cudaMallocPitch（），在分配的時候動態設定邊界大小，但是參考呼叫之後其優化的效果不是很明顯，沒有原作者所說的三倍效能提升，可能和本人的相關知識掌握不足有關。