2. 程式人生 > >C++矩陣乘法計算 || GPU && CPU 實現

C++矩陣乘法計算 || GPU && CPU 實現

阿新 • • 發佈：2019-01-30

前言

矩陣乘法運算是機器學習的基礎。比如，卷積神經網路通過矩陣化輸入資料，然後通過矩陣乘法計算獲得結果。而效能對於演算法是至關重要的事情，所以本文主要介紹c++呼叫普通的矩陣乘法庫進行計算，以及通過cuda計算矩陣乘法。C++常用cblas庫加速cpu上的矩陣乘法運算。為了將速度提高更高，GPU版本矩陣乘法運算則通過cublas庫進行操作，在cublas庫中，使用cublasSgemv()和cublasSgemm()分別進行矩陣向量間的乘法運算與矩陣矩陣間的乘法運算。本文將具體的解釋上述兩個函式的引數以及具體的應用例子。參照官方解釋



矩陣乘法函式解釋

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

相關推薦

C++矩陣乘法計算 || GPU && CPU 實現

前言 矩陣乘法運算是機器學習的基礎。比如，卷積神經網路通過矩陣化輸入資料，然後通過矩陣乘法計算獲得結果。而效能對於演算法是至關重要的事情，所以本文主要介紹c++呼叫普通的矩陣乘法庫進行計算，以及通過cuda計算矩陣乘法。C++常用cblas庫加速cpu上的矩陣

線性代數教程之一——矩陣乘法計算、理解及程式碼實現

參考了《深度學習》鉅作，以下是矩陣篇的目錄。 1 矩陣的乘法 設矩陣A為m×n矩陣，B為n×p矩陣，則它們的乘法公式為： 相關程式碼實現： # 矩陣滴乘法運算 # 注意：需要傳入np.matrix型別資料 def Matrix_Mul(a,b):

C語言 乘法計算

#include <stdio.h> #include <stdlib.h> /* run this program using the console pauser or

【華為機試070】矩陣乘法計算量估算

題目描述： 矩陣乘法的運算量與矩陣乘法的順序強相關。 例如：     A是一個50×10的矩陣，B是10×20的矩陣，C是20×5的矩陣   計算A*B*C有兩種順序：（（AB）C）或者（A（BC）），前者需要計算15000

大矩陣乘法運算map reduce實現思路

實現思路： 儲存： 大矩陣很多都是稀疏矩陣，並且有可能有上百萬的行和上百萬的列。 那麼矩陣可以存在類似HBase面向列的分散式資料庫中。 假設HTable中有兩個表A和表B分別儲存兩個巨型矩陣a和b。表A和表B都是隻有一個列族。列名都是1開始計數。 那麼表A和表B所儲存的矩

MPI矩陣乘法的兩種實現方法

MPI矩陣乘法 去年學習了平行計算，接觸了MPI、Pthreads和OpenMP等常用的並行方法實現了並行的矩陣乘法，本章在此總結一下MPI的矩陣乘法使用。 使用簡單的MPI_Send和MPI_Recv實現 使用較高階的MPI_Scatter和MPI_Ga

C#矩陣乘法

/// <summary> /// <summary> /// 矩陣乘法 /// </summary>

華為OJ——矩陣乘法計算量估算

矩陣乘法計算量估算 題目描述 矩陣乘法的運算量與矩陣乘法的順序強相關。例如：     A是一個50×10的矩陣，B是10×20的矩陣，C是20×5的矩陣 計算A*B*C有兩種順序：（（AB）C）或者（A（BC）），前者需要計算15000次乘法，後者只需要3500次。 編寫程

C語言之兩矩陣乘法的實現

首先我們要清楚矩陣乘法實現需要滿足的條件， 矩陣相乘最重要的方法是一般矩陣乘積。它只有在第一個矩陣的列數（column）和第二個矩陣的行數（row）相同時才有意義[1] 。一般單指矩陣乘積時，指的便是一般矩陣乘積。一個m×n的矩陣就是m×n個數排成m行n列的一

C語言矩陣乘法（指標實現）

這是C和指標書上的一道題，充分體現了指標實現陣列操作的過程 void matrix_multiply( int *m1, int *m2, int *r, int x, int y, int z ) { regist

遞推式計算與矩陣乘法 C++

以斐波那契數列為例，計算斐波那契數列的矩陣方法程式碼如下： （快速冪、矩陣乘法）《挑戰程式設計競賽第二版》P199 #include<iostream> #include<vect

一些演算法的MapReduce實現——矩陣分塊乘法計算(1)

矩陣分塊簡介 一個分塊矩陣(分段矩陣)就是將矩陣分割出較小的矩形矩陣，這些較小的矩陣就稱為區塊。換個方式來說，就是以較小的矩陣組合成一個矩陣。通過將大的矩陣通過分塊的方式劃分，並將每個分塊(稱為子塊)看做另一個矩陣的元素，這樣之後再參與運算，通常可以簡化運算。例如，有的大

模板C++ 02數論算法 4矩陣乘法

矩陣 快速冪 行數 正方形 eof str memset isp images 矩陣乘法：用來求某種 遞推關系。 矩陣相乘只有在第一個矩陣的列數和第二個矩陣的行數相同時才有意義。 定義 設A為A*M的矩陣，B為M*B的矩陣，那麽矩陣C為矩陣A與B的乘積，其中矩陣C中的第i行

分治法實現矩陣乘法

name cout namespace size cas put 分治 ade add 整體的思路就是分，加&乘，拼 #include <iostream> #include <cstddef> #include <cstdlib&g

(Matlab)GPU計算及CPU計算能力的比較

new start atl border ble 大小 uri win courier %%首先以200*200的矩陣做加減乘除 做比較 t = zeros(1,100); A = rand(200,200);B = rand(200,200);C = rand(200,

矩陣乘法來加速遞推式計算

span code 分享 pow quic mage src .com image Codevs1281: 給你6個數，m, a, c, x0, n, g Xn+1 = ( aXn + c ) mod m，求Xn 計算遞推式，運用矩陣來進行計算加速 然後註意用類似快速冪的

C++ 方矩陣乘法 + Strassen矩陣

  這幾天看演算法導論，看到矩陣一章，就實現了一下。 下面是普通的矩陣乘法，複雜度為：n^3。 template<unsigned M,unsigned N, unsigned Q> void Square_matrix_multiply(int(&A)[M][

矩陣乘法模板C/C++

由於矩陣相乘是諸多演算法中的基礎，就像乘法之於快速冪一樣，所以有必要熟練掌握，以下模板是我個人覺得用起來比較方便的，但是隨著時間遷移肯定是需要對其進行精進的。 程式碼示例： struct Matrix{ int n,m; int v[maxn][maxn]; Matrix(int

神經網路高效能運算 卷積計算優化 openblas GEMM 矩陣乘法優化 ncnn mobileNet-ssd shueezeNet-ssd

HighPerformanceComputing 高效能運算(High performance computing， 縮寫HPC) 指通常使用很多處理器（作為單個機器的一部分） 或者某一叢集中組織的幾臺計算機（作為單個計 算資源操作）的計算系統和環境。 有許多型別的HP

java實現兩個矩陣乘法 有個錯誤希望有大佬幫忙

//java實現兩個矩陣相乘 有個錯誤在下邊 有沒有哪個大佬幫我看看 十分感謝 package 實驗五; import java.util.Scanner; public class Matrix { private int rows; private int cols;