Linux 桌面玩家指南：16. 使用 CUDA 發揮顯卡的計算性能

原文: Linux 桌面玩家指南：16. 使用 CUDA 發揮顯卡的計算性能

特別說明：要在個人隨筆後寫評論的小夥伴們請注意了，個人博客開啓了 MathJax 數學公式支持，MathJax 使用$標記數學公式的開始和結束。若是某條評論中出現了兩個$，MathJax 會將兩個$之間的內容按照數學公式進行排版，從而致使評論區格式混亂。若是你們的評論中用到了$，可是又不是爲了使用數學公式，就請使用\$轉義一下，謝謝。

想從頭閱讀該系列嗎？下面是傳送門：

• Linux 桌面玩家指南：01. 玩轉 Linux 系統的方法論 【約 1.1 萬字，22 張圖片】
• Linux 桌面玩家指南：02. 以最簡潔的方式打造實用的 Vim 環境 【約 0.97 萬字，7 張圖片】
• Linux 桌面玩家指南：03. 針對 Gnome 3 的 Linux 桌面進行美化 【約 0.58 萬字，32 張圖片】
• Linux 桌面玩家指南：04. Linux 桌面系統字體配置要略 【約 1.2 萬字，34 張圖片】
• Linux 桌面玩家指南：05. 發博客必備的圖片處理和視頻錄製神器 【約 0.25 萬字，14 張圖片】
• Linux 桌面玩家指南：06. 優雅地使用命令行及 Bash 腳本編程語言中的美學與哲學 【約 1.4 萬字，16 張圖片】
• Linux 桌面玩家指南：07. Linux 中的 Qemu、KVM、VirtualBox、Xen 虛擬機體驗 【約 0.95 萬字，31 張圖片】
• Linux 桌面玩家指南：08. 使用 GCC 和 GNU Binutils 編寫能在 x86 實模式運行的 16 位代碼 【約 0.6 萬字，17 張圖片】
• Linux 桌面玩家指南：09. X Window 的奧祕 【約 0.5 萬字，14 張圖片】
• Linux 桌面玩家指南：10. 沒有 GUI 的時候應該怎麼玩 【約 0.5 萬字，32 張圖片】
• Linux 桌面玩家指南：11. 在同一個硬盤上安裝多個 Linux 發行版以及爲 Linux 安裝 Nvidia 顯卡驅動 【約 0.4 萬字，22 張圖片】
• Linux 桌面玩家指南：12. 優秀的文本化編輯思想大碰撞（Markdown、LaTeX、MathJax） 【約 0.4 萬字，8 張圖片】
• Linux 桌面玩家指南：13. 使用 Git 及其 和 Eclipse 的集成 【約 0.6 萬字，31 張圖片】
• Linux 桌面玩家指南：14. 數值計算和符號計算 【約 0.6 萬字，15 張圖片】
• Linux 桌面玩家指南：15. 深度學習能夠這樣玩 【約 0.2 萬字，10 張圖片】

前言

科學計算碰到數據量很大的時候，每每很是消耗時間。使用多核進行並行計算是加快計算速度的主要方法，而顯卡天生具備成百上千個計算核心，因此使用 GPU 進行計算也就愈來愈流行。得益於 Nvidia 提供的 CUDA，咱們編寫利用 GPU 進行計算的程序愈來愈方便。那麼，在 Linux 系統下，使用 CUDA 究竟效果如何呢？這一篇將爲你揭曉答案。

我測試的是兩個 2000×2000 矩陣相乘所耗費的時間，測試環境是個人 Dell XPS 15，這是一款 2015 年的產品了，CPU 是 Intel 的 6700H，4 核 8 線程，顯卡是 Nvidia 的 GTX 960M，和如今 2019 年的主流比起來，已經落後幾代了。前段時間「核彈」廠的老黃又發佈了新的甜點級的顯卡，RTX 2060，流處理器又多了幾倍，還支持光線追蹤，害得我內心又長了草。

