使用cublasSgemm实现行主序矩阵乘法

7

我正在尝试使用cublasSgemm来计算两个非方阵以行主序存储的矩阵的乘积。我知道这个函数有一个参数可以指定是否要转置矩阵(CUBLAS_OP_T),但结果是以列主序存储的,而我需要以行主序存储。

此外,我的代码无法使用参数CUBLAS_OP_T来计算非方阵的乘积,只能计算方阵或非方阵的乘积且不进行转置。

另外,我知道声明矩阵为列主序也是一种选择。

define IDX2C(i,j,ld) (((j)*(ld))+(i)) 

但这不是一种选择,因为我必须使用的矩阵将在其他程序中设置。

我想互联网上有很多信息,但我找不到任何关于我的问题的答案。

我的代码如下:


    int m = 2;
    int k = 3;
    int n = 4;
    int print = 1;
    cudaError_t cudaStat; // cudaMalloc status
    cublasStatus_t stat; // CUBLAS functions status
    cublasHandle_t handle; // CUBLAS context

    int i,j;

    float *a, *b,*c;

    //malloc for a,b,c...

    // define a mxk matrix a row by row
    int ind =11;                               
    for(j=0;j<m*k;j++){                     
        a[j]=(float)ind++;                                    
    }                                        

    // define a kxn matrix b column by column
    ind =11;                               
    for(j=0;j<k*n;j++){                     
        b[j]=(float)ind++;                                    
    }                                        


    // DEVICE
    float *d_a, *d_b, *d_c; 

    //cudaMalloc for d_a, d_b, d_c...

    stat = cublasCreate(&handle); // initialize CUBLAS context

    stat = cublasSetMatrix(m,k, sizeof(*a), a, m, d_a, m);
    stat = cublasSetMatrix(k,n, sizeof(*b), b, k, d_b, k); 

    float al =1.0f;
    float bet =0.0f; 

    stat=cublasSgemm(handle,CUBLAS_OP_T,CUBLAS_OP_T,m,n,k,&al,d_a,m,d_b,k,&bet,d_c,m);

    stat = cublasGetMatrix (m,n, sizeof (*c) ,d_c ,m,c,m); // cp d_c - >c

    if(print == 1) {
    printf ("\nc after Sgemm :\n");
        for(i=0;i<m*n;i++){
                printf ("%f ",c[i]); // print c after Sgemm
        }
    }

    cudaFree (d_a); 
    cudaFree (d_b); 
    cudaFree (d_c);
    cublasDestroy (handle); // destroy CUBLAS context
    free (a); 
    free (b); 
    free (c); 

    return EXIT_SUCCESS ;
}

输出结果是矩阵 A * B 的转置,即:(A * B)T。

但我想要的是以行优先顺序计算的 C = A * B。

我希望有人可以帮助我。


2
这个可能会有兴趣:https://peterwittek.com/cublas-matrix-c-style.html - Robert Crovella
1个回答

15
如您所述,cuBLAS将矩阵解释为列优先顺序,因此当您执行cublasSgemm(handle,CUBLAS_OP_T,CUBLAS_OP_T,m,n,k,&al,d_a,m,d_b,k,&bet,d_c,m)时,您正确地对每个输入进行了转置(这些输入是以行优先形式创建的),以准备按列优先解释。问题在于,cuBLAS还以列主要顺序转储结果。
我们将欺骗cuBLAS计算,它将以列主要顺序输出,因此在以行优先顺序巧妙解释时,它将看起来像。因此,我们不是计算AB = C,而是计算 = 。幸运的是,已由创建A和B时采用行主要顺序的动作获得,因此我们可以通过CUBLAS_OP_N简单地绕过转置。所以将该行更改为cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,n,m,k,&al,d_b,n,d_a,k,&bet,d_c,n)
您提供的示例代码应计算。

经过更新后调用 cublasSgemm,我们得到了以上结果:

c after Sgemm :
548.000000 584.000000 620.000000 656.000000 683.000000 728.000000 773.000000 818.000000 

1
当BLAS实现允许行/列主要说明符时,这是大多数BLAS实现在幕后执行的操作吗? - hillard28

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接