# OpenCL 學習筆記 （三）： Local Memory 之特性

###### tags: `OpenCL`

接下來我們要講解如何利用 local memory 設計一個有效率的 kernel。首先我們先展示基本的程式碼，請到[這裡](https://drive.google.com/drive/folders/1JsKaUHX-Ek-h3j6GpmtrDeyTsOAgPbGp?usp=share_link)下載程式碼，這是一個幾本的矩陣計算，我們將利用 local memory 的特性改進之 。

## 程式碼簡介

這是矩陣計算的 kernel 部份，M、N 和 K 為矩陣 A、B 和 C 的大小，直接看程式碼可能難以理解，我們來看看示意圖

__kernel void myGEMM(const int M, const int N, const int K,
                     const __global float* A,
                     const __global float* B,
                     __global float* C) {
    
    // Thread identifiers
    const int globalRow = get_global_id(0); // Row ID of C (0..M)
    const int globalCol = get_global_id(1); // Col ID of C (0..N)

    // Compute a single element (loop over K)
    float acc = 0.0f;
    for (int k=0; k<K; k++) {
        acc += A[k*M + globalRow] * B[globalCol*K + k];
    }
     
    // Store the result
    C[globalCol*M + globalRow] = acc;
}

</br>

圖為 A 和 B 的矩陣相乘示意圖，M 和 K 為矩陣 A 的大小，K 和 N 為矩陣 B 的大小，M 和 N 為兩矩陣相乘後矩陣 C 的大小。


（圖片[來源](https://github.com/iVishalr/GEMM)）

</br>

編譯 kernel 時，可以通過 ```build()``` 加入特定 Macro 到 kernel 內，下列部份為加入 ``` -DTS=32``` 之 Macro 到 OpenCL 內，加入後 kernel 就可以直接使用 ```TS```，無須再次定義。此參數等等改進時會用到。

    void Initialize(std::string ker, std::string name) {
        program_ = cl::Program(context_, ker);
        std::ostringstream ss;
        ss << " -DTS=" << TS;

        auto args = ss.str();
        if (program_.build(args.c_str()) != CL_SUCCESS) {
            throw std::runtime_error("Fail to build the kernel.");
        }
        '''
    }               

</br>

此部份程式碼決定 global size 和 local size 的大小，每一個 work group 大小為 ```TS``` x ```TS```， 而每個 work group 會共用相同的 local memory，此特性等等改進時會用到。

    void Run(int M, int N, int K, float *A, float *B, float *C) {

        '''
        
        // Run the kernel.
        queue_.enqueueNDRangeKernel(
            kernel_, cl::NullRange,
            cl::NDRange(M, N),
            cl::NDRange(TS, TS)
        );
        
        '''
    }

</br>

## 改進 kernel 的缺陷

之前有提到過，如果直接多次讀取 global memory 是沒有效率的行為。這裡先回想一下矩陣相乘的方法，是否矩陣乘法會多次用到相同元素，很明顯，矩陣相乘會使矩陣內同一元素會被讀取多次，例如，矩陣 A 的第一列會被重複讀取 N 次 (矩陣 B 的寬度)，因此相乘前先將部份元素寫入 local memory，之後就多次讀取 local memory，如此一定程度上能減少記憶體延遲。但具體上要怎樣做呢？

大學或高中的時候，應該有學到一個矩陣計算的小技巧，就是將一個矩陣劃分成多個區塊，之後看成一個新的矩陣，而最後相乘結果等於此區塊相乘結果，如下圖，將矩陣 A 劃分成 a、b、 c 和 d 四個小矩陣，而相乘後的 C 等效於劃分矩陣的相乘。

（圖片[來源](https://hackmd.io/@siahuat0727/HJscCbWgV?type=view)）

</br>

我們可以利用這個特性達到加速效果，假設我要求 $C_{0,0}$ 的值，首先先將 a 和 e 矩陣讀進 local memory ，之後計算 ae 矩陣位於 (0, 0) 的值並紀錄下來，假設稱為 $C^{0}_{0,0}$，接下來將 b 和 g 矩陣讀進 local memory ，之後計算 bg 矩陣位於 (0, 0) 的值並紀錄下來，假設稱為 $C^{1}_{0,0}$，最後所求 $C_{0,0} = C^{0}_{0,0} + C^{1}_{0,0}$。如此可以減少讀取的延遲時間。

下為改進的 kernel ，其中 ```__local``` 代表宣告的為 local memory，```Asub``` 和 ```Bsub``` 為劃分的小矩陣，我們就是要將 global memory 的值複製到這裡，再進行讀取。而```barrier(CLK_LOCAL_MEM_FENCE)``` 代表同步 work group 內所有的 work items（或稱 threads）。

__kernel void myGEMM(const int M, const int N, const int K,
                     const __global float* A,
                     const __global float* B,
                     __global float* C) {
    
    // Thread identifiers
    const int row = get_local_id(0); // Local row ID (max: TS)
    const int col = get_local_id(1); // Local col ID (max: TS)
    const int globalRow = TS*get_group_id(0) + row; // Row ID of C (0..M)
    const int globalCol = TS*get_group_id(1) + col; // Col ID of C (0..N)
 
    // Local memory to fit a tile of TS*TS elements of A and B
    __local float Asub[TS][TS];
    __local float Bsub[TS][TS];
 
    // Initialise the accumulation register
    float acc = 0.0f;
    
    // Loop over all tiles
    const int numTiles = K/TS;
    for (int t=0; t<numTiles; t++) {
 
        // Load one tile of A and B into local memory
        const int tiledRow = TS*t + row;
        const int tiledCol = TS*t + col;
        Asub[col][row] = A[tiledCol*M + globalRow];
        Bsub[col][row] = B[globalCol*K + tiledRow];
 
        // Synchronise to make sure the tile is loaded
        barrier(CLK_LOCAL_MEM_FENCE);
 
        // Perform the computation for a single tile
        for (int k=0; k<TS; k++) {
            acc += Asub[k][row] * Bsub[col][k];
        }
 
        // Synchronise before loading the next tile
        barrier(CLK_LOCAL_MEM_FENCE);
    }
 
    // Store the final result in C
    C[globalCol*M + globalRow] = acc;
}

</br>

## 改進前後對比

在我的電腦上直接執行改進前後的程式碼，執行時間為循環 10000 的結果，我的 GPU 為 RTX 2070s，可以看到性能有一定的提昇

改進前 ```Sec: 8.58197```
改進後 ```Sec: 7.92214```

但由於此程式碼，包含複製資料到 GPU 的時間，此時間在我的電腦上大概為 5 秒，扣掉複製的時間後，可以看到如此簡單的改進就有不少性能提昇。

改進前 ```Sec: 3.58197```
改進後 ```Sec: 2.92214```

## 後記

此 kernel 來自於[這篇](https://cnugteren.github.io/tutorial/pages/page1.html)教學，想深入了解其它加技術的讀者可以閱讀這系列文章，此作者非常厲害，世界最快的 OpenCL Blas 庫是這位先生寫的。

## Linking

[上一篇](https://hackmd.io/@yrHb-fKBRoyrKDEKdPSDWg/BJCUHajju)
下一篇