CUBLAS库入门教程（从环境配置讲起）

---

前言

CUBLAS库是NVIDIA CUDA用于线性代数计算的库。使用CUBLAS库的原因是我不想去直接写核函数。
（当然，你还是得学习核函数该怎么写。但是人家写好的肯定比我自己写的更准确！）

---

一、搭建环境

1. 安装CUDA库，具体可以看我上一篇文章：在C++项目中集成CUDA程序加速（从环境配置讲起）
2. 如果你是装在默认路径下，那么 CUBLAS库的头文件就在：C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\include 路径下面，以cublas开头的.h文件。
3. 所以，还是按照步骤1的文章进行环境配置，然后只需要多在添加依赖项中增加一个cublas.lib就可以了。

二、简单介绍

1. CUBLAS Introdution 是官方文档。（全英文的，还有不少数学公式。大家有不理解的可以直接留言区问相关API，我们一起讨论学习。）
2. CUBLAS Samples 是官方示例，所有API都有。
3. 对于API名称，都是cublasl<t>...，其中有下述类型选择：
   
4. CUBLAS库的矩阵是列向量的，跟glm一致。
5. CUBLAS对于矩阵或者向量的index是从1开始的。所以，如果有函数的返回结果是个index（比如查找矩阵中的最大值），记得要index - 1才是我们要的。

---

三、 具体例子

下面我以矩阵与向量相乘的函数进行举例，看看是怎么用的。

1. 首先，通过查找官方文档，知道是如下的函数：

cublasStatus_t cublasDgemv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n,
                           const double          *alpha,
                           const double          *A, int lda,
                           const double          *x, int incx,
                           const double          *beta,
                           double          *y, int incy)
/*
* handle		: CUBLAS的句柄，用以管理CUBLAS库的上下文和资源
* CUBLAS_OP_N	: 指定矩阵操作模式。CUBLAS_OP_N代表正常模式（列向量）；CUBLAS_OP_T代表转置模式（行向量）
* m				: 矩阵A的行数
* n				: 矩阵A的列数
* alpha			: 与矩阵A相乘的标量
* A				: 指向存储在device上面的矩阵数据指针
* lda			: 矩阵的行数，代表矩阵在内存中的存储方式
* x				: 向量X
* incx			: 向量x中相邻两个元素的index间隔，一般为1
* beta			: 与向量y相乘的标量
* y				: 向量y
* incy			: 向量y中相邻两个元素的index间隔，一般为1
*/

具体计算公式如下：

2. 如果我们只是想计算矩阵和向量相乘，那么我们只需要令 α = 1.0, β = 0.0，然后传入我们要的A和x就行了。
3. 最后，具体代码如下：

/// MyCublas.cuh

#pragma once

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cublas_v2.h"

extern "C" void MatrixMulVectorCublas(
	const double* matrix, const int row, const int col, 
	const double* vector, double* result
);

/// MyCublas.cu

#include "MyCublas.cuh"
#include "CublasUtility.h"
#include <vector>

void MatrixMulVectorCublas(
    const double* matrix, const int row, const int col,
    const double* vector, double* result
)
{
    /* step 1: create cublas handle, bind a stream */
    cublasHandle_t cublasH = NULL;
    cudaStream_t stream = NULL;

    CUBLAS_CHECK(cublasCreate(&cublasH));

    CUDA_CHECK(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));
    CUBLAS_CHECK(cublasSetStream(cublasH, stream));

    /* step 2: copy data to device */
    double* dev_matrix = nullptr;
    double* dev_vector = nullptr;
    double* dev_result = nullptr;

    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&dev_matrix), sizeof(double) * row * col));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&dev_vector), sizeof(double) * col));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&dev_result), sizeof(double) * row));

    CUDA_CHECK(cudaMemcpyAsync(dev_matrix, matrix, sizeof(double) * row * col, cudaMemcpyHostToDevice,
        stream));
    CUDA_CHECK(cudaMemcpyAsync(dev_vector, vector, sizeof(double) * col, cudaMemcpyHostToDevice,
        stream));

    /* step 3: compute */
    cublasOperation_t transa = CUBLAS_OP_N;
    const int lda = row;
    const double alpha = 1.0;
    const double beta = 0.0;
    const int incx = 1;
    const int incy = 1;

    CUBLAS_CHECK(
        cublasDgemv(cublasH, transa, row, col, &alpha, dev_matrix, lda, dev_vector, incx, &beta, dev_result, incy));

    /* step 4: copy data to host */
    CUDA_CHECK(cudaMemcpyAsync(result, dev_result, sizeof(double) * row, cudaMemcpyDeviceToHost,
        stream));

    CUDA_CHECK(cudaStreamSynchronize(stream));

    /* step 5: free resources */
    CUDA_CHECK(cudaFree(dev_matrix));
    CUDA_CHECK(cudaFree(dev_vector));
    CUDA_CHECK(cudaFree(dev_result));

    CUBLAS_CHECK(cublasDestroy(cublasH));
    CUDA_CHECK(cudaStreamDestroy(stream));
    CUDA_CHECK(cudaDeviceReset());
}

/// CublasUtility.h

#pragma once
#include <string>
#include <stdexcept>
// CUDA API error checking
#define CUDA_CHECK(err)                                                                            \
    do {                                                                                           \
        cudaError_t err_ = (err);                                                                  \
        if (err_ != cudaSuccess) {                                                                 \
            std::printf("CUDA error %d at %s:%d\n", err_, __FILE__, __LINE__);                     \
            throw std::runtime_error("CUDA error");                                                \
        }                                                                                          \
    } while (0)

// cublas API error checking
#define CUBLAS_CHECK(err)                                                                          \
    do {                                                                                           \
        cublasStatus_t err_ = (err);                                                               \
        if (err_ != CUBLAS_STATUS_SUCCESS) {                                                       \
            std::printf("cublas error %d at %s:%d\n", err_, __FILE__, __LINE__);                   \
            throw std::runtime_error("cublas error");                                              \
        }                                                                                          \
    } while (0)

/// main.cpp

#include "MyCublas.cuh"

#include <iostream>

int main()
{
	/*
	 *   matrix   = | 1.0, 2.0, 3.0, 4.0  |
	 *              | 5.0, 6.0, 7.0, 8.0  |
	 *              | 9.0, 10.0,11.0,12.0 |
     * 
	 *   vector   = | 1.0, 2.0, 3.0, 4.0  |
     * 
     *   result   = | 30.0, 70.0, 11.0 |
	 */

    const int row = 3;
    const int col = 4;

    double matrix[row * col] = {
        1.0, 5.0, 9.0,
        2.0, 6.0, 10.0,
        3.0, 7.0, 11.0,
        4.0, 8.0, 12.0
    };

    double vector[col] = { 1.0, 2.0, 3.0, 4.0 };

    double result[row] = { 0.0 };

    MatrixMulVectorCublas(matrix, row, col, vector, result);

    for (int i = 0; i < row; ++i)
    {
        printf("%.1f, ", result[i]);
    }
	return 0;
}

四、注意

1. 如果你的matrix是行向量（row * col），那么你如果使用我上面的代码，你就要先将matrix写成列向量的形式，但是输入的row和col 仍然保持行向量的row和col。

五、疑问

对于上述代码，我还有以下的疑问：

1. 我在运行下面这句的时候，VS显示我的进程内存会到2.2GB左右，难道真的需要这么大吗？

CUBLAS_CHECK(cublasCreate(&handle));