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

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/qq_31261509/article/details/82354015

前言


矩阵乘法运算是机器学习的基础。比如,卷积神经网络通过矩阵化输入数据,然后通过矩阵乘法计算获得结果。而性能对于算法是至关重要的事情,所以本文主要介绍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的列数

  • 计算结果: 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的列数,与转置无关

运算式子: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)

这个函数实现的是矩阵向量乘法运算:
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)

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

最终的输出结果为:

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 
展开阅读全文

没有更多推荐了,返回首页