cublasSgemm详解

来源：互联网发布：下载种子最快的软件编辑：程序博客网时间：2024/05/21 07:12

cublasSgemm函数详解

cublasSgemm是CUDA的cublas库的矩阵相乘函数，由于cublas中矩阵的存储是列优先，所以cublasSgemm函数的参数容易让人误解，防止忘记，留笔记如下。

首先，在一般的C/C++程序中，我们开辟一段连续的内存，放入1,2,3,4,5,6,7,8,9，指定矩阵行和列均为3，则可表示矩阵[1,2,3 ; 4,5,6 ; 7,8,9]，然而，在使用cublas时，这样表示出来的矩阵应该是[1,4,7; 2,5,8; 3,6,9]。cublas的cublasSgemm函数完成C=αop(A)op(B)+βC的计算，当需要计算C=AB时，显然可以直接设置α=1,β=0。其中op操作对决定矩阵是否转置，即决定该矩阵是按照行优先还是列优先。当我们选择CUBLAS_OP_N时表示不转置，按列优先存储；当我们选择CUBLAS_OP_T时表示需要转置，按行优先存储。
这里写图片描述
直接看例子吧，注释中记录参数意义。
看到一个总结很不错，引用一下,：
如果前边的参数是’CUBLAS_OP_T’，那么leading dimesion 就是矩阵的列数，因为此时的矩阵是按照C语言以行优先的方式来存储的；反之如果前边的参数是’CUBLAS_OP_N’，那么leading dimesion 就是矩阵的行数，此时的矩阵保持CUBLAS的列优先存储方式。

按转置方式进行求解C=AB

// CUDA runtime 库 + CUBLAS 库#include "cuda_runtime.h"#include "cublas_v2.h"#include <iostream>#include <stdlib.h>using namespace std;// 定义测试矩阵的维度int const A_ROW = 5;int const A_COL = 6;int const B_ROW = 6;int const B_COL = 7;int main(){  // 定义状态变量  cublasStatus_t status;  float *h_A,*h_B,*h_C;   //存储于内存中的矩阵  h_A = (float*)malloc(sizeof(float)*A_ROW*A_COL);  //在内存中开辟空间  h_B = (float*)malloc(sizeof(float)*B_ROW*B_COL);  h_C = (float*)malloc(sizeof(float)*A_ROW*B_COL);  // 为待运算矩阵的元素赋予 0-10 范围内的随机数  for (int i=0; i<A_ROW*A_COL; i++) {    h_A[i] = (float)(rand()%10+1);  }  for(int i=0;i<B_ROW*B_COL; i++) {    h_B[i] = (float)(rand()%10+1);  }  // 打印待测试的矩阵  cout << "矩阵 A :" << endl;  for (int i=0; i<A_ROW*A_COL; i++){    cout << h_A[i] << " ";    if ((i+1)%A_COL == 0) cout << endl;  }  cout << endl;  cout << "矩阵 B :" << endl;  for (int i=0; i<B_ROW*B_COL; i++){    cout << h_B[i] << " ";    if ((i+1)%B_COL == 0) cout << endl;  }  cout << endl;  float *d_A,*d_B,*d_C;    //存储于显存中的矩阵  cudaMalloc((void**)&d_A,sizeof(float)*A_ROW*A_COL); //在显存中开辟空间  cudaMalloc((void**)&d_B,sizeof(float)*B_ROW*B_COL);  cudaMalloc((void**)&d_C,sizeof(float)*A_ROW*B_COL);  cublasHandle_t handle;  cublasCreate(&handle);  cudaMemcpy(d_A,h_A,sizeof(float)*A_ROW*A_COL,cudaMemcpyHostToDevice); //数据从内存拷贝到显存  cudaMemcpy(d_B,h_B,sizeof(float)*B_ROW*B_COL,cudaMemcpyHostToDevice);  float a = 1, b = 0;  cublasSgemm(          handle,          CUBLAS_OP_T,   //矩阵A的属性参数，转置，按行优先          CUBLAS_OP_T,   //矩阵B的属性参数，转置，按行优先          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的列数          d_B,            //右矩阵，为B          B_COL,          //B的leading dimension，此时选择转置，按行优先，则leading dimension为B的列数          &b,             //beta的值          d_C,            //结果矩阵C          A_ROW           //C的leading dimension，C矩阵一定按列优先，则leading dimension为C的行数  );  //此时得到的结果便是C=AB,但由于C是按列优先，故此时得到的C应该是正确结果的转置  std::cout << "计算结果的转置 ( (A*B)的转置 )：" << std::endl;  cudaMemcpy(h_C,d_C,sizeof(float)*A_ROW*B_COL,cudaMemcpyDeviceToHost);  for(int i=0;i<A_ROW*B_COL;++i) {    std::cout<<h_C[i]<<" ";    if((i+1)%B_COL==0) std::cout<<std::endl;  }  cudaFree(d_A);  cudaFree(d_B);  cudaFree(d_C);  free(h_A);  free(h_B);  free(h_C);  return 0;}

