在一些诸如特征比对等等的应用场景中,特征值其实就是一个个的浮点数(float)组成的矩阵,这些浮点数每一个都是4个字节(32位),然后对这两个矩阵进行乘法计算,A * B = C,得出的C矩阵就是能够代表两个特征矩阵相似性的值。毫无疑问矩阵乘法可以在CPU中进行计算,但是,在一些需要实时性的场合下,CPU的运算速度可能就不能满足要求,而通过GPU去做矩阵乘法反而再合适不过,GPU处理浮点数运算的能力非常强,这里记录两种接口来使用GPU计算矩阵乘法。
首先,两个矩阵相乘有意义必须满足:矩阵A为m行n列,矩阵B为n行k列,则矩阵C=AB即C为m行k列。使用cuda的runtime库和cublas库在GPU中进行计算。现在CPU中使用随机浮点数产生两个矩阵h_A和h_B,再使用cublasSetVector函数或cudaMemcpy函数把h_A和h_B从CPU内存中复制到GPU中相应的d_A和d_B中,如果是float型,可以直接使用cublasSegmm函数进行矩阵乘法运算,结果置于d_C中,若为int8型,则可使用cublasGemmEx函数,再使用cublasGetVector函数或cudaMemcpy函数将d_C从GPU内存中复制到CPU内存中。
cublasSegmm函数做FP32型矩阵乘法实现:
A:M行N列矩阵 B:N行M列矩阵 C:M行M列矩阵
#include "cuda_runtime.h"
#include "cublas_v2.h"
#include <iostream>
using namespace std;
#define M 5
#define N 3
int main()
{
// 定义状态变量,矩阵乘法接口的返回值
cublasStatus_t status;
// 在 CPU内存 中为将要计算的矩阵开辟空间
float *h_A = (float*)malloc (N*M*sizeof(float));
float *h_B = (float*)malloc (N*M*sizeof(float));
// 在 CPU内存 中为将要存放运算结果的矩阵开辟空间
float *h_C = (float*)malloc (M*M*sizeof(float));
// 为待运算矩阵的元素赋予 0-10 范围内的随机数,实际使用时,A矩阵需要做转置,B矩阵不需要转 置,可以手动转置也可以算法内转置
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 // 连续元素之间的存储间隔
);
//注意:当矩阵过大时,使用cudaMemcpy是更好地选择:
//cudaMemcpy(d_A, h_A, sizeof(float)*N*M, cudaMemcpyHostToDevice);
cublasSetVector (
N*M,
sizeof(float),
h_B,
1,
d_B,
1
);
//cudaMemcpy(d_B, h_B, sizeof(float)*N*M, cudaMemcpyHostToDevice);
// 同步函数
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 // 连续元素之间的存储间隔
);
//或使用cudaMemcpy(h_C, d_C, sizeof(float)*M*M, cudaMemcpyDeviceToHost);
// 打印运算结果
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);
try
{
cudaFree (d_A);
cudaFree (d_B);
cudaFree (d_C);
}
catch(...)
{
cout << "cudaFree Error!" << endl;
// 释放 CUBLAS 库对象
}
cublasDestroy (handle);
return 0;
}
但是GPU的内存一般比较有限,实际项目中,为了能够做更大的矩阵乘法,使用FP32型数据占用空间会很大,所以可以压缩为int8型数据进行矩阵乘法运算,这时使用的接口模型为:
cublasStatus_t cublasGemmEx(cublasHandle_t handle, //句柄
cublasOperation_t transa, //a矩阵转置选项
cublasOperation_t transb, //b矩阵转置选项
int m, //a矩阵行数
int n, //b矩阵列数
int k, //a矩阵列数兼b矩阵行数
const void *alpha, //乘法因子alpha
const void *A, //A矩阵
cudaDataType_t Atype, //A矩阵的数据类型
int lda, //A矩阵的列数
const void *B, //B矩阵
cudaDataType_t Btype, //B矩阵的数据类型
int ldb, //B矩阵的行数
const void *beta, //乘法因子beta
void *C, //C结果矩阵
cudaDataType_t Ctype, //C结果矩阵数据类型
int ldc, //C矩阵的行数
cudaDataType_t computeType, //计算模式
cublasGemmAlgo_t algo)
该矩阵接口和cublasSgemm几乎相同,但多了几个参数:cudaDataType_t Atype A矩阵数据类型、cudaDataType_t Btype B矩阵数据类型、cudaDataType_t Ctype C矩阵数据类型以及cudaDataType_t computeType 计算模式和cublasGemmAlgo_t参数。首先来说,这个接口可以完成的运算表达式如下:
C=αop(A)op(B)+βC
其中,α和β为标量,A、B、C分别是m×k、k×n、m×n的矩阵,当α = 1, β = 0时,即为:
C=op(A)op(B)
也就是做矩阵乘法。α和β分别对应cublasGemmEx中的const void *alpha、const void *beta,这一点与前面cublasSegmm做FP32类型的接口是一样的。
在接口中,cublasOperation_t transa矩阵转置选项有如下三种:
transa == CUBLAS_OP_N表示A矩阵,无任何转置。
transa == CUBLAS_OP_T表示A^T,A的转置。
transa == CUBLAS_OP_C表示A^H,A的共轭转置。
同理,B矩阵也是一样。
在接口中,新增的cudaDataType_t Atype A矩阵数据类型、cudaDataType_t Btype B矩阵数据类型、cudaDataType_t Ctype C矩阵数据类型以及cudaDataType_t computeType这几个参数的取值情况如下表所示:
computeType | Atype/Btype | Ctype |
---|---|---|
CUDA_R_16F | CUDA_R_16F | CUDA_R_16F |
CUDA_R_32I | CUDA_R_8I | CUDA_R_32I |
CUDA_R_32F | CUDA_R_16F | CUDA_R_16F |
CUDA_R_8I | CUDA_R_32F | |
CUDA_R_16F | CUDA_R_32F | |
CUDA_R_32F | CUDA_R_32F | |
CUDA_R_64F | CUDA_R_64F | CUDA_R_64F |
CUDA_C_32F | CUDA_C_8I | CUDA_C_32F |
CUDA_C_32F | CUDA_C_32F | |
CUDA_C_64F | CUDA_C_64F | CUDA_C_64F |
我们这里选择的int8型矩阵乘法,最好选择Atype/Btype为CUDA_R_8I型,计算模式CUDA_R_32I,乘出的结果Ctype为CUDA_R_32I类型。最后的cublasGemmAlgo_t参数可以选择的值如下表所示:
CublasGemmAlgo_t | 含义 |
---|---|
CUBLAS_GEMM_DEFAULT | Apply Heuristics to select the GEMM algorithm |
CUBLAS_GEMM_ALGO0 to CUBLAS_GEMM_ALGO23 | Explicitly choose an algorithm |
CUBLAS_GEMM_DEFAULT_TENSOR_OP | Apply Heuristics to select the GEMM algorithm while allowing the use of Tensor Core operations if possible |
CUBLAS_GEMM_ALGO0_TENSOR_OP to CUBLAS_GEMM_ALGO15_TENSOR_OP | Explicitly choose a GEMM algorithm allowing it to use Tensor Core operations if possible, otherwise falls back to cublas<t>gemmBatched based on computeType |
这里可以选择CUBLAS_GEMM_ALGO0~7。
cublasGemmEX函数做int8型矩阵乘法实现:
#include "cuda_runtime.h"
#include "cublas_v2.h"
#include <iostream>
#include <math.h>
using namespace std;
#define M 5
#define N 3
#define BYTE 128
int main()
{
// 定义状态变量,矩阵乘法接口的返回值
cublasStatus_t status;
// 在 CPU内存 中为将要计算的矩阵开辟空间
char *h_A = (char *)malloc (N * M * sizeof(char));
char *h_B = (char *)malloc (N * M * sizeof(char));
// 在 CPU内存 中为将要存放运算结果的矩阵开辟空间
int *h_C = (int *)malloc ( M * M * sizeof(int));
float *f_A = (float *)malloc(N * M * sizeof(float));
float *f_B = (float *)malloc(N * M * sizeof(float));
// 为待运算矩阵的元素赋予 0-10 范围内的随机数,实际使用时,A矩阵需要做转置,B矩阵不需要转置,可以手动转置也可以算法内转置
for (int i = 0; i < N * M; i++)
{
f_A[i] = (float)(rand()%10+1);
f_B[i] = (float)(rand()%10+1);
}
//FP32转int8型,有精度损失
for(int i = 0; i < N * M; i++)
{
h_A[i] = (char)(round(f_A[i] * BYTE));
h_B[i] = (char)(round(f_B[i] * BYTE));
}
// 打印待测试的矩阵
/*cout << "矩阵 A :" << endl;
for (int i = 0; i < N * M; i++){
cout << (int)h_A[i] << " ";
if ((i + 1) % N == 0)
cout << endl;
}
cout << endl;
cout << "矩阵 B :" << endl;
for (int i = 0; i < N * M; i++){
cout << (int)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;
}
char *d_A, *d_B;
int *d_C;
// 在 显存 中为将要计算的矩阵开辟空间
cudaMalloc (
(void **)&d_A, // 指向开辟的空间的指针
N * M * sizeof(char) // 需要开辟空间的字节数
);
cudaMalloc (
(void **)&d_B,
N * M * sizeof(char)
);
// 在 显存 中为将要存放运算结果的矩阵开辟空间
cudaMalloc (
(void **)&d_C,
M * M * sizeof(int)
);
// 将矩阵数据传递进 显存 中已经开辟好了的空间
cublasSetVector (
N * M, // 要存入显存的元素个数
sizeof(char), // 每个元素大小
h_A, // 主机端起始地址
1, // 连续元素之间的存储间隔
d_A, // GPU 端起始地址
1 // 连续元素之间的存储间隔
);
//注意:当矩阵过大时,使用cudaMemcpy是更好地选择:
//cudaMemcpy(d_A, h_A, sizeof(char)*N*M, cudaMemcpyHostToDevice);
cublasSetVector (
N * M,
sizeof(char),
h_B,
1,
d_B,
1
);
//cudaMemcpy(d_B, h_B, sizeof(char)*N*M, cudaMemcpyHostToDevice);
// 同步函数
cudaThreadSynchronize();
// 传递进矩阵相乘函数中的参数,具体含义请参考函数手册。
float a = 1.0;
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
);
cublasGemmEx (handle, //句柄
CUBLAS_OP_T, //矩阵 A 属性参数
CUBLAS_OP_T, //矩阵 B 属性参数
M, //A, C 的行数
M, //B, C 的列数
N, //A 的列数和 B 的行数
&a, //运算式的 α 值
d_A, //A矩阵
CUDA_R_8I, //A矩阵计算模式,int8型
N, //A矩阵的列数
d_B, //B矩阵
CUDA_R_8I, //B矩阵计算模式,int8型
M, //B矩阵的行数
&b, //乘法因子beta
d_C, //C结果矩阵
CUDA_R_32I, //C矩阵计算模式,int32型
M, //C矩阵的行数
CUDA_R_32I, //计算模式,int32模式
CUBLAS_GEMM_ALGO0 //算法参数
)
// 同步函数
cudaThreadSynchronize();
// 从 显存 中取出运算结果至 内存中去
cublasGetVector (
M * M, // 要取出元素的个数
sizeof(int), // 每个元素大小
d_C, // GPU 端起始地址
1, // 连续元素之间的存储间隔
h_C, // 主机端起始地址
1 // 连续元素之间的存储间隔
);
//或使用cudaMemcpy(h_C, d_C, sizeof(int)*M*M, cudaMemcpyDeviceToHost);
// 打印运算结果
cout << "计算结果的转置 ( (A*B)的转置 ):" << endl;
for (int i = 0;i < M * M; i++)
{
cout << h_C[i] << " ";
if ((i+1)%M == 0)
cout << endl;
}
//注意,这里需要进行归一化操作,乘出来的矩阵需要除以128*128,以还原原来的大小。在此就省略这一步。
// 清理掉使用过的内存
free (h_C);
free (h_B);
free (h_A);
free (f_B);
free (f_A);
try
{
cudaFree (d_A);
cudaFree (d_B);
cudaFree (d_C);
}
catch(...)
{
cout << "cudaFree Error!" << endl;
// 释放 CUBLAS 库对象
}
cublasDestroy (handle);
return 0;
}
以上代码编译命令可以是:g++ GPURetrievalInt8.cpp -o GPURetrievalDemo -L/usr/local/cuda/lib64 -lcudart -lcuda
当运行有误时,可以打印cublasGemmEx的返回值来确认是什么错误,官方定义的错误码如下表所示:
Error Value | Meaning |
---|---|
CUBLAS_STATUS_SUCCESS | the operation completed successfully |
CUBLAS_STATUS_NOT_INITIALIZED | the library was not initialized |
CUBLAS_STATUS_ARCH_MISMATCH | cublasGemmEx is only supported for GPU with architecture capabilities equal or greater than 5.0 |
CUBLAS_STATUS_NOT_SUPPORTED | the combination of the parameters Atype, Btype and Ctype or the algorithm, algo is not supported |
CUBLAS_STATUS_INVALID_VALUE | the parameters m,n,k<0 |
CUBLAS_STATUS_EXECUTION_FAILED | the function failed to launch on the GPU |
成功时会返回0,也就是CUBLAS_STATUS_SUCCESS。多数情况会返回CUBLAS_STATUS_INVALID_VALUE,这就是m,n,k参数传入不正确所致。
经测试,GPU矩阵乘法可以达到大约1.6亿次/秒的速度,只要提前将向量存入GPU显存中,待比对向量以1个或n个的形式拷贝到GPU显存中,乘法之后选用合适的排序方法如堆排序等,就可以从大到小将相似度最高到最低排列出来。应用前景广泛。