自己写的cuda函数和cublas和ispc的对比（均支持非方阵的计算）

Posted 帅的发光发亮

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了自己写的cuda函数和cublas和ispc的对比（均支持非方阵的计算）相关的知识，希望对你有一定的参考价值。

自己写的cuda函数和cublas和ispc的对比

！！！代码在文末 ！！！

开文废话

憋了这么久，终于开始写了。

预备知识

因为cublas的数据存储是按照列优先的，而c/c++是按行存储的。

行优先还是列优先

首先了解“行优先”和“列优先”的知识，这两种方式在数学上的直观描述如下，给定如下矩阵：

矩阵在逻辑上表达为2维的矩阵，M行K列，但存储到内存的时候都是按一维布局，其中按行优先存储和按列优先存储的差异如上图所示

如上图所示，当矩阵按行优先存储然后又按相反的列优先读取的话，就会得到原矩阵转置的结果；同理适用于按列优先存储然后按行优先读取。

例 cublasSgemm 函数

cublasStatus_t cublasSgemm(cublasHandle_t handle,
                           cublasOperation_t transa, cublasOperation_t transb,
                           int m, int n, int k,
                           const float           *alpha,
                           const float           *A, int lda,
                           const float           *B, int ldb,
                           const float           *beta,
                           float           *C, int ldc)


