第二個 CUDA 程式

洛欣發表於2009-11-17

前面介紹的計算平方和的程式，似乎沒有什麼實用價值。所以我們的第二個 CUDA 程式，要做一個確實有（某些）實用價值的程式，也就是進行矩陣乘法。而且，這次我們會使用浮點數。

雖然矩陣乘法有點老套，不過因為它相當簡單，而且也可以用來介紹一些有關 CUDA 的有趣性質。

矩陣乘法

為了單純起見，我們這裡以方形的矩陣為例子。基本上，假設有兩個矩陣 A 和 B，則計算 AB = C 的方法如下：

    for(i = 0; i < n; i++) {
        for(j = 0; j < n; j++) {
            C[i][j] = 0;
            for(k = 0; k < n; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

一開始，我們先準備好產生資料、設定 CUDA 等等的工作：

    int main()
    {
        float *a, *b, *c, *d;
        int n = 1000;
        if(!InitCUDA()) return 0;
        a = (float*) malloc(sizeof(float) * n * n);
        b = (float*) malloc(sizeof(float) * n * n);
        c = (float*) malloc(sizeof(float) * n * n);
        d = (float*) malloc(sizeof(float) * n * n);
        srand(0);
        matgen(a, n, n);
        matgen(b, n, n);
        clock_t time = matmultCUDA(a, n, b, n, c, n, n);
        matmult(a, n, b, n, d, n, n);
        compare_mat(c, n, d, n, n);
        double sec = (double) time / CLOCKS_PER_SEC;
        printf("Time used: %.2f (%.2lf GFLOPS)\n", sec,
            2.0 * n * n * n / (sec * 1E9));
        return 0;
    }

InitCUDA 函式和第一個 CUDA 程式一樣，可以直接參考前面的文章。以下是上面用到的一些其它的函式：

產生矩陣：

    void matgen(float* a, int lda, int n)
    {
        int i, j;
        for(i = 0; i < n; i++) {
            for(j = 0; j < n; j++) {
                a[i * lda + j] = (float) rand() / RAND_MAX +
                    (float) rand() / (RAND_MAX * RAND_MAX);
            }
        }
    }

這個函式只是利用亂數產生器把矩陣填滿 0 ~ 1 之間的數字。特別注意到因為 C 語言中無法宣告變動大小的二維矩陣，所以我們使用 i * lda + j 的方式。

進行矩陣乘法：

    void matmult(const float* a, int lda, const float* b, int ldb,
        float* c, int ldc, int n)
    {
        int i, j, k;
        for(i = 0; i < n; i++) {
            for(j = 0; j < n; j++) {
                double t = 0;
                for(k = 0; k < n; k++) {
                    t += a[i * lda + k] * b[k * ldb + j];
                }
                c[i * ldc + j] = t;
            }
        }
    }

這是以 CPU 進行矩陣乘法、用來進行驗證答案正確與否的程式。特別注意到它用 double 來儲存暫時的計算結果，以提高精確度。

驗證結果：

    void compare_mat(const float* a, int lda,
        const float* b, int ldb, int n)
    {
        float max_err = 0;
        float average_err = 0;
        int i, j;
        for(i = 0; i < n; i++) {
            for(j = 0; j < n; j++) {
                if(b[i * ldb + j] != 0) {
                    float err = fabs((a[i * lda + j] -
                        b[i * ldb + j]) / b[i * ldb + j]);
                    if(max_err < err) max_err = err;
                    average_err += err;
                }
            }
        }
        printf("Max error: %g Average error: %g\n",
            max_err, average_err / (n * n));
    }

這個函式計算兩個矩陣的最大相對誤差和平均相對誤差，並把結果印出來。

最後是 CUDA 的矩陣乘法的部份：

    #define NUM_THREADS 256
    clock_t matmultCUDA(const float* a, int lda,
        const float* b, int ldb, float* c, int ldc, int n)
    {
        float *ac, *bc, *cc;
        clock_t start, end;
        start = clock();
        cudaMalloc((void**) &ac, sizeof(float) * n * n);
        cudaMalloc((void**) &bc, sizeof(float) * n * n);
        cudaMalloc((void**) &cc, sizeof(float) * n * n);
        cudaMemcpy2D(ac, sizeof(float) * n, a, sizeof(float) * lda,
            sizeof(float) * n, n, cudaMemcpyHostToDevice);
        cudaMemcpy2D(bc, sizeof(float) * n, b, sizeof(float) * ldb,
            sizeof(float) * n, n, cudaMemcpyHostToDevice);
        int blocks = (n + NUM_THREADS - 1) / NUM_THREADS;
        matMultCUDA<<>>
            (ac, n, bc, n, cc, n, n);
        cudaMemcpy2D(c, sizeof(float) * ldc, cc, sizeof(float) * n,
        sizeof(float) * n, n, cudaMemcpyDeviceToHost);
        cudaFree(ac);
        cudaFree(bc);
        cudaFree(cc);
        end = clock();
        return end - start;
    }

這個函式相當單純，就是在顯示記憶體中配置存放矩陣的記憶體，然後把主記憶體中的矩陣資料複製到顯示記憶體上。不過，因為我們的矩陣乘法函式可以指定 pitch（即 lda、ldb、和 ldc），所以如果用一般的 cudaMemcpy 函式來複製記憶體的話，會需要每個 row 都分開複製，那會需要呼叫很多次 cudaMemcpy 函式，會使效率變得很差。因此，在這裡我們用了一個新的 cudaMemcpy2D 函式，它是用來複製二維陣列，可以指定陣列的 pitch。這樣就可以透過一次函式呼叫就可以了。

進行計算的 kernel 如下：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        const int tid = threadIdx.x;
        const int bid = blockIdx.x;
        const int idx = bid * blockDim.x + tid;
        const int row = idx / n;
        const int column = idx % n;
        int i;
        if(row < n && column < n) {
            float t = 0;
            for(i = 0; i < n; i++) {
                t += a[row * lda + i] * b[i * ldb + column];
            }
            c[row * ldc + column] = t;
        }
    }

這個函式一開始先從 bid 和 tid 計算出這個 thread 應該計算的 row 和 column，在判斷 row 和 column 在範圍內之後，就直接進行計算，並把結果寫到 c 矩陣中，是非常單純的函式。

在 GeForce 8800GT 上實際執行的結果如下：

Max error: 2.01484e-006 Average error: 3.36637e-007
Time used: 1.1560 (1.73 GFLOPS)

可以看到兩個問題：

很明顯的，執行效率相當低落。
最大相對誤差偏高。理想上應該要低於 1e-6。

計算結果的誤差偏高的原因是，在 CPU 上進行計算時，我們使用 double（即 64 bits 浮點數）來累進計算過程，而在 GPU 上則只能用 float（32 bits 浮點數）。在累加大量數字的時候，由於累加結果很快會變大，因此後面的數字很容易被捨去過多的位數。

由於 CUDA 的浮點數運算，在進行加、減、乘法時是符合 IEEE 754 規定的精確度的，因此，我們可以利用 Kahan's Summation Formula 來提高精確度。把程式改成：

    if(row < n && column < n) {
        float t = 0;
        float y = 0;
        for(i = 0; i < n; i++) {
            float r;
            y -= a[row * lda + i] * b[i * ldb + column];
            r = t - y;
            y = (r - t) + y;
            t = r;
        }
    }

修改後的程式的執行結果是：

Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 1.1560 (1.73 GFLOPS)

可以看到相對誤差有很大的改善，效率則沒什麼變化。

由於 Kahan's Summation Formula 需要的運算量提高，但是效率卻沒有什麼改變，可以看出這個 kernel 主要的瓶頸應該是在記憶體的存取動作上。這是因為有大量的記憶體讀取是重覆的。例如，矩陣 a 的一個 row 在每次進行計算時都被重覆讀入，但這是相當浪費的。這樣的計算方式，總共需要讀取 2*n³ 次記憶體。如果讓一個 row 只需要讀入一次的話，就可以減到為 n³+n² 次。

第一個改良

和我們的第一個 CUDA 程式一樣，我們可以利用 shared memory 來儲存每個 row 的資料。不過，因為只有同一個 block 的 threads 可以共用 shared memory，因此現在一個 row 只能由同一個 block 的 threads 來進行計算。另外我們也需要能存放一整個 row 的 shared memory。因此，把先把呼叫 kernel 的部份改成：

matMultCUDA<<>>
(ac, n, bc, n, cc, n, n);

kernel 的部份則改成：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        extern __shared__ float data[];
        const int tid = threadIdx.x;
        const int row = blockIdx.x;
        int i, j;
        for(i = tid; i < n; i += blockDim.x) {
            data[i] = a[row * lda + i];
        }
        __syncthreads();
        for(j = tid; j < n; j += blockDim.x) {
            float t = 0;
            float y = 0;
            for(i = 0; i < n; i++) {
                float r;
                y -= data[i] * b[i * ldb + j];
                r = t - y;
                y = (r - t) + y;
                t = r;
            }
            c[row * ldc + j] = t;
        }
    }

第一個部份先把整個 row 讀到 shared memory 中，而第二個部份則進行計算，並沒有太大的變化。主要的差別是現在一個 row 只由一個 block 進行計算。

在 GeForce 8800GT 上，執行的結果是：

Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.4220 (4.74 GFLOPS)

很明顯的，計算的結果並沒有改變，不過速度則提高了超過一倍。雖然如此，但是這樣的效率仍不盡理想，因為理論上 GeForce 8800GT 有超過 300GFLOPS 的運算效能。即使是把 Kahan's Summation Formula 所需要的額外運算考慮進去，這樣的效率仍然連理論最大值的十分之一都不到。

會有這樣的結果，原因其實還是同樣的：對記憶體的存取次數太多了。雖然現在 A 矩陣的 row 的資料已經不再需要重覆讀取，但是 B 矩陣的 column 的資料仍然一直被重覆讀取。

另一個問題比較不是那麼明顯：對 B 矩陣的讀取，雖然看起來不連續，但實際上它是連續的。這是因為不同的 thread 會讀取不同的 column，因此同時間每個 thread 讀取的各個 column 加起來，就是一個連續的記憶體區塊。那麼，為什麼效率還是不佳呢？這是因為，GPU 上的記憶體控制器，從某個固定的倍數位址開始讀取，才會有最高的效率（例如 16 bytes 的倍數）。由於矩陣大小並不是 16 的倍數（這裡使用的是 1000x1000 的矩陣），所以造成效率不佳的情形。

要解決這個問題，我們可以在 cudaMalloc 的時候稍微修改一下，讓寬度變成適當的倍數就可以了。但是，適當的倍數是多少呢？幸運的是，我們並不需要知道這些細節。CUDA 提供了一個 cudaMallocPitch 的函式，可以自動以最佳的倍數來配置記憶體。因此，我們可以把 cudaMalloc 的部份改成：

    size_t pitch_a, pitch_b, pitch_c;
    cudaMallocPitch((void**) &ac, &pitch_a, sizeof(float) * n, n);
    cudaMallocPitch((void**) &bc, &pitch_b, sizeof(float) * n, n);
    cudaMallocPitch((void**) &cc, &pitch_c, sizeof(float) * n, n);

cudaMallocPitch 函式會以適當的倍數配置記憶體，並把配置的寬度傳回。因此，在把矩陣複製到顯示記憶體上時，要使用它傳回的寬度：

    cudaMemcpy2D(ac, pitch_a, a, sizeof(float) * lda,
        sizeof(float) * n, n, cudaMemcpyHostToDevice);
    cudaMemcpy2D(bc, pitch_b, b, sizeof(float) * ldb,
        sizeof(float) * n, n, cudaMemcpyHostToDevice);

呼叫 kernel 的部份也需要修改：

    matMultCUDA<<>>
        (ac, pitch_a / sizeof(float), bc, pitch_b / sizeof(float),
        cc, pitch_c / sizeof(float), n);

同樣的，把計算結果複製回到主記憶體時，也要使用傳回的寬度值：

cudaMemcpy2D(c, sizeof(float) * ldc, cc, pitch_c,
sizeof(float) * n, n, cudaMemcpyDeviceToHost);

這樣就完成了。Kernel 部份則不需要修改。

這樣的修改有多大的效果呢？在 GeForce 8800GT 上執行，結果如下：

Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.1250 (16.00 GFLOPS)

可以看到，執行速度又再大幅提高了三倍多！而這只是把記憶體的配置方式稍微修改一下而已。

雖然執行速度提高了很多，但是，和前面提到的理論值相比，其實還是有相當的差距。這是因為，前面也提到過，這樣的做法需要 n³+n²次的記憶體讀取，和 n²次的記憶體寫入動作。由於 n = 1000，每個數字的大小是 32 bits，所以總共的記憶體存取資料量約為 4GB。除以實際執行的時間 0.125 秒，得到的頻寬數值是約 32GB/s，這已經接近 GeForce 8800GT 顯示記憶體的頻寬了。由於我們計算時間的時候，把配置記憶體、以及資料的複製動作也計算進去，因此實際上花費在 kernel 的時間是更短的（約 0.09 秒）。因此，可以很明顯的看出，這個程式的效率，是受限於記憶體頻寬的。

進一步的改良

上一節的結論顯示出，矩陣乘法的程式，效率是受限於記憶體頻寬的。那有沒有辦法降低記憶體的存取次數呢？答案當然是有的，不然就不會有這一節了 :)

