使用 CUBLAS 库给矩阵运算提速

前言

　　编写 CUDA 程序真心不是个简单的事儿，调试也不方便，很费时。那么有没有一些现成的 CUDA 库来调用呢？

　　答案是有的，如 CUBLAS 就是 CUDA 专门用来解决线性代数运算的库。

　　本文将大致介绍如何使用 CUBLAS 库，同时演示一个使用 CUBLAS 库进行矩阵乘法的例子。

CUBLAS 内容

　　CUBLAS 是 CUDA 专门用来解决线性代数运算的库，它分为三个级别：

　　Lev1. 向量相乘

　　Lev2. 矩阵乘向量

　　Lev3. 矩阵乘矩阵

　　同时该库还包含状态结构和一些功能函数。

CUBLAS 用法

　　大体分成以下几个步骤：

　　1. 定义 CUBLAS 库对象

　　2. 在显存中为待运算的数据以及需要存放结果的变量开辟显存空间。( cudaMalloc 函数实现 )

　　3. 将待运算的数据传输进显存。( cudaMemcpy，cublasSetVector 等函数实现 )

　　3. 调用 CUBLAS 库函数 ( 根据 CUBLAS 手册调用需要的函数 )

　　4. 从显存中获取结果变量。( cudaMemcpy，cublasGetVector 等函数实现 )

　　5. 释放申请的显存空间以及 CUBLAS 库对象。( cudaFree 及 cublasDestroy 函数实现 )

代码示例

　　如下程序使用 CUBLAS 库进行矩阵乘法运算，请仔细阅读注释，尤其是 API 的参数说明：

// CUDA runtime 库 + CUBLAS 库 
#include "cuda_runtime.h"
#include "cublas_v2.h"

#include <time.h>
#include <iostream>

using namespace std;

// 定义测试矩阵的维度
int const M = 5;
int const N = 10;

int main() 
{   
    // 定义状态变量
    cublasStatus_t status;

    // 在 内存 中为将要计算的矩阵开辟空间
    float *h_A = (float*)malloc (N*M*sizeof(float));
    float *h_B = (float*)malloc (N*M*sizeof(float));
    
    // 在 内存 中为将要存放运算结果的矩阵开辟空间
    float *h_C = (float*)malloc (M*M*sizeof(float));

    // 为待运算矩阵的元素赋予 0-10 范围内的随机数
    for (int i=0; i<N*M; i++) {
        h_A[i] = (float)(rand()%10+1);
        h_B[i] = (float)(rand()%10+1);
    
    }
    
    // 打印待测试的矩阵
    cout << "矩阵 A :" << endl;
    for (int i=0; i<N*M; i++){
        cout << h_A[i] << " ";
        if ((i+1)%N == 0) cout << endl;
    }
    cout << endl;
    cout << "矩阵 B :" << endl;
    for (int i=0; i<N*M; i++){
        cout << h_B[i] << " ";
        if ((i+1)%M == 0) cout << endl;
    }
    cout << endl;
    
    /*
    ** GPU 计算矩阵相乘
    */

    // 创建并初始化 CUBLAS 库对象
    cublasHandle_t handle;
    status = cublasCreate(&handle);
    
    if (status != CUBLAS_STATUS_SUCCESS)
    {
        if (status == CUBLAS_STATUS_NOT_INITIALIZED) {
            cout << "CUBLAS 对象实例化出错" << endl;
        }
        getchar ();
        return EXIT_FAILURE;
    }

    float *d_A, *d_B, *d_C;
    // 在 显存 中为将要计算的矩阵开辟空间
    cudaMalloc (
        (void**)&d_A,    // 指向开辟的空间的指针
        N*M * sizeof(float)    //　需要开辟空间的字节数
    );
    cudaMalloc (
        (void**)&d_B,    
        N*M * sizeof(float)    
    );

    // 在 显存 中为将要存放运算结果的矩阵开辟空间
    cudaMalloc (
        (void**)&d_C,
        M*M * sizeof(float)    
    );

    // 将矩阵数据传递进 显存 中已经开辟好了的空间
    cublasSetVector (
        N*M,    // 要存入显存的元素个数
        sizeof(float),    // 每个元素大小
        h_A,    // 主机端起始地址
        1,    // 连续元素之间的存储间隔
        d_A,    // GPU 端起始地址
        1    // 连续元素之间的存储间隔
    );
    cublasSetVector (
        N*M, 
        sizeof(float), 
        h_B, 
        1, 
        d_B, 
        1
    );

    // 同步函数
    cudaThreadSynchronize();

    // 传递进矩阵相乘函数中的参数，具体含义请参考函数手册。
    float a=1; float b=0;
    // 矩阵相乘。该函数必然将数组解析成列优先数组
    cublasSgemm (
        handle,    // blas 库对象 
        CUBLAS_OP_T,    // 矩阵 A 属性参数
        CUBLAS_OP_T,    // 矩阵 B 属性参数
        M,    // A, C 的行数 
        M,    // B, C 的列数
        N,    // A 的列数和 B 的行数
        &a,    // 运算式的 α 值
        d_A,    // A 在显存中的地址
        N,    // lda
        d_B,    // B 在显存中的地址
        M,    // ldb
        &b,    // 运算式的 β 值
        d_C,    // C 在显存中的地址(结果矩阵)
        M    // ldc
    );
    
    // 同步函数
    cudaThreadSynchronize();

    // 从 显存 中取出运算结果至 内存中去
    cublasGetVector (
        M*M,    //  要取出元素的个数
        sizeof(float),    // 每个元素大小
        d_C,    // GPU 端起始地址
        1,    // 连续元素之间的存储间隔
        h_C,    // 主机端起始地址
        1    // 连续元素之间的存储间隔
    );
    
    // 打印运算结果
    cout << "计算结果的转置 ( (A*B)的转置 )：" << endl;

    for (int i=0;i<M*M; i++){
            cout << h_C[i] << " ";
            if ((i+1)%M == 0) cout << endl;
    }
    
    // 清理掉使用过的内存
    free (h_A);
    free (h_B);
    free (h_C);
    cudaFree (d_A);
    cudaFree (d_B);
    cudaFree (d_C);

    // 释放 CUBLAS 库对象
    cublasDestroy (handle);

    getchar();
    
    return 0;
}

