[平行運算] 用CUDA實現矩陣相乘

CUDA架構在device方面儲存資料的方式有兩種, 分別是linear memory和CUDA arrays, 在這裡我們使用linear memory, 因此在device的時候, 我們必須把二維矩陣轉為一維矩陣來看待, 其實也蠻簡單的, 官網上的文件是用row-major, 意思是以row為主, 一次把一個row依序塞進一維陣列, 所以要access的時候就會是 M(row, col) = *(陣列初始位址 + row * M.width + col); 也有人是用column major, 其實就是相反過來, M(row, col) = *(陣列初始位址 + row+ col*M.height);

在這裡其實有一個還沒提到的概念, 被稱作shared memory, 這是CUDA架構下的一種記憶體架構, 可以提升計算速度, 部過這個概念要等到下一篇文章我才會寫, 因此下面的程式碼是沒有shared memory的版本!

此外, 這裡的陣列大小選取都是blocksize的倍數, 因為若不是blocksize的倍數, 記憶體在分配上會有技術問題要處理, 此處也先省略

#include <iostream>
#include <cstdlib> 
#include <ctime>

#include <cuda_runtime.h>

#include "device_launch_parameters.h"  
#define BLOCK_SIZE 2
#define A_row 2
#define A_col 3
#define B_row A_col
#define B_col 2
#define C_row A_row
#define C_col B_col

using namespace std;

// 宣告Matrix class
class Matrix{
    
    public:
        int width;
        int height;
        float* elements;
};

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)
{
    float ans = 0;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    for (int p = 0; p < A.width; ++p)
        ans += A.elements[row * A.width + p]
                * B.elements[p * B.width + col];
    C.elements[row * C.width + col] = ans;
}


int main(void)
{
    float A[A_row * BLOCK_SIZE][A_col * BLOCK_SIZE];
    float B[B_row * BLOCK_SIZE][B_col * BLOCK_SIZE];
    float C[A_row * BLOCK_SIZE][B_col * BLOCK_SIZE]={0};
    
    //隨機分配值給矩陣
    for(int i=0;i< A_row * BLOCK_SIZE;i++)
    {
        for(int j=0;j< A_col * BLOCK_SIZE;j++)
        {
            A[i][j]=rand();        
        }
        
    }
    
    for(int i=0;i< B_col * BLOCK_SIZE;i++)
    {
        for(int j=0;j< B_row * BLOCK_SIZE;j++)
        {
            B[i][j]=rand();
        }
    }
    
    Matrix h_A,h_B,h_C;
    
    //把給定的矩陣轉為Matrix類別
     
    h_A.height = A_row * BLOCK_SIZE;
    h_B.height = B_row * BLOCK_SIZE;
    h_C.height = C_row * BLOCK_SIZE;
     
    h_A.width = A_col * BLOCK_SIZE;
    h_B.width = B_col * BLOCK_SIZE;
    h_C.width = C_col * BLOCK_SIZE;

    h_A.elements = (float*)malloc(A_row*A_col*BLOCK_SIZE*BLOCK_SIZE*sizeof(float));

    h_B.elements = (float*)malloc(B_row*B_col*BLOCK_SIZE*BLOCK_SIZE*sizeof(float));  

    h_C.elements = (float*)malloc(C_row*C_col*BLOCK_SIZE*BLOCK_SIZE*sizeof(float));     

    for(int i=0;i< A_row * BLOCK_SIZE;i++)
    {
        for(int j=0;j< A_col * BLOCK_SIZE;j++)
        {
            *(h_A.elements + i*h_A.width + j) = A[i][j];
        }
        
    }


    for(int i=0;i< B_row * BLOCK_SIZE;i++)
    {
        for(int j=0;j< B_col * BLOCK_SIZE;j++)
        {
            *(h_B.elements + i*h_B.width + j) = B[i][j];
        }
    }
    
    // 把h_A 和 h_B 複製給 device
    Matrix d_A;
    d_A.width = h_A.width; d_A.height = h_A.height;
    size_t size = h_A.width * h_A.height * sizeof(float);
    cudaMalloc(&d_A.elements, size);
    cudaMemcpy(d_A.elements, h_A.elements, size, cudaMemcpyHostToDevice);
    
    Matrix d_B;
    d_B.width = h_B.width; d_B.height = h_B.height;
    size = h_B.width * h_B.height * sizeof(float);
    cudaMalloc(&d_B.elements, size);
    cudaMemcpy(d_B.elements, h_B.elements, size, cudaMemcpyHostToDevice);

    // 在 device memory 裡創出 d_C 的記憶體空間
    Matrix d_C;
    d_C.width = h_C.width; d_C.height = h_C.height;
    size = h_C.width * h_C.height * sizeof(float);
    cudaMalloc(&d_C.elements, size);

    // 實際計算
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid(h_B.width / dimBlock.x, h_A.height / dimBlock.y);
    MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
    
    // 把算好的結果從d_C 複製回 h_C
    cudaMemcpy(h_C.elements, d_C.elements, size,
               cudaMemcpyDeviceToHost);
    
    // 驗算
    int ok=1;
    
    for(int i=0;ok && i< A_row * BLOCK_SIZE ;i++)
    {
        for(int j=0;ok && j< B_col * BLOCK_SIZE;j++)
        {
            for(int k=0;k< A_col * BLOCK_SIZE;k++)
            {
                C[i][j]+=A[i][k]*B[k][j];
            }
            if(C[i][j]!=*(h_C.elements+i*h_C.width+j))
            {
                cout<<"fail";
                ok=0;
                break;
            }
        }
    }
    
    if(ok==1)
    {
        cout<<"success"<<endl;
    }

    //釋放記憶體空間
    cudaFree(d_A.elements);
    cudaFree(d_B.elements);
    cudaFree(d_C.elements);
   

    system("pause");
}

reference:
1. Nvidia developer zone