繼續說回矩陣相乘，我先測試了一下本身用 C 語言寫一個簡單的利用三重循環計算矩陣相乘的程序，而後使用 gcc 編譯，運行，進行計時。其次，我本身寫一段 CUDA 代碼，利用顯卡計算兩個矩陣相乘，使用 nvcc 編譯運行，進行計時。最後，使用 Nvidia 提供的 cuBLAS 庫中的函數直接計算兩個矩陣相乘，並進行計時。其中我本身用 C 語言寫的三重循環是徹底沒有優化的，因此計算時間確定比較慢，若是進行充分優化，利用 CPU 的 SSE、AVX 等向量化指令，並優化內存訪問以充分利用 CPU 的緩存，能夠將性能提高大約 10 倍。在雷鋒網上有一篇 OpenBLAS 項目創始人和主要維護者張先軼介紹的OpenBLAS項目與矩陣乘法優化，能夠一看。我用 Python 的 NumPy 測試了一下，確實可提速 8 到 10 倍，由於 NumPy 的底層調用了 OpenBLAS 庫。並且這還只是利用了 CPU 的一個核，若是要利用 CPU 進行並行計算，能夠考慮 OpenMP 或 MPI。

先貼代碼和結論

完整的 C 語言代碼我放在cublas_test.cu文件中，大家沒看錯，擴展名爲.cu，由於 CUDA 的編譯器 nvcc 要求這樣，nvcc 對代碼進行初步處理後，仍是調用的 gcc 進行編譯鏈接生成可執行文件。完整代碼以下：

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <sys/time.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>

#define M 2000
#define N 2000
#define K 2000

#define idx(i, j, m) ((j) * (m) + (i))

typedef struct _matrix {
      size_t height;
      size_t width;
      float *elements;
} Matrix;

Matrix empty(size_t m, size_t n){
    Matrix temp = {m, n, NULL};
    temp.elements = (float*)malloc( m*n*sizeof(float) );
    return temp;
}

Matrix zeros(size_t m, size_t n){
    Matrix temp = empty(m, n);
    memset(temp.elements, 0, sizeof( m*n*sizeof(float) ));
    return temp;
}

Matrix rands(size_t m, size_t n){
    Matrix temp = empty(m, n);
    for(size_t i=0; i<m*n; i++){
        temp.elements[i] = (float)rand()/RAND_MAX;
    }
    return temp;
}

void deleteMatrix(Matrix *mat){
    mat->height = 0;
    mat->width = 0;
    free(mat->elements);
    mat->elements = NULL;
}


void printRow(Matrix mat, size_t i){
    size_t m = mat.height;
    size_t n = mat.width;
    if(n < 7){
        for(size_t j=0; j<n; j++){
            printf(" %7.4f ", mat.elements[idx(i, j, m)]);
        }
    }else{
        for(size_t j=0; j<3; j++){
            printf(" %7.4f ", mat.elements[idx(i, j, m)]);
        }
        printf("   ...   ");
        for(size_t j=n-3; j<n; j++){
            printf(" %7.4f ", mat.elements[idx(i, j, m)]);
        }
    }
    printf("\n");
}

void printMatrix(Matrix mat){
    size_t m = mat.height;
    if(m < 7){
        for(size_t i=0; i<m; i++){
            printRow(mat, i);
        }
    }else{
        for(size_t i=0; i<3; i++){
            printRow(mat, i);
        }
        printf("   ...   \n");
        for(size_t i=m-3; i<m; i++){
            printRow(mat, i);
        }
    }
    printf("\n");
}

Matrix matMul(Matrix A, Matrix B){
    size_t m, n, k;
    Matrix C = {0, 0, NULL};
    if(A.width == B.height){
        m = A.height;
        n = B.width;
        k = B.height;
        C = zeros(m, n);
        for(size_t i=0; i<m; i++){
            for(size_t j=0; j<n; j++){
                float sum = 0;
                for(size_t p=0; p<k; p++){
                    sum += A.elements[idx(i, p, m)] * B.elements[idx(p, j, k)];
                }
                C.elements[idx(i, j, m)] = sum;
            }
        }
    }else{
        printf("Matrix shape error!\n");
    }
    return C;
}