cublasSgemm的官方API说明文档链接cublas

  • 根据文档说可以知道，cublasSgemm完成了 C = alpha * op ( A ) * op ( B ) + beta * C 的矩阵乘、加运算。

  • 其中alpha和beta是标量， A、 B、 C是以列优先存储的矩阵，A称为乘法左矩阵、B称为乘法右矩阵、C称为结果矩阵，当alpha = 1.0f 并且 beta =0.0f 的时候 cublasSgemm完成了计算： 结果矩阵= op (乘法左矩阵) * op ( 乘法右矩阵) 。

  • cublasOperation_t 该类型表明输入的密集矩阵的形式，其值有 CUBLAS_OP_N（非转置）;CUBLAS_OP_T（转置）; CUBLAS_OP_C（共轭转置）。该函数对应于BLAS（FORTRAN版）的变量字符’N’或’n’(非转置,即正常形式的矩阵)，T’或’t’(转置矩阵);‘C’或’c’(共轭转置矩阵，对应的是复数矩阵。

求解C=AxB

其中(A为A_ROW行A_COL列 B为B_ROW行B_COL列 所以 C为A_ROW行B_COL列）

不使用cublasSgemm transa与transb参数

由于C/C++程序中输入的A和B是按行存储，所以在的情况下，cublas其实读取到的是A和B的转置矩阵AT和BT.

根据线性代数的规则可知CT = (A x B)T = BT x AT 所以cublasSgemm API中几个参数设置如下：

  • 设置了cublasSgemm的transa与transb参数=CUBLAS_OP_N
  • 乘法左矩阵为BT=参数设置为B，乘法右矩阵为AT=参数设置为A
  • 结果矩阵的行数为CT的行数（正常c/c++矩阵的列数）= 参数设置为:B_COL
  • 结果矩阵的列数为CT的列数（正常c/c++矩阵的行数）= 参数设置为:A_ROW
  • 乘法左矩阵列与乘法右矩阵的行=参数设置为:B_ROW
  • 按列优先读取乘法左矩阵B的主维度（即BT有几行、正常c/c++矩阵的列数）=参数设置为B_COL
  • 按列优先读取乘法右矩阵A的主维度（即AT有几行、正常c/c++矩阵的列数）=参数设置为A_COL
  • 结果矩阵存储在参数C中，它的主维度（即结果矩阵c有几行）= 参数设置为B_COL

cublasSgemm(
handle,
CUBLAS_OP_N, //矩阵A的属性参数，不转置，按列优先
CUBLAS_OP_N, //矩阵B的属性参数，不转置，按列优先
B_COL, //矩阵BT、CT的行数
A_ROW, //矩阵AT、CT的列数
B_ROW, //BT的列数，AT的行数，此处也可为A_COL,一样的
&a, //alpha的值
d_B, //左矩阵，为B^T
B_COL, //B^T的leading dimension，按列优先，则leading dimension为B^T的行数(B的列数)
d_A, //右矩阵，为A^T
A_COL, //A^T的leading dimension，按列优先，则leading dimension为A^T的行数(A的列数)
&b, //beta的值
d_C, //结果矩阵C
B_COL //C^T的leading dimension，C^T矩阵一定按列优先，则leading dimension为C^T的行数(C的列数)
);//此处的h_C是按列存储的CT

按上面的参数调用cublasSgemm API (矩阵A按行存储在指针d_A, 矩阵B按行存储在指针d_B, 矩阵C的存储空间指针d_C) 最后结果d_C传输到host端的h_C,从结果矩阵的存储空间h_C中按行读取到的就是C=AxB的结果，整个cublasSgemm的计算过程如下图所示。

部分示例代码：

template <typename T1, typename T2>
void cublas_kernel()

    // 定义状态变量
    cublasHandle_t handle;
    cublasCreate(&handle);
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    //存储于内存中的矩阵
    T1* h_A, * h_B;
    T2* h_C, * h_CC;

    //在内存中开辟空间
    cudaHostAlloc((void**)&h_A, sizeof(T1) * A_ROW * A_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_B, sizeof(T1) * B_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_C, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_CC, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);

    MatrixINIT<T1>(A_ROW, A_COL, h_A);
    MatrixINIT<T1>(B_ROW, B_COL, h_B);

    // 打印待测试的矩阵
#if defined(USE_FLOAT_T)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);

#elif defined(USE_FLOAT_N)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);

#elif defined(USE_DOUBLE_T)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);

#elif defined(USE_DOUBLE_N)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);
#endif

    //存储于显存中的矩阵
    T1* d_A, * d_B;
    T2* d_C;

    cudaMalloc((void**)&d_A, sizeof(T1) * A_ROW * A_COL);
    cudaMalloc((void**)&d_B, sizeof(T1) * B_ROW * B_COL);
    cudaMalloc((void**)&d_C, sizeof(T2) * A_ROW * B_COL);
  
    /*
    cudaHostAlloc((void**)&d_A, sizeof(T1) * A_ROW * A_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&d_B, sizeof(T1) * B_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&d_C, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);
    */


    //创建流对象，用于任务级(Grid)同步
    cudaStream_t stream;
    cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);

    cublasSetStream(handle, stream);

    const T2 a = 1.0f, b = 0.0f;

    //计时开始
    TIMER_START(_X);
    //cudaEventRecord(start, 0);

    for (int i = 0; i < N  ; i++)
    
        //数据从Host端拷贝到Device端、 cubals方式
        /*
        cublasSetMatrix(A_ROW, A_COL, sizeof(*h_A), h_A, A_ROW, d_A, A_ROW);
        cublasSetMatrix(B_ROW, B_COL, sizeof(*h_B), h_B, B_ROW, d_B, B_ROW);
        */

        cublasSetMatrixAsync(A_ROW, A_COL, sizeof(*h_A), h_A, A_ROW, d_A, A_ROW, stream);
        cublasSetMatrixAsync(B_ROW, B_COL, sizeof(*h_B), h_B, B_ROW, d_B, B_ROW, stream);

        //数据从Host端拷贝到Device端、 传统方式
        //cudaMemcpy(d_A, H_A, sizeof(T1) * A_ROW * A_COL, cudaMemcpyHostToDevice); 
        //cudaMemcpy(d_B, H_B, sizeof(T1) * B_ROW * B_COL, cudaMemcpyHostToDevice);

        //单独计算核函数运算时间
        cudaEventRecord(start, 0);

#if defined(USE_FLOAT_N)
        cublasSgemm(
            handle,
            CUBLAS_OP_N,   //矩阵A的属性参数，不转置，按列优先
            CUBLAS_OP_N,   //矩阵B的属性参数，不转置，按列优先
            B_COL,          //矩阵B^T、C^T的行数
            A_ROW,          //矩阵A^T、C^T的列数
            B_ROW,          //B^T的列数，A^T的行数，此处也可为A_COL,一样的
            &a,             //alpha的值
            d_B,            //左矩阵，为B^T
            B_COL,          //B^T的leading dimension，按列优先，则leading dimension为B^T的行数(B的列数)
            d_A,            //右矩阵，为A^T
            A_COL,          //A^T的leading dimension，按列优先，则leading dimension为A^T的行数(A的列数)
            &b,             //beta的值
            d_C,            //结果矩阵C
            B_COL           //C^T的leading dimension，C^T矩阵一定按列优先，则leading dimension为C^T的行数(C的列数)
        );//此处的h_C是按列存储的C^T


#elif defined(USE_DOUBLE_N)
        cublasDgemm(
            handle,
            CUBLAS_OP_N,   //矩阵A的属性参数，不转置，按列优先
            CUBLAS_OP_N,   //矩阵B的属性参数，不转置，按列优先
            B_COL,          //矩阵B^T、C^T的行数
            A_ROW,          //矩阵A^T、C^T的列数
            B_ROW,          //B^T的列数，A^T的行数，此处也可为A_COL,一样的
            &a,             //alpha的值
            d_B,            //左矩阵，为B^T
            B_COL,          //B^T的leading dimension，按列优先，则leading dimension为B^T的行数(B的列数)
            d_A,            //右矩阵，为A^T
            A_COL,          //A^T的leading dimension，按列优先，则leading dimension为A^T的行数(A的列数)
            &b,             //beta的值
            d_C,            //结果矩阵C
            B_COL           //C^T的leading dimension，C^T矩阵一定按列优先，则leading dimension为C^T的行数(C的列数)
        );//此处的h_C是按列存储的C^T

#endif

        //计时结束
        cudaDeviceSynchronize();
        cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);
        cudaEventElapsedTime(&elapsedTime, start, stop);

        //TIMER_STOP(_X);
        //cout << "GPU耗费了: " << TIMER_MSEC(_X) << " ms " << "\\n";
        // 
        //将Device端计算完的结果传输会Host端  cublas方式
        //cublasGetMatrix(A_ROW, B_COL, sizeof(*h_C), d_C, A_ROW, h_C, A_ROW);
        cublasGetMatrixAsync(A_ROW, B_COL, sizeof(*h_C), d_C, A_ROW, h_C, A_ROW, stream);
        //传统方式
        //cudaMemcpy(H_C, d_C, sizeof(T2) * A_ROW * B_COL, cudaMemcpyDeviceToHost);
    
    TIMER_STOP(_X);
    /*
    cudaDeviceSynchronize();
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime_cublas, start, stop);
    */
    //打印结果
    cout << "cublas_kernel GPU传输、计算花费了:  " << TIMER_MSEC(_X) << " ms " << "\\n";
    //std::cout<< "GPU传输、计算花费了：" << elapsedTime_cublas << " ms" << std::endl;
    std::cout << "cublas_kernel GPU计算花费了：" << elapsedTime * N<< " ms" << std::endl<< std::endl;

