C++矩阵乘法计算 || GPU && CPU 实现_c++矩阵乘法库

前言

---

矩阵乘法运算是机器学习的基础。比如，卷积神经网络通过矩阵化输入数据，然后通过矩阵乘法计算获得结果。而性能对于算法是至关重要的事情，所以本文主要介绍c++调用普通的矩阵乘法库进行计算，以及通过cuda计算矩阵乘法。C++常用cblas库加速cpu上的矩阵乘法运算。为了将速度提高更高，GPU版本矩阵乘法运算则通过cublas库进行操作，在cublas库中，使用cublasSgemv()和cublasSgemm()分别进行矩阵向量间的乘法运算与矩阵矩阵间的乘法运算。本文将具体的解释上述两个函数的参数以及具体的应用例子。参照官方解释

矩阵乘法函数解释

---

如何安装cpu版本和gpu版本的相应矩阵乘法库，可以参考这篇文章Ubuntu16.04 + Caffe + CUDA9.0 + cudnn7.0 的配置详细教程&& Ubuntu18.04 可用，虽然说是安装caffe框架的，但caffe则使用了两者的矩阵乘法库都用到了

CPU版本矩阵矩阵乘法

cblas有两个函数实现矩阵乘法，一个是cblas_sgemm(),另一个是cblas_dgemm()两者的不同点在于传入参数一个是float型，一个是double型。

void cblas_sgemm(
        OPENBLAS_CONST enum CBLAS_ORDER Order,         // 矩阵存储形式，行优先或者列优先
        OPENBLAS_CONST enum CBLAS_TRANSPOSE TransA,    // 进行矩阵乘运算前，A是否转置
        OPENBLAS_CONST enum CBLAS_TRANSPOSE TransB,    // 进行矩阵运算前，B是否转置
        OPENBLAS_CONST blasint M,                      // A的行数
        OPENBLAS_CONST blasint N,                      // B的列数
        OPENBLAS_CONST blasint K,                      // A的列数 <==> B的行数 
		OPENBLAS_CONST double alpha,                   // 比例因子
		OPENBLAS_CONST double *A,                      // A的首地址
		OPENBLAS_CONST blasint lda,                    // A的列数，与是否转置无关
		OPENBLAS_CONST double *B,                      // B的首地址
		OPENBLAS_CONST blasint ldb,                    // B的列数，与转置无关
		OPENBLAS_CONST double beta,                    // 比例因子
		double *C,                                     // C的首地址
		OPENBLAS_CONST blasint ldc);                   // C的列数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

• 计算结果: C = alpha * TransA(A) x TransB(B) + beta * C
  其中, A为M x K矩阵, B为 K x N矩阵，C 为 M x N 矩阵

CPU版本矩阵向量乘法

void cblas_sgemv(
        OPENBLAS_CONST enum CBLAS_ORDER order,       // A的数据存储型式  
        OPENBLAS_CONST enum CBLAS_TRANSPOSE trans,   // A进行矩阵运算前，是否转置
        OPENBLAS_CONST blasint M,                    // 转置后, A的行维度
        OPENBLAS_CONST blasint N,                    // 转置后, A的列维度
		OPENBLAS_CONST float alpha,                  // 比例因子
		OPENBLAS_CONST float  *A,                    // A的首地址
		OPENBLAS_CONST blasint lda,                  // A的列数，与转置无关
		OPENBLAS_CONST float  *X,                    // X的首地址
		OPENBLAS_CONST blasint incx,                 // X的首地址大小
		OPENBLAS_CONST float beta,                   // 比例因子
		float  *Y,                                   // Y的首地址 
		OPENBLAS_CONST blasint incy);                // Y的列数，与转置无关

1
2
3
4
5
6
7
8
9
10
11
12
13
14

运算式子：Y = alpha*trans(A) x X + beta * C
其中，A为M x N的矩阵，而X为N x 1向量，C为M x 1 的向量

• 可以看出，函数需要传入的参数比较多，第一个参数是矩阵的存储型式(可选行优先或者列优先)，

GPU版本矩阵乘法

仅仅使用单精度(float)进行解释，若想尝试双精度(double),可以进入官网查看相应函数.

cublasStatus_t cublasSgemv(cublasHandle_t handle,            //传入的cublas句柄 
                           cublasOperation_t trans,          
                           int m, int n,
                           const float           *alpha,
                           const float           *A, int lda,
                           const float           *x, int incx,
                           const float           *beta,
                           float           *y, int incy)
1
2
3
4
5
6
7
8

这个函数实现的是矩阵向量乘法运算:
y = alpha*op(A) x + beta *y
A是一个mxn维的矩阵，以列优先的形式存储在内存中，x,y分别是向量，alpha,beta分别规模因子，也就是倍数关系，而对于矩阵A来说 op(A) = A 如果第二个参数 transa == CUBLAS_OP_N，如果第二个参数为transa == CUBLAS_OP_T，那么op(A)结果是A的转置，而如果transa==cublas_op_H，那么op(A)的结果是A的共轭转置.
而lda是A矩阵主轴的维度，而incx与incy分别是x,y各个元素之间的间隔。