long getTimer(struct timeval start, struct timeval end){
  return (end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec);
}

void printTimer(long timer){
  printf("%ld,%ld,%ld us\n", timer/1000000, (timer/1000)%1000, timer%1000);
}
    

bool initCUBLAS(cublasHandle_t *handle){
    int count;
    cudaGetDeviceCount(&count);
    if(count == 0){
        printf("There is no device.\n");
        return false;
    }else{
        printf("There is %d device.\n", count);
    }
    int i;
    for(i=0; i<count; i++){
        cudaDeviceProp prop;
        if(cudaGetDeviceProperties(&prop, i) == cudaSuccess){
            if(prop.major >= 1){
                break;
            }
        }
    }

    if(i == count){
        printf("There is no device supporting CUDA 1.x.\n");
        return false;
    }
    cudaSetDevice(i);

    cublasStatus_t stat;
    stat = cublasCreate(handle);
    if(stat!=CUBLAS_STATUS_SUCCESS){
        printf("CUBLAS initialization failed\n");
        return false;
    }
    printf("CUBLAS initialized.\n");

    return true;
}

__global__ void cudaMatMul(Matrix devA, Matrix devB, Matrix devC){
  size_t m = devA.height;
  size_t k = devB.height;
  size_t n = devB.width;

  size_t j = blockIdx.y * blockDim.y + threadIdx.y;
  size_t i = blockIdx.x * blockDim.x + threadIdx.x;

  if(i<m && j<n){
    float sum = 0.0;
    for(size_t p=0; p<k; p++){
      sum += devA.elements[idx(i, p, m)] * devB.elements[idx(p, j, k)];
    }
    devC.elements[idx(i, j, m)] = sum;
  }
}