#if defined(USE_FLOAT_N)
    //按行读取h_C相当于做了CTT=C的结果
    Matrixshow<T2>("计算结果C的值 ( C^T = (B^T*A^T) = (B*A)^T )", A_ROW, B_COL, h_C, 1, 0);
    cout << endl;

#if defined(USE_CPU_COST)
    cpu_matrix_mult<T1, T2>(h_A, h_B, A_ROW, A_COL, B_COL, h_C, h_CC, 0);
#endif 

#elif defined(USE_DOUBLE_N)
    //按行读取h_C相当于做了CTT=C的结果
    Matrixshow<T2>("计算结果C的值 ( C^T = (B^T*A^T) = (B*A)^T )", A_ROW, B_COL, h_C, 1, 0);
    cout << endl;

#if defined(USE_CPU_COST)
    cpu_matrix_mult<T1, T2>(h_A, h_B, A_ROW, A_COL, B_COL, h_C, h_CC, 0);
#endif 

#endif

    //释放内存
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    /*
    cudaFreeHost(d_A);
    cudaFreeHost(d_B);
    cudaFreeHost(d_C);
    */
    cudaFreeHost(h_A);
    cudaFreeHost(h_B);
    cudaFreeHost(h_C);
    cudaFreeHost(h_CC);
    cublasDestroy(handle);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    cudaStreamDestroy(stream);

