特別說明：要在我的隨筆後寫評論的小夥伴們請註意了，我的博客開啟了 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郵箱：[email protected]

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