int main(int argc, char** argv){
    struct timeval start;
    struct timeval end;

    //第一步，建立兩個隨機矩陣，使用 CPU 計算它們相乘，並計時
    gettimeofday(&start, NULL);
    Matrix A = rands(M, K);
    printf("Matrix A, shape: %ldx%ld, address in memory:%ld\n", A.height, A.width, (size_t)A.elements);
    printMatrix(A);

    Matrix B = rands(K, N);
    printf("Matrix B, shape: %ldx%ld, address in memory:%ld\n", B.height, B.width, (size_t)B.elements);
    printMatrix(B);

    Matrix C = matMul(A, B);
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);
    gettimeofday(&end, NULL);
    printf("CPU 計算矩陣相乘，用時:");
    long timer1 = getTimer(start, end);
    printTimer(timer1);

    //第二步，將矩陣 A 和 B 拷貝到顯卡中，利用顯卡計算矩陣乘法，再將結果考回 C 中，並計時
    cublasHandle_t handle;
    if(!initCUBLAS(&handle)){
        return EXIT_FAILURE;
    }
    gettimeofday(&start, NULL);
    size_t m = A.height;
    size_t n = B.width;
    size_t k = B.height;
    Matrix devA = {m, k, NULL};
    Matrix devB = {k, n, NULL};
    Matrix devC = {m, n, NULL};
    cudaMalloc(&devA.elements, m*k*sizeof(float));
    cudaMemcpy(devA.elements, A.elements, m*k*sizeof(float), cudaMemcpyHostToDevice);
    cudaMalloc(&devB.elements, k*n*sizeof(float));
    cudaMemcpy(devB.elements, B.elements, k*n*sizeof(float), cudaMemcpyHostToDevice);
    cudaMalloc(&devC.elements, m*n*sizeof(float));
    cudaMemset(devC.elements, 0, m*n*sizeof(float));
    //調用 GPU 進行計算
    size_t blockSize = 32;
    size_t gridWidth = (n + blockSize -1)/blockSize;
    size_t gridHeight = (m + blockSize -1)/blockSize;
    cudaMatMul<<<dim3(gridHeight, gridWidth), dim3(blockSize, blockSize)>>>(devA, devB, devC);
    //從 GPU 取回數據
    deleteMatrix(&C);
    C = zeros(m, n);
    printf("從顯卡取回數據以前的矩陣 C\n");
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);
    cudaMemcpy(C.elements, devC.elements, m*n*sizeof(float), cudaMemcpyDeviceToHost);
    printf("從顯卡取回數據以後的矩陣 C\n");
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);

    gettimeofday(&end, NULL);
    printf("本身寫 CUDA 代碼，利用 GPU 計算矩陣相乘，用時:");
    long timer2 = getTimer(start, end);
    printTimer(timer2);

    //第三步，直接利用 cublas 庫計算矩陣乘法，並計時
    gettimeofday(&start, NULL);
    float scalar = 1.0;
    Matrix devA2 = {m, k, NULL};
    Matrix devB2 = {k, n, NULL};
    Matrix devC2 = {m, n, NULL};
    cudaMalloc(&devA2.elements, m*k*sizeof(float));
    cudaMalloc(&devB2.elements, k*n*sizeof(float));
    cudaMalloc(&devC2.elements, m*n*sizeof(float));

    Matrix C2 = zeros(m, n);

    cublasSetMatrix(m, k, sizeof(float), A.elements, m, devA2.elements, m);
    cublasSetMatrix(k, n, sizeof(float), B.elements, k, devB2.elements, k);
    cublasSetMatrix(m, n, sizeof(float), C2.elements, m, devC2.elements, m);

    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &scalar, devA2.elements, m, devB2.elements, k, &scalar, devC2.elements, m);

    printf("從顯卡取回數據以前的矩陣 C2\n");
    printf("Matrix C2, shape: %ldx%ld, address in memory:%ld\n", C2.height, C2.width, (size_t)C2.elements);
    printMatrix(C2);

    cublasGetMatrix(m, n, sizeof(float), devC2.elements, m, C2.elements, m);

    printf("從顯卡取回數據以後的矩陣 C2\n");
    printf("Matrix C2, shape: %ldx%ld, address in memory:%ld\n", C2.height, C2.width, (size_t)C2.elements);
    printMatrix(C2);
    gettimeofday(&end, NULL);
    long timer3 = getTimer(start, end);
    printf("直接使用 Nvidia 提供的 cuBlas 庫進行矩陣乘法計算，用時：");
    printTimer(timer3);

    printf("\n");
    printf("本身寫 CUDA 代碼利用 GPU 計算矩陣乘法，速度是 CPU 的 %f 倍，利用 cuBlass 庫進行計算，速度是 CPU 的 %f 倍。\n",
            (float)timer1/timer2, (float)timer1/timer3);
    return 0;
}

而後，使用nvcc cublas_test.cu -lcublas -o cublas_test編譯，使用./cublas_test運行，就能夠看到結果了，以下兩圖：

能夠看到，C 語言的三重循環，用時 47.763 秒，而本身寫 CUDA 代碼用顯卡進行計算，用時 222 毫秒，是 CPU 的 214 倍。固然，我本身寫的 CUDA 代碼確定是沒有優化的，若是調用 Nvidia 官方實現的 cuBLAS 庫，則只用時 32 毫秒，是 CPU 的 1489 倍。這計算速度的提高，仍是至關巨大的。

我使用 Python 的 NumPy 測試了一下 OpenBLAS 優化後的性能，以下圖：

用時 6.23 秒，是 C 語言三重循環的 8 倍左右。 Python 的代碼很是簡單，只有四行，以下：

import numpy as np
A = np.random.randn(2000, 2000).astype('float32')
B = np.random.randn(2000, 2000).astype('float32')
%timeit np.dot(A, B)

其中，%timeit是 IPython 或 Jupyter notebook 的魔術指令，專用於測性能。

開發環境的安裝

上面直接給出了代碼和測試結果，那麼在 Linux 下寫 CUDA 程序方便嗎？安裝開發環境容易嗎？個人回答是：至關的方便，至關的容易。

首先，須要一臺帶 Nvidia 顯卡的電腦。其次，須要安裝好 Nvidia 的顯卡驅動，在 Ubuntu 中，這就是一條sudo apt-get install nvidia-384或sudo apt-get install nvidia-390命令的事，我以前的隨筆中有介紹。

