自己写的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的对比（均支持非方阵的计算）