要進一步降低記憶體頻寬的使用，可以注意到，在上一節的方法中，雖然 A 矩陣的存取次數被減至最低，但是 B 矩陣的存取次數並沒有減少。這是因為我們只將 A 矩陣的 row 載入到 shared memory 中，但是 B 矩陣的 column 也是有被重覆使用的。理想上應該也可以避免重覆載入才對。不過，由於 B 矩陣的 column 使用的時機，和 A 矩陣的 row 是不同的，所以並不能直接這樣做。

解決方法是 "blocking"。也就是把整個矩陣乘法的動作，切割成很多小矩陣的乘法。例如，要計算 C 矩陣的 (0, 0) ~ (15, 15) 的值，可以把它想成：

A(0~15, 0~15) * B(0~15, 0~15) + A(0~15,16~31) * B(16~31, 0~15)
+ A(0~15, 32~47) * B(32~47, 0~15) + ...

這樣一來，我們就可以把兩個小矩陣載入到 shared memory，則小矩陣本身的乘法就不需要再存取任何外部的記憶體了！這樣一來，假設小矩陣的大小是 k，則實際上需要的記憶體存取次數就會變成約 2k²(n/k)³ = 2n³/k。

由於目前 CUDA 每個 block 的 thread 數目最多是 512，因此 k = 16 似乎是一個相當理想的數字（共 256 個 threads）。因此，對於一個 n = 1000 的矩陣來說，我們可以把記憶體存取的量減少到約 500MB，也就是上一節的存取量的 1/8。理論上，這樣應該可以讓效率提高八倍（假設沒有遇到別的瓶頸）。

