cuBLAS的使用

cuBLAS API主要特性

1. 错误码：所有cuBLAS库函数的都返回cublasStatus_t类型的错误码。
2. cuBLAS上下文：调用cublasCreate()函数初始化cuBLAS库的上下文，完成计算后调用cublasDestroy()释放上下文相关资源。
3. 线程安全：cuBLAS库是线程安全的，其函数可以应用在多线程场景。
4. 结果可复现：同一个cuBLAS版本，在GPU上使用相同架构和相同数量SM时，每次运行生成的结果应该都相同，不同版本的cuBLAS不保证结果相同。
5. 标量参数：只有两个类别的函数使用到了标量参数。获取α或者β引用作为比例因数的函数，如gemm；返回标量结果的函数，如amax()/amin()/asum()/rotg()/rotmg()/dot()/nrm2()。
6. Stream并行：如果应用程序支持多个独立的计算任务，可以使用CUDA的Stream技术来并行执行这些任务。
7. 批处理内核：可以使用Stream来批处理每个小内核函数，例如当一个应用程序需要计算许多小型独立的矩阵乘时的稠密矩阵，就可以使用批处理内核提升性能。
8. 缓存配置：缓存配置可以直接使用CUDA运行时函数cudaDeviceSetCacheConfig设置。cudaFuncSetCacheConfig函数是专门为一些函数设置缓存配置。
9. 支持静态库：cuBLAS库支持静态库的链接方式，其静态库为libculibos.a。
10. GEMM算法数值行为：一些GEMM算法会根据K维度来分割计算以提高GPU的使用率。对于cublas<t>gemmEx和cubalsGemmEx函数，当计算输入类型大于输出类型时，分割块的总和可能会有一些溢出导致最终的结果溢出。该表现可以使用cublasSerMathMode()函数计算精确模式CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION来规避。
11. 使用Tensor代码：从cuBLAS 11.0.0版本开始，cuBLAS库会自动地尽可能地使用Tensor特性，以加速矩阵乘法。
12. 支持CUDA Graphs：在大多数情况下，cuBLAS可以被CUDA Grapes stream使用。

cuBLAS库的使用方法

代码引用头文件"cublas.h"或者"cublas_v2.h"，编译时链接cuBLAS的动态库（Linux系统下动态cublas.so）。

示例代码：

//cublas_example.c, Application Using C and cuBLAS: 0-based indexing
//-----------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#include "cublas_v2.h" // 包含cublas头文件
#define M 6
#define N 5
#define IDX2C(i,j,ld) (((j)*(ld))+(i)) // 定义0-based数组的访问方式
static __inline__ void modify (cublasHandle_t handle, float *m, int ldm, int n, int p, int q, float alpha, float beta){
     cublasSscal (handle, n-q, &alpha, &m[IDX2C(p,q,ldm)], ldm);
     cublasSscal (handle, ldm-p, &beta, &m[IDX2C(p,q,ldm)], 1); 
}

int main (void){
     cudaError_t cudaStat;
     cublasStatus_t stat;
     cublasHandle_t handle;
     int i, j;
     float* devPtrA;
     float* a = 0;
     a = (float *)malloc (M * N * sizeof (*a)); // 申请CPU数组内存
     if (!a) {
         printf ("host memory allocation failed");
         return EXIT_FAILURE;
     }
     for (j = 0; j < N; j++) {
         for (i = 0; i < M; i++) {
             a[IDX2C(i,j,M)] = (float)(i * N + j + 1);
         }
     }
     cudaStat = cudaMalloc ((void**)&devPtrA, M*N*sizeof(*a)); // 申请GPU数组内存
     if (cudaStat != cudaSuccess) {
         printf ("device memory allocation failed");
         return EXIT_FAILURE;
     }
     stat = cublasCreate(&handle); // 创建cublas上下文
     if (stat != CUBLAS_STATUS_SUCCESS) {
         printf ("CUBLAS initialization failed\n");
         return EXIT_FAILURE;
     }
     stat = cublasSetMatrix (M, N, sizeof(*a), a, M, devPtrA, M); // 矩阵赋值
     if (stat != CUBLAS_STATUS_SUCCESS) {
         printf ("data download failed");
         cudaFree (devPtrA);
         cublasDestroy(handle);
         return EXIT_FAILURE;
     }
     modify (handle, devPtrA, M, N, 1, 2, 16.0f, 12.0f); // 矩阵计算
     stat = cublasGetMatrix (M, N, sizeof(*a), devPtrA, M, a, M);
     if (stat != CUBLAS_STATUS_SUCCESS) {
         printf ("data upload failed");
         cudaFree (devPtrA);
         cublasDestroy(handle);
         return EXIT_FAILURE;
     }
     cudaFree (devPtrA); // 释放GPU内存
     cublasDestroy(handle); // 销毁cublas上下文句柄
     for (j = 0; j < N; j++) {
         for (i = 0; i < M; i++) {
             printf ("%7.0f", a[IDX2C(i,j,M)]); // 打印计算结果
         }
         printf ("\n");
     }
     free(a); // 释放CPU内存
     return EXIT_SUCCESS; 

编译示例代码：

nvcc cublas_example.c -lcublas -o cublas_example

示例代码执行结果：

      1      6     11     16     21     26
      2      7     12     17     22     27
      3   1536    156    216    276    336
      4    144     14     19     24     29
      5    160     15     20     25     30

官网参考链接：https://docs.nvidia.com/cuda/cublas/index.html#introduction