运行测试

　　

　　PS：矩阵元素是随机生成的

小结

　　1. 使用 CUDA 库固然方便，但也要仔细的参阅函数手册，其中每个参数的含义都要很清晰才不容易出错。

　　2. 如果程序仅使用 CUDA 库的话，用 .cpp 源码文件即可 (不用 .cu)

 

Cublas是一个可以与cuda一同使用的函数库，它提供了多种矩阵运算的API，但是它列主序的存储方式却让人十分疑惑，今天我就以cublas中的矩阵乘法运算简单说一下我的理解。

Cublas中的矩阵乘法运算函数有5个，分别是cublasSgemm、cublasDgemm、cublasCgemm、cublasZgemm、cublasHgemm，分别包括了不同数据类型的计算，比如单精度浮点、双精度浮点、已经包括了复数等的运算。

我们查看一个具体的函数：

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);

 

在文档中对该函数的介绍：

 

 

该函数实现了矩阵运算C= αop(A)op(B)+βC，op就是指使用原矩阵运算还是使用转置之后的矩阵运算。其中m、k为第一个矩阵的行列数，k、n为第二个矩阵行列数。

初一看，并没有什么奇怪的地方，但是需要注意上边的一点：column-major format，列主序。他这个列主序是什么意思呢，下面我通过一个例子来解释一下。

我们待计算的矩阵有

 

A4*3 = [1,2,3,

     4, 5,6,

     7,8, 9,

     10,11,12]

B3*2 = [1,2,

             3,4,

             5,6]

 

若我们令m=4,k=3,n=2，令transa和transb均为CUBLAS_OP_N，那它会不会进行如下图的矩阵乘法？

图1

 

答案是不会的，我们需要注意一下前边提到的列主序。

对于一个数组A4*3= [1,2,3,4,5,6,7,8,9,10,11,12]，cublas在填充矩阵时，是按照列的顺序来填充的，比如对于数组A，若将其填充到一个4*3的矩阵中，它是这样填充的：

图2

 

所以，若参数按前边那样的话，实际进行的运算是这样的：

图3

 

那么，如果我们还是想进行如图1的运算，该怎么办呢？这就用到了transa和transb这2个参数。将这2个参数设置为CUBLAS_OP_T，它会做一个矩阵的转置，最终的结果就是如图1所示。需要注意的是，参数m、k、n是转置后矩阵的行列数。

 

上边图3中的矩阵计算结果如下图：

 

那么，该矩阵是如何存储的呢？

计算结果矩阵也遵循列主序的方式，所以，若将上图矩阵从GPU拷贝回主机会发现实际是以下面的数组的方式存放：

C = [38, 44, 50, 56, 83, 98,113,128]。