為了方便進行區塊的計算，我們讓每個 block 有 16x16 個 threads，再建立 (n/16)x(n/16) 個 blocks。把呼叫 kernel 的地方改成：

    int bx = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    dim3 blocks(bx, bx);
    dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
    matMultCUDA<<>>(ac, pitch_a / sizeof(float),
        bc, pitch_b / sizeof(float), cc, pitch_c / sizeof(float), n);

BLOCK_SIZE 則是定義成 16。dim3 是 CUDA 的一種資料型態，表示一個 3D 的向量。在這裡，我們透過 dim3 來建立 16x16 個 threads 的 block，和 (n/16)x(n/16) 個 blocks。

Kernel 程式的部份，則改成：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        __shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
        const int tidc = threadIdx.x;
        const int tidr = threadIdx.y;
        const int bidc = blockIdx.x * BLOCK_SIZE;
        const int bidr = blockIdx.y * BLOCK_SIZE;
        int i, j;
        float results = 0;
        float comp = 0;
        for(j = 0; j < n; j += BLOCK_SIZE) {
            if(tidr + bidr < n && tidc + j < n) {
                matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];
            }
            else {
                matA[tidr][tidc] = 0;
            }
            if(tidr + j < n && tidc + bidc < n) {
                matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc];
            }
            else {
                matB[tidr][tidc] = 0;
            }
            __syncthreads();
            for(i = 0; i < BLOCK_SIZE; i++) {
                float t;
                comp -= matA[tidr][i] * matB[i][tidc];
                t = results - comp;
                comp = (t - results) + comp;
                results = t;
            }
            __syncthreads();
        }
        if(tidr + bidr < n && tidc + bidc < n) {
            c[(tidr + bidr) * ldc + tidc + bidc] = results;
        }
    }