cublasStatus_t cublasSgemm(cublasHandle_t handle,       //cublas句柄
                           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)
1
2
3
4
5
6
7
8
9

C = alpha * op(A)op(B) + beta*C，其中A,B,C都为矩阵，而A维度为mxk，B维度为kxn，而C的维度为mxk，lda,ldb,ldc分别是矩阵A,B,C主轴上的维度。

#实例

---

#include "./common/common.h"
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include "cublas_v2.h"
#include <iostream>
using namespace std;

/*
 * A simple example of performing matrix-vector multiplication using the cuBLAS
 * library and some randomly generated inputs.
 */

/*
 * M = # of rows
 * N = # of columns
 */
int M = 5;
int N = 2;

/*
 * Generate a vector of length N with random single-precision floating-point
 * values between 0 and 100.
 */
void generate_random_vector(int N, float **outX)
{
    int i;
    double rMax = (double)RAND_MAX;
    float *X = (float *)malloc(sizeof(float) * N);

    for (i = 0; i < N; i++)
    {
        int r = rand();
        double dr = (double)r;
        X[i] = (dr / rMax) * 100.0;
    }

    *outX = X;
}

/*
 * Generate a matrix with M rows and N columns in column-major order. The matrix
 * will be filled with random single-precision floating-point values between 0
 * and 100.
 */
void generate_random_dense_matrix(int M, int N, float **outA)
{
    int i, j;
    double rMax = (double)RAND_MAX;
    float *A = (float *)malloc(sizeof(float) * M * N);

    // For each column
    for (j = 0; j < N; j++)
    {
        // For each row
        for (i = 0; i < M; i++)
        {
            double dr = (double)rand();
            A[j * M + i] = (dr / rMax) * 100.0;
        }
    }

    *outA = A;
}

int main(int argc, char **argv)
{
    int i;
    float *A, *dA;
    float *X, *dX;
    float *Y, *dY;
    float beta;
    float alpha;
    cublasHandle_t handle = 0;

    alpha = 3.0f;
    beta = 4.0f;

    // Generate inputs
    srand(9384);
    generate_random_dense_matrix(M, N, &A);
    cout<<"A: "<<endl;
    for(int i = 0; i < M; i ++){
        for(int j = 0; j < N; j ++){
            cout << A[M * j + i] <<" ";
        }
        cout<< endl;
    }


    generate_random_vector(N, &X);

    cout<<"X: "<<endl;
    for(int i = 0; i < N; i ++){
        cout<<X[i]<<" ";
    }
    cout<< endl;


    generate_random_vector(M, &Y);

    cout<<"Y: "<<endl;
    for(int i = 0; i < M; i++){
        cout<<Y[i]<<" ";
    }
    cout<<endl;

    // Create the cuBLAS handle
    CHECK_CUBLAS(cublasCreate(&handle));

    // Allocate device memory
    CHECK(cudaMalloc((void **)&dA, sizeof(float) * M * N));
    CHECK(cudaMalloc((void **)&dX, sizeof(float) * N));
    CHECK(cudaMalloc((void **)&dY, sizeof(float) * M));


    // Transfer inputs to the device
    CHECK_CUBLAS(cublasSetVector(N, sizeof(float), X, 1, dX, 1));
    CHECK_CUBLAS(cublasSetVector(M, sizeof(float), Y, 1, dY, 1));
    CHECK_CUBLAS(cublasSetMatrix(M, N, sizeof(float), A, M, dA, M));

    // Execute the matrix-vector multiplication
    CHECK_CUBLAS(cublasSgemv(handle, CUBLAS_OP_N, M, N, &alpha, dA, M, dX, 1,
                             &beta, dY, 1));
    cout<< "aplha: "<<alpha<< " beta: "<<beta<<endl;
    // Retrieve the output vector from the device
    CHECK_CUBLAS(cublasGetVector(M, sizeof(float), dY, 1, Y, 1));

    cout<<"result: "<<endl;
    for (i = 0; i < M; i++)
    {
        printf("%2.2f\n", Y[i]);
    }

    printf("...\n");

    free(A);
    free(X);
    free(Y);

    CHECK(cudaFree(dA));
    CHECK(cudaFree(dY));
    CHECK_CUBLAS(cublasDestroy(handle));

    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

最终的输出结果为:

A: 
51.3409 68.9939 
27.9576 78.9213 
52.4111 97.6222 
36.8833 74.3454 
39.9056 4.77248 
B: 
82.3203 91.5174 84.3425 
30.7926 37.0584 81.2048 
C: 
34.009 36.5475 79.0337 
7.12583 76.2508 76.6652 
71.0715 39.7865 79.2842 
3.62919 88.8859 93.6551 
78.67 60.1165 20.3579 
X: 
30.6251 21.6128 
72.769 67.5084 61.5184 41.7629 46.4297 
Y: 
9481.49 7955.78 11391 8376.15 4161.5 
C: 
19188.7 21912.4 30114.7 
14223.5 16754.9 26607.1 
22245.9 25401.9 37360.8 
15991.1 18747.3 27818.7 
10610.7 11727.2 11341.3 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26