最後，就是安裝 CUDA 的開發環境了，也很簡單，一條sudo apt-get install nvidia-cuda-toolkit命令就搞定。而後，就可使用nvcc命令了，nvcc背後調用的是gcc。而 Python 環境，在 Ubuntu 中是天生的方便，前面隨筆中有介紹，這裏就再也不羅嗦。

那用什麼工具寫代碼和調試代碼呢？因爲這裏用的是 C 語言，因此我首先推薦的是 JetBrains 全家桶中的 CLion，以下圖：

只是有點小貴，一年 199 刀。我是不同意使用破解版的，我只能在它提供的 30 天試用期內進行使用。土豪請隨意。

而 Nvidia 也提供了一款基於 Eclipse 的工具 Nsight。安裝完nvidia-cuda-toolkit軟件包後，該工具就有了，不須要另外下載。在命令行輸入nsight命令就能夠啓動它，以下：

而我這裏，就寫一個文件而以，區區兩百多行代碼，要啥 IDE 啊，我用 Vim 足以，以下圖：

關於 Vim 的配置，我前面的隨筆也有介紹。

代碼解釋

首先，我定義了一些輔助的結構和函數。Matrix用來定義一個矩陣，而矩陣的元素我沒有使用二維數組，只是在內存中分配一塊連續空間，以一維數組的形式進行訪問，爲了和 Fortran 兼容，使用列優先存儲，因此定義了宏#define idx(i, j, m) ((j) * (m) + (i))用來根據矩陣的行i和列j定位矩陣的元素在一維數組中的下標。另外，定義的empty(m, n)、zeros(m, n)、rands(m, n)函數分別用來初始化一個空矩陣、零矩陣和隨機矩陣，deleteMatrix()用來刪除一個矩陣並釋放內存。printMatrix()函數用來打印矩陣，便於查看結果。而矩陣相乘的函數matMul()，裏面就是三重循環，很少解釋。

對程序的計時，使用的是 Linux 中的gettimeofday()函數，須要包含<sys/time.h>頭文件。該函數計時的精度能夠到微秒級別，夠用了。同時，定義printTimer()輔助函數用於輸出時間數據。

重點在於對 CUDA 代碼的解釋。

首先，CUDA 把代碼和數據都分爲兩部分，一部分代碼是在 CPU 上執行的，叫 Runtime 函數，一部分代碼是在 GPU 上執行的，叫 Kernel 函數，一部分數據存儲在內存中，咱們稱之爲 Host 上的數據，一部分數據儲存在顯存中，咱們稱之爲 Device 數據。因此，咱們編寫 CUDA 的流程是這樣的，先調用 Runtime 函數初始化 CUDA 環境，最重要的是檢測系統中有沒有 Nvidia 的顯卡，並調用 setDevice()設置使用哪一個顯卡進行計算。好羨慕那些有好幾個顯卡的人啊。而後，在 Host 中分配空間、初始化數據、在 Device 中分配空間，將 Host 中的數據拷貝到 Device 中。這個過程須要用到malloc()、cudaMalloc()、cudaMemcpy()函數，從 Host 中向 Device 中複製數據，須要用到cudaMemcpyHostToDevice常量。

其次，就是設計 Kernel 函數和調用 Kernel 函數了。在這方面，nvcc對 C 語言進行了擴展，提供了__global__關鍵字，用於定義 Kernel 函數。每個 Kernel 函數都在一個 GPU 流處理器上執行，而 GPU 每每有成百上千個流處理器，因此該 Kernel 函數至關於被成百上千個線程同時執行，進行並行計算。定義好 Kernel 函數後，就能夠在 C 程序中調用，使用的是func_name<<<gridDim, blockDim>>>(args)這樣的語法，而gridDim和blockDim能夠指定 GPU 中的線程組織成什麼形狀。多個 thread 能夠組織成一個 block，多個 block 能夠組織成一個 grid。而 grid 和 block 能夠是一維的、二維的、三維的，這樣組織，很是方便咱們設計程序。