注意到因為我們現在使用 16x16 的 threads，因此 threadIdx 變數可以取得 threadIdx.x 和 threadIdx.y，範圍分別是 0 ~ 15。blockIdx.x 和 blockIdx.y 變數也是同樣的情形，範圍分別是 0 ~ n/16。

在程式中，因為矩陣的大小不一定會是 16 的倍數，因此需要使用 if 判斷式檢查是否超出矩陣範圍。

這個版本在 GeForce 8800GT 上的執行結果如下：

Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.0780 (25.64 GFLOPS)

速度雖然提高了，但是似乎並沒有達到預期中的八倍。當然，前面提到過，我們在計算時間時，把一些複製記憶體、配置記憶體的動作也計算在內，這些動作的時間並不會縮短。實際上 kernel 的執行時間，大約是 0.053 秒左右（約略相當於 38GFLOPS），比上一節的版本快了將近一倍。

如果這一版程式已經不再限於記憶體頻寬，那為什麼沒有達到預期的效率呢？這是因為這一版程式已經是限於運算速度了。除了使用 Kahan's Summation Formula 會需要更多的運算之外，程式中也有大量計算矩陣位址的乘法等等，這都會需要花費運算資源。另外，那些用來判斷超出矩陣範圍的 if 判斷式，也會有一定的影響。