不按转置方式直接求解

前一种方法求得的C是正确结果的转置，那么我们可以直接求CT，不就是正确结果了吗？CT=(AB)T=BTAT，此时参与运算的是BT和AT，那么我们都不需要像上一种方法一样选择CUBLAS_OP_T参数来对AB进行转置了，指定了op(A)和op(B)的维数，按列填充得到的就是BT和AT。所以此时直接选择CUBLAS_OP_N就好了。代码如下：

// CUDA runtime 库 + CUBLAS 库#include "cuda_runtime.h"#include "cublas_v2.h"#include <iostream>#include <stdlib.h>using namespace std;// 定义测试矩阵的维度int const A_ROW = 5;int const A_COL = 6;int const B_ROW = 6;int const B_COL = 7;int main(){  // 定义状态变量  cublasStatus_t status;  float *h_A,*h_B,*h_C;   //存储于内存中的矩阵  h_A = (float*)malloc(sizeof(float)*A_ROW*A_COL);  //在内存中开辟空间  h_B = (float*)malloc(sizeof(float)*B_ROW*B_COL);  h_C = (float*)malloc(sizeof(float)*A_ROW*B_COL);  // 为待运算矩阵的元素赋予 0-10 范围内的随机数  for (int i=0; i<A_ROW*A_COL; i++) {    h_A[i] = (float)(rand()%10+1);  }  for(int i=0;i<B_ROW*B_COL; i++) {    h_B[i] = (float)(rand()%10+1);  }  // 打印待测试的矩阵  cout << "矩阵 A :" << endl;  for (int i=0; i<A_ROW*A_COL; i++){    cout << h_A[i] << " ";    if ((i+1)%A_COL == 0) cout << endl;  }  cout << endl;  cout << "矩阵 B :" << endl;  for (int i=0; i<B_ROW*B_COL; i++){    cout << h_B[i] << " ";    if ((i+1)%B_COL == 0) cout << endl;  }  cout << endl;  float *d_A,*d_B,*d_C;    //存储于显存中的矩阵  cudaMalloc((void**)&d_A,sizeof(float)*A_ROW*A_COL); //在显存中开辟空间  cudaMalloc((void**)&d_B,sizeof(float)*B_ROW*B_COL);  cudaMalloc((void**)&d_C,sizeof(float)*A_ROW*B_COL);  cublasHandle_t handle;  cublasCreate(&handle);  cudaMemcpy(d_A,h_A,sizeof(float)*A_ROW*A_COL,cudaMemcpyHostToDevice); //数据从内存拷贝到显存  cudaMemcpy(d_B,h_B,sizeof(float)*B_ROW*B_COL,cudaMemcpyHostToDevice);  float a = 1, b = 0;  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的列数)  );  //此时得到的结果便是C=AB,但由于C是按列优先，故此时得到的C应该是正确结果的转置  std::cout << "计算结果的转置 ( (A*B)的转置 )：" << std::endl;  cudaMemcpy(h_C,d_C,sizeof(float)*A_ROW*B_COL,cudaMemcpyDeviceToHost);  for(int i=0;i<A_ROW*B_COL;++i) {    std::cout<<h_C[i]<<" ";    if((i+1)%B_COL==0) std::cout<<std::endl;  }  cudaFree(d_A);  cudaFree(d_B);  cudaFree(d_C);  free(h_A);  free(h_B);  free(h_C);  return 0;}

阅读全文

0 0