第二個 CUDA 程式
 

前面介紹的計算平方和的程式,似乎沒有什麼實用價值。所以我們的第二個 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<<<blocks * n, NUM_THREADS>>>
            (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)

可以看到兩個問題:

  1. 很明顯的,執行效率相當低落。
  2. 最大相對誤差偏高。理想上應該要低於 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*n3 次記憶體。如果讓一個 row 只需要讀入一次的話,就可以減到為 n3+n2 次。

第一個改良

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

        matMultCUDA<<<n, NUM_THREADS, sizeof(float) * n>>>
            (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<<<n, NUM_THREADS, sizeof(float) * n>>>
        (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)

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

雖然執行速度提高了很多,但是,和前面提到的理論值相比,其實還是有相當的差距。這是因為,前面也提到過,這樣的做法需要 n3+n2 次的記憶體讀取,和 n2 次的記憶體寫入動作。由於 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,則實際上需要的記憶體存取次數就會變成約 2k2(n/k)3 = 2n3/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<<<blocks, threads>>>(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.xthreadIdx.y,範圍分別是 0 ~ 15。blockIdx.xblockIdx.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 的大小。

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