要把那些 if 判斷式去掉，有一個方法是，在配置記憶體時，就配置成 16 的倍數，並在複製矩陣到顯示記憶體之前，先將它清為 0。如下所示：

    int newn = ((n + BLOCK_SIZE - 1) / BLOCK_SIZE) * BLOCK_SIZE;
    cudaMallocPitch((void**) &ac, &pitch_a,
        sizeof(float) * newn, newn);
    cudaMallocPitch((void**) &bc, &pitch_b,
        sizeof(float) * newn, newn);
    cudaMallocPitch((void**) &cc, &pitch_c,
        sizeof(float) * newn, newn);
    cudaMemset(ac, 0, pitch_a * newn);
    cudaMemset(bc, 0, pitch_b * newn);

這樣一來，我們就可以把 kernel 中的 if 判斷式都移除了：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        __shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
        const int tidc = threadIdx.x;
        const int tidr = threadIdx.y;
        const int bidc = blockIdx.x * BLOCK_SIZE;
        const int bidr = blockIdx.y * BLOCK_SIZE;
        int i, j;
        float results = 0;
        float comp = 0;
        for(j = 0; j < n; j += BLOCK_SIZE) {
            matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];
            matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc];
            __syncthreads();
            for(i = 0; i < BLOCK_SIZE; i++) {
                float t;
                comp -= matA[tidr][i] * matB[i][tidc];
                t = results - comp;
                comp = (t - results) + comp;
                results = t;
            }
            __syncthreads();
        }
        c[(tidr + bidr) * ldc + tidc + bidc] = results;
    }