使用cublasSgemm transa与transb参数

由于C/C++程序中输入的A和B是按行存储，所以在的情况下，cublas其实读取到的是A和B的转置矩阵A^T
和 B^T
设置了cublasSgemm的transa与transb参数后，可以在做矩阵运算前对读取到的A^T、
B^T矩阵做一次转置，获得A和B根据线性代数的规则可知C = A x B 所以cublasSgemm API中几个参数设置如下:

  • 设置了cublasSgemm的transa与transb参数=CUBLAS_OP_T，在进行矩阵运算前对读取的矩阵做一次转置
  • 乘法左矩阵为A=参数设置为A，乘法右矩阵为B=参数设置为B
  • 结果矩阵的行数为C的行数=参数设置为 A_ROW
  • 结果矩阵的列数为C的列数=参数设置为 B_COL
  • 乘法左矩阵列与乘法右矩阵的行=参数设置为 A_COL
  • 按列优先读取乘法左矩阵A的主维度(就变成了类似CUBLAS_OP_N的参数情况)（即A^T有几行）=参数设置为 A_COL
  • 按列优先读取乘法右矩阵B的主维度(就变成了类似CUBLAS_OP_N的参数情况)（即B^T有几行）=参数设置为 B_COL
  • 结果矩阵存储在参数C中，它的主维度（即结果矩阵c有几行）= 参数设置为 A_ROW

cublasSgemm(
handle,
CUBLAS_OP_T, //矩阵A的属性参数，还是按列优先读取，但是在计算前，转置，变成正常c/c++的方式
CUBLAS_OP_T, //矩阵B的属性参数，还是按列优先读取，但是在计算前，转置，变成正常c/c++的方式
A_ROW, //矩阵A、C的行数
B_COL, //矩阵B、C的列数
A_COL, //A的列数，B的行数，此处也可为B_ROW一样的
&a, //alpha的值
d_A, //左矩阵，为A
A_COL, //A的leading dimension，按列优先，则leading dimension为（A^T的行数）A的列数
d_B, //右矩阵，为B
B_COL, //B的leading dimension，按列优先，则leading dimension为（B^T的行数）A的列数
&b, //beta的值
d_C, //结果矩阵C
A_ROW //C的leading dimension，C矩阵一定按列优先，则leading dimension为C的行数
);//此处的h_C是按列存储的C

红色的参数标记出与“不使用cublasSgemm transa与transb参数”例子中的不同，按上面的参数调用cublasSgemm API (矩阵A按行存储在指针d_a, 矩阵B按行存储在指针d_b, 矩阵C的存储空间指针d_c) 最后从结果矩阵的存储空间d_c中按行读取到的就是C=AxB后C^T
的结果，所以在C/C++程序中还需要对读取的结果C^T做一次矩阵转置操作才能获得最终正确的C。整个cublasSgemm的计算过程如下图所示:

结果：

部分示例代码