最後，就是把顯卡計算後的結果數據再從 Device 中拷貝到 Host 中，仍是用的cudaMemcpy函數，只不過用的是cudaMemcpyDeviceToHost常數。

若是把 GPU 中的線程組織成二維的，就是這樣：

在本例中，咱們須要計算一個 M×K 的矩陣和一個 K×N 的矩陣相乘，結果是一個 M×N 的矩陣，咱們就能夠啓用 M×N 個線程，並把它組織成如上形狀，每個線程只計算結果矩陣中的一個元素，以下圖：

CUDA 中每一個 block 中的 thread 數是有上限的，在個人電腦上，該上限是 1024。因此我定義blockSize = 32，而後 block 的形狀就是(blockSize, blockSize)，也就是 block 的大小爲 32×32。而後再把 grid 的形狀定義爲((m + blockSize -1)/blockSize, (n + blockSize -1)/blockSize)，這樣，這個 grid 中的總線程數就是 M×N 了，並且線程排列的形狀和結果矩陣的形狀徹底同樣，因此每一個線程計算結果矩陣中的一個元素，極其方便。

下面來看看完整的流程。

先在 Host 上創建矩陣 A、B、C，並初始化：

Matrix A = rands(M, K);
    printf("Matrix A, shape: %ldx%ld, address in memory:%ld\n", A.height, A.width, (size_t)A.elements);
    printMatrix(A);

    Matrix B = rands(K, N);
    printf("Matrix B, shape: %ldx%ld, address in memory:%ld\n", B.height, B.width, (size_t)B.elements);
    printMatrix(B);

    Matrix C = matMul(A, B);
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);

而後建立三個變量用來保存在 Device 上的矩陣的形狀和數據指針 devA、devB、devC，這三個變量剛開始其中的數據爲空，以下：

size_t m = A.height;
    size_t n = B.width;
    size_t k = B.height;
    Matrix devA = {m, k, NULL};
    Matrix devB = {k, n, NULL};
    Matrix devC = {m, n, NULL};

而後在 Device 中分配空間，並把 Host 中的數據拷貝到 Device 中，這時，devA、devB、devC中的數據指針指向的是 Device 中的空間，以下：

cudaMalloc(&devA.elements, m*k*sizeof(float));
    cudaMemcpy(devA.elements, A.elements, m*k*sizeof(float), cudaMemcpyHostToDevice);
    cudaMalloc(&devB.elements, k*n*sizeof(float));
    cudaMemcpy(devB.elements, B.elements, k*n*sizeof(float), cudaMemcpyHostToDevice);
    cudaMalloc(&devC.elements, m*n*sizeof(float));
    cudaMemset(devC.elements, 0, m*n*sizeof(float));

再而後，調用 Kernel 函數cudaMatMul()，指定啓用多少個線程，並指定 grid 和 block 的形狀，並把 devA、devB、devC 傳遞給它作參數。以下：

size_t blockSize = 32;
    size_t gridWidth = (n + blockSize -1)/blockSize;
    size_t gridHeight = (m + blockSize -1)/blockSize;
    cudaMatMul<<<dim3(gridHeight, gridWidth), dim3(blockSize, blockSize)>>>(devA, devB, devC);

Kernel 函數cudaMatMul()是怎麼定義的呢？以下：

__global__ void cudaMatMul(Matrix devA, Matrix devB, Matrix devC){
  size_t m = devA.height;
  size_t k = devB.height;
  size_t n = devB.width;

  size_t j = blockIdx.y * blockDim.y + threadIdx.y;
  size_t i = blockIdx.x * blockDim.x + threadIdx.x;

  if(i<m && j<n){
    float sum = 0.0;
    for(size_t p=0; p<k; p++){
      sum += devA.elements[idx(i, p, m)] * devB.elements[idx(p, j, k)];
    }
    devC.elements[idx(i, j, m)] = sum;
  }
}