這個版本的執行結果是：

Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.0780 (25.64 GFLOPS)

似乎沒有改善。不過，實際上 kernel 的執行時間已經減少到 0.042 秒（約略相當於 48GFLOPS）。

結論

有些讀者可能會想，如果把 block 再變得更大（例如 32x32）是否會有幫助呢？當然，由於最後的程式已經不再是受限於記憶體頻寬（在 0.042 秒記憶體取 500MB 的資料約相當於 12GB/s 的頻寬），所以把 block 再加大並不會有幫助了。而且，由於一個 block 內的 thread 數目最多隻能到 512 個，將 block 變大也會造成很多額外負擔。而且 shared memory 的大小也有限制（GeForce 8800GT 的 shared memory 大小限制是 16384 bytes），所以也不能任意增加 block 的大小。

最後一版程式的完整檔案可以從這裡下載。

原文

來自 “ ITPUB部落格 ” ，連結：http://blog.itpub.net/22785983/viewspace-619830/，如需轉載，請註明出處，否則將追究法律責任。

CUDA C 程式設計權威指南學習筆記：第二章 CUDA程式設計模型
2018-03-02
程式設計筆記模型
第二篇：CUDA 並行程式設計簡介
2017-01-19
並行行程程式設計
CUDA程式設計模式
2010-01-25
程式設計設計模式
cuda程式設計與gpu平行計算（四）：cuda程式設計模型
2020-12-31
程式設計GPU模型
CUDA學習筆記-1: CUDA程式設計概覽
2021-08-04
筆記程式設計
【CUDA學習】核心程式除錯
2013-07-25
除錯
cuda
2024-03-10
CUDA 8的混合精度程式設計
2020-12-29
程式設計
CUDA 學習筆記之程式棧
2010-04-30
筆記
CUDA程式優化心得之序列優化
2010-05-20
優化
cuda 流
2024-03-31
第二個OpenGL程式，矩形（VAO VBO）_後續 EBO
2024-04-04
ubuntu 14.04 安裝cuda 7.5/CUDA 8.0
2016-04-14
Ubuntu
CUDA程式設計模型【中科院課件】
2017-11-08
程式設計模型
CUDA程式優化心得之錯誤處理
2010-05-14
優化
自己用CUDA做的簡單的加密程式
2009-11-13
加密
cuda和cudatoolkit
2020-11-28
CUDA優化
2010-05-28
優化
CUDA簡介
2009-12-28
CUDA架構
2010-01-05
架構
GPU高效能程式設計CUDA實戰(二)
2019-04-01
GPU程式設計
第三篇：CUDA 標準程式設計模式
2017-01-19
程式設計設計模式
第五篇：CUDA 並行程式中的同步
2017-01-19
並行行程
CUDA程式設計（4.1）—— 宣告符（global、device、host等）
2024-03-23
程式設計dev
CUDA進階第三篇：CUDA計時方式
2018-03-06
第二章 - 程式
2023-03-19
CUDA 開發包安裝環境搭建程式開發
2009-11-16
使用Visual Studio 2005 撰寫CUDA 程式
2009-11-16
CUDA學習指南
2020-10-04
平行計算cuda
2017-11-29
CUDA執行模式
2010-03-01
模式
CUDA C語言
2010-01-10
C語言
cuda 核函式
2024-03-09
函式
cuda安裝教程
2024-04-18
CUDA的問題
2024-04-19
5.第二個Activity
2020-12-10
[譯]在CUDA C/C++中如何衡量程式碼效能
2019-05-12
C++
簡述CUDA執行緒及求CUDA中執行緒索引
2022-10-13
執行緒索引

第二個 CUDA 程式

相關文章