template <typename T1, typename T2>
void cublas_kernel()

    // 定义状态变量
    cublasHandle_t handle;
    cublasCreate(&handle);
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    //存储于内存中的矩阵
    T1* h_A, * h_B;
    T2* h_C, * h_CC;

    //在内存中开辟空间
    cudaHostAlloc((void**)&h_A, sizeof(T1) * A_ROW * A_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_B, sizeof(T1) * B_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_C, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_CC, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);

    MatrixINIT<T1>(A_ROW, A_COL, h_A);
    MatrixINIT<T1>(B_ROW, B_COL, h_B);

    // 打印待测试的矩阵
#if defined(USE_FLOAT_T)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 1);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 1);

#elif defined(USE_FLOAT_N)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);

#elif defined(USE_DOUBLE_T)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);

#elif defined(USE_DOUBLE_N)
    Matrixshow<T1>("A", A_ROW, A_COL, h_A, 0);
    Matrixshow<T1>("B", B_ROW, B_COL, h_B, 0);
#endif

    //存储于显存中的矩阵
    T1* d_A, * d_B;
    T2* d_C;

    cudaMalloc((void**)&d_A, sizeof(T1) * A_ROW * A_COL);
    cudaMalloc((void**)&d_B, sizeof(T1) * B_ROW * B_COL);
    cudaMalloc((void**)&d_C, sizeof(T2) * A_ROW * B_COL);
  
    /*
    cudaHostAlloc((void**)&d_A, sizeof(T1) * A_ROW * A_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&d_B, sizeof(T1) * B_ROW * B_COL, cudaHostAllocDefault);
    cudaHostAlloc((void**)&d_C, sizeof(T2) * A_ROW * B_COL, cudaHostAllocDefault);
    */


    //创建流对象，用于任务级(Grid)同步
    cudaStream_t stream;
    cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);

    cublasSetStream(handle, stream);

    const T2 a = 1.0f, b = 0.0f;

    //计时开始
    TIMER_START(_X);
    //cudaEventRecord(start, 0);

    for (int i = 0; i < N  ; i++)
    
        //数据从Host端拷贝到Device端、 cubals方式
        /*
        cublasSetMatrix(A_ROW, A_COL, sizeof(*h_A), h_A, A_ROW, d_A, A_ROW);
        cublasSetMatrix(B_ROW, B_COL, sizeof(*h_B), h_B, B_ROW, d_B, B_ROW);
        */

        cublasSetMatrixAsync(A_ROW, A_COL, sizeof(*h_A), h_A, A_ROW, d_A, A_ROW, stream);
        cublasSetMatrixAsync(B_ROW, B_COL, sizeof(*h_B), h_B, B_ROW, d_B, B_ROW, stream);

        //数据从Host端拷贝到Device端、 传统方式
        //cudaMemcpy(d_A, H_A, sizeof(T1) * A_ROW * A_COL, cudaMemcpyHostToDevice); 
        //cudaMemcpy(d_B, H_B, sizeof(T1) * B_ROW * B_COL, cudaMemcpyHostToDevice);

        //单独计算核函数运算时间
        cudaEventRecord(start, 0);


#if defined(USE_FLOAT_T)
        cublasSgemm(
            handle,
            CUBLAS_OP_T,   //矩阵A的属性
以上是关于自己写的cuda函数和cublas和ispc的对比（均支持非方阵的计算）的主要内容，如果未能解决你的问题，请参考以下文章
 

 自己写的cuda函数和cublas和ispc的对比（均支持非方阵的计算）


 自己写的cuda函数和cublas和ispc的对比（均支持非方阵的计算）


 CUDA 内核可以调用 cublas 函数吗？


 加速图像处理的神器: INTEL ISPC编译器 迁移图像旋转算法 - 从C代码双精度到 ISPC双精度


 是否可以从全局或设备函数调用 CUDA CUBLAS 函数


 cublas 函数调用 cublasSgemv