矩陣乘法是有實用價值的程式，我們會使用浮點數。

雖然矩陣乘法有點老套，不過因為它相當簡單，而且也可以用來介紹一些有關 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)

　　可以看到兩個問題：

　　很明顯的，執行效率相當低落。

　　最大相對誤差偏高。理想上應該要低於 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<<>>

　　(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)

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

　　雖然執行速度提高了很多，但是，和前面提到的理論值相比，其實還是有相當的差距。這是因為，前面也提到過，這樣的做法需要 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.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 的大小。