cuda 矩陣乘法函數之cublasSgemm

，可以考慮使用, 例如cublasSgeam()(矩陣加法), 進行一次1.0 * AT + 0.0 * B的參數設定, 利用内置的轉置功能(注意這裡的1和0), 來進行将A轉換成AT. 在使用CUDA的cuBLAS庫中矩陣乘法函數cublasSgemm時，注意到cuda其中的二維矩陣的儲存是“按列儲存”，一天都處于蒙蔽狀态，查了很多資料，按所得結果情況，總結出如下幾條。

一、獲得按行存儲的結果

由博文：http://blog.csdn.net/xfortius/article/details/9225799收到啟發：

比如，我們想求C=A*B這個矩陣運算，其中A={{1,1},{2,2},{3,3}};B={{1},{1}};C={{4},{5},{6}}，而對于A、B、C進行一維數組表示有A={1,1,2,2,3,3}，B={1,1}，C={4,5,6}；這個在c/c++是和前面表示一樣，但是在cublasSgemm中就完全不對了，那麼這個一維的A其實表示的是{{1,2,3},{1,2,3}};可以看到兩個矩陣其實剛好是轉置關系。那麼我們要求C=A*B,按照一維資料輸入的話結果A表示的是AT，B表示的是BT，是以我們要輸入的是AT和BT，這樣在公式中得到的才是A*B.假設這是得到矩陣C=A*B，但是C也是按列存儲的，我們要的是CT,而CT=BT*AT，而這裡的BT其實就是原矩陣B，AT其實就是原矩陣A。可見，我們通過交換AB的順序就可以得到按行存儲的C。

這裡還有一點，就是cublasSgemm的參數，把自己繞的有點暈。這裡我們輸入的是B*A即（2*1）*（3*2），但是真正在函數中執行的是BT*AT=CT(1*2)*(2*3=(1*3)).是以主要不要把參數搞錯了。。

#include <cublas_v2.h> //cuda自帶庫函數  
#include <helper_cuda.h>  
#include <stdio.h>  
int main(void)  
{  
    float alpha=;  
    float beta=;  
    float h_A[]={,,,,,};  
    float h_B[]={,};  
    float h_C[];  
    float *d_a,*d_b,*d_c;  
    checkCudaErrors(cudaMalloc((void**)&d_a,*sizeof(float)));  
    checkCudaErrors(cudaMalloc((void**)&d_b,*sizeof(float)));  
    checkCudaErrors(cudaMalloc((void**)&d_c,*sizeof(float)));  
    checkCudaErrors(cudaMemcpy(d_a,&h_A,*sizeof(float),cudaMemcpyHostToDevice));  
    checkCudaErrors(cudaMemcpy(d_b,&h_B,*sizeof(float),cudaMemcpyHostToDevice));  
    checkCudaErrors(cudaMemset(d_c,,*sizeof(float)));  
    cublasHandle_t handle;  
    cublasCreate(&handle);  
    cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,,,,&alpha,d_b,,d_a,,&beta,d_c,);  
    checkCudaErrors(cudaMemcpy(h_C,d_c,*sizeof(float),cudaMemcpyDeviceToHost));  
    for(int i=;i<;i++)  
    {  
        printf("%f\n",h_C[i]);  
    }  
    printf("\n");  
    return ;  
}  
           

總結：若想獲得按行存儲的結果矩陣C=A*B（可複制到host）,可以使用代碼：

cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,m,n,k,B,ldb,A,lda,C,ldc);
           

說明：

1. 輸入的參數為CUBLAS_OP_N ， CUBLAS_OP_N
2. 相比原參數，A和B的位置調換（C=A*B —>CT=(BT*AT) (其中CT按列存儲，實際上就是按行存儲的C ！)）
3. 其中參數m為col of B ,n 為row of A , k為A和B相同的(row of B or col of A) ,參數ldb , lda , ldc 分别為col of B,A,C ;

經過驗證，結果正确。

二、獲得按列存儲結果（存儲在GPU顯存中）

有些矩陣為中間結果，為了避免重複使用上述方法，可以使用下面的方法，其中，GPU顯存中的為按列存儲，Host記憶體中的為按行存儲。
           

靈感來源博文：http://www.xuebuyuan.com/1966642.html

由于cublas為了更大的适應Fortan語言，二維資料的存儲采用以列優先的方式，這與C/C++中，行優先的存儲方式不同。由于本人的研究是資料的來源是C代碼得到的，為了加速矩陣的運算效率，利用cublas來完成。本文檔提出了一種有效的解決方案。

為了更好的說明，以函數cublasSgemm的實作C= A*B為例。接口cublasSgemm 實作的功能為C = alpha*A*B + beta*C，為了完成C= A*B 的功能，令alpha= 1.0f，beta = 0.0f
           

利用cublasSgemm的參數 transa（transb） 和 lda（ldb）的設定來共同解決存儲方式改變的問題。

• 第一種：輸入矩陣A,B均在CPU上

  cublasSgemm(‘t’,’t’,結果矩陣C的行,結果矩陣C的列，左矩陣A的列, alpha，左矩陣,左矩陣的列,右矩陣,右矩陣的列, beta,結果矩陣,結果矩陣的行);

  如果求AT*B，隻需将A的t改為n（t的轉置為t）,并且關于A的參數（行和列）是轉置以前的。

• 第二種：輸入矩陣A在顯存上以column-major的方式存儲，B在CPU上

  cublasSgemm(‘n’,’t’,結果矩陣C的行,結果矩陣C的列，左矩陣A的列, alpha，左矩陣A,左矩陣A的行,右矩陣B,右矩陣B的列, beta,結果矩陣C,結果矩陣C的行);

• 第三種：輸入矩陣B在顯存上以column-major的方式存儲，A在CPU上

  cublasSgemm(‘t’,’n’,結果矩陣C的行,結果矩陣C的列，左矩陣A的列, alpha，左矩陣A,左矩陣A的列,右矩陣B,右矩陣B的行, beta,結果矩陣C,結果矩陣C的行);

• 第四種：輸入矩陣A，B均在顯存上，以column-major的方式存儲

  cublasSgemm(‘n’,’n’,結果矩陣C的行,結果矩陣C的列，左矩陣A的列, alpha，左矩陣A,左矩陣A的行,右矩陣B,右矩陣B的行, beta,結果矩陣C,結果矩陣C的行);

得到的結果C都是按列存儲的（拷貝出host記憶體之前轉置一下，可以考慮使用, 例如cublasSgeam()(矩陣加法), 進行一次1.0 * AT + 0.0 * B的參數設定, 利用内置的轉置功能(注意這裡的1和0), 來進行将A轉換成AT.）。

經過驗證，結果正确。

總結：

如果前邊的參數是’t’，那麼leading dimesion 就是矩陣的列數，因為此時的矩陣是按照C語言以行優先的方式來存儲的；反之如果前邊的參數是’n’，那麼leading dimesion 就是矩陣的行數，此時的矩陣保持CUBLAS的列優先存儲方式。