能夠看到，定義 Kernel 函數時，用到了__global__關鍵字，而且使用了threadIdx.x、threadIdx.y、blockIdx.x、blockIdx.y、blockDim.x和blockDim.y這些內建變量，根據前面的介紹，咱們很是容易猜到這些變量的意義。而後，很是容易算出當前線程計算的是結果矩陣的哪一個元素 C[i, j]，再而後，三重循環變成一重循環，計算起來不要過輕鬆，效率主要取決於訪問內存的速度。

最後，將結果矩陣從 Device 拷貝到 Host 並顯示，都是常規操做。以下：

deleteMatrix(&C);
    C = zeros(m, n);
    printf("從顯卡取回數據以前的矩陣 C\n");
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);
    cudaMemcpy(C.elements, devC.elements, m*n*sizeof(float), cudaMemcpyDeviceToHost);
    printf("從顯卡取回數據以後的矩陣 C\n");
    printf("Matrix C, shape: %ldx%ld, address in memory:%ld\n", C.height, C.width, (size_t)C.elements);
    printMatrix(C);

若是不想本身寫 CUDA 代碼，可使用 Nvidia 提供的 cuBLAS 庫。除了 cuBLAS 庫以外，還有 cuSPARSE 用來計算稀疏矩陣、cuFFT 用來進行傅立葉變化、cuSOLVER 用來解線性方程組、cuDNN 用來加速神經網絡的計算等等。這些庫優化得好，速度比咱們本身寫 CUDA 代碼要快一些。另外，若是直接用 cuBLAS，則只須要用到 Runtime 函數，不須要本身定義 Kernel 函數，因此文件的擴展名能夠不是.cu，編譯器能夠不用nvcc而直接用gcc，只須要包含正確的頭文件和鏈接正確的庫就能夠了。固然，使用nvcc也能夠。

cuBLAS 使用<cublas_v2.h>頭文件，編譯時使用-lcublas鏈接相應的庫。爲何頭文件中有個v2呢？這一由於這個新版本是線程安全的，因此每一個 cuBLAS 函數都須要一個handle做爲參數，這個handle可使用cublasCreate()函數建立，不一樣的線程可使用不一樣的handle。我把建立handle()的工做也放到 CUDA 初始化的代碼裏面了，在initCUBLAS()函數中。

分配內存和拷貝數據的過程是同樣的，在 cuBLAS 中，同樣使用cudaMalloc()和cudaMemcpy()函數，可是推薦使用cublasSetMatrix()和cublasGetMatrix()函數，由於這兩個函數能夠拷貝子矩陣的數據。

而後就是調用cublasSgemm()函數進行計算了，爲何這個函數名這麼奇怪呢？這是由於 BLAS 庫中的函數名原本就極度簡化，好比mm就表明兩個矩陣相乘，mv就表明一個矩陣和一個向量相乘。ge可能指的就是普通矩陣，除此以外，還有帶狀矩陣、對稱矩陣等，分別用b、s指示。而中間的大寫S，指的是數據類型，S是單精度浮點數，D是雙精度浮點數，C爲單精度複數，Z爲雙精度複數。

具體的函數，你們直接閱讀 Nvidia 的官方文檔吧。

總結

使用顯卡進行科學計算，得益於 GPU 大量的流處理器，計算速度能夠獲得成百上千倍的提高。在 Linux 環境下，安裝 Nvidia 的顯卡驅動和 CUDA toolkit，也是極其方便的事，並且 Nvidia 提供完善的文檔和庫，真的是咱們的福音。

求打賞

我對此次寫的這個系列要求是很是高的：首先內容要有意義、夠充實，信息量要足夠豐富；其次是每個知識點要講透徹，不能模棱兩可含糊不清；最後是包含豐富的截圖，讓那些不想裝 Linux 系統的朋友們也能夠領略到 Linux 桌面的風采。若是個人努力獲得你們的承認，能夠掃下面的二維碼打賞一下：

版權申明

該隨筆由京山遊俠在2019年01月14日發佈於博客園，引用請註明出處，轉載或出版請聯繫博主。QQ郵箱：1841079@qq.com