Put It Altogether! : Matrix Multiplication with libdispatch + Clang SIMD optimization

先前介紹過 Clang OpenCL vector 與 libdispatch
這裡以 Matrix Multiplication 為例將這兩者結合, 來顯示這樣作法的強大

相關程式碼可以透過 github 取得

git clone https://github.com/champyen/gcd_vector.git 

首先在 matmul8x8.h 中我們定義了需要 Matrix Multiplication 8x8 tiling 的介面

void matmul8x8(float *m1, int p1, float *m2, int p2, float *mo); 

其中 m1, m2 為需要相乘的兩個 8x8 矩陣, 而 mo 為 output
而 p1 為 m1 的 pitch size, p2 為 m2 的 pitch size

因此在 test.c 中, 我們透過這介面實作了Matrix Multiplication 的主 flow

    for(int y = 0; y < H1; y+=8){
        for(int x = 0; x < W2; x+=8){
            matmul8x8(m1+y*W1, W1, m2+x, W2, mo+y*W2+x);
        }
    }

接著實作了三個版本, 分別是 naive, tile 與 vec
naive.c: 這個版本單純是直接計算出 mo 於該位置的 output

void matmul8x8(float *m1, int p1, float *m2, int p2, float *mo)
{
    for(int y = 0; y < 8; y++){
        for(int x = 0; x < 8; x++){
            float val = 0.0f;
            for(int idx = 0; idx < p1; idx++){
                val += (*(m1+y*p1+idx)) * (*(m2+idx*p2+x));
            }
            *(mo+y*p2+x) = val;
        }
    }
}

tile.c: 這個版本中, 以 8x8 矩陣相乘為單位的方式累加出 mo 的結果

void matmul8x8(float *m1, int p1, float *m2, int p2, float *mo)
{
    for(int y = 0; y < 8; y++){
        for(int x = 0; x < 8; x++){
            *(mo+y*p2+x) = 0.0f;
        }
    }

    for(int idx = 0; idx < p1; idx+= 8){
        for(int y = 0; y < 8; y++){
            for(int x = 0; x < 8; x++){
                for(int idx2 = 0; idx2 < 8; idx2++){
                    *(mo+y*p2+x) += (*(m1+y*p1+(idx+idx2))) * (*(m2+(idx+idx2)*p2+x));
                }
            }
        }
    }
}

 vec.c: 這是個 tile.c 再加上使用 Clang OpenCL vector

typedef float float8 __attribute__((ext_vector_type(8)));

#define PF_DIST 8
void matmul8x8(float *m1, int p1, float *m2, int p2, float *mo)
{
    float8 vrow[8];
    for(int y = 0; y < 8; y++){
        vrow[y] = (float8)0.0f;
    }

    for(int idx = 0; idx < p1; idx+= 8){
        float8 vm2[8];
        int ofs = idx*p2;
        int ofs2 = PF_DIST*p2;
        for(int y = 0; y < 8; y++, ofs+=p2){
            vm2[y] = *((float8*)(m2+ofs));
            __builtin_prefetch(m2+ofs+ofs2);
        }
        ofs = 0;
        for(int y = 0; y < 8; y++, ofs+=p1){
            float8 vm1 = *((float8*)(m1+ofs+idx));
            __builtin_prefetch(m1+ofs+idx+PF_DIST);
            for(int idx2 = 0; idx2 < 8; idx2++)
                vrow[y] += vm1[idx2]*vm2[idx2];
        }
    }
    int ofs = 0;
    for(int y = 0; y < 8; y++, ofs+=p2){
        *((float8*)(mo+ofs)) = vrow[y];
    }
}

接著我們導入使用 libdispatch 來做多核運算加速
於 test_gcd 中我們改寫了矩陣相乘的 flow

    for(int y = 0; y < H1; y+=8){
        for(int x = 0; x < W2; x+=8){
            dispatch_group_async(
                group,
                queue,
               ^{
                    matmul8x8(m1+y*W1, W1, m2+x, W2, mo+y*W2+x);
                }
            );
        }
    }
    dispatch_group_wait(group, DISPATCH_TIME_FOREVER);

最後是實際的測試, 這裡測試的是兩個 512x512 矩陣相乘
由於 libdispatch 與 Clang OpenCL vector 的可移植性, 所以我們能夠在多樣的平台上測試
於 Raspberry Pi 3 上, 透過 SIMD + GCD 甚至可以得到超過百倍的加速而架構間的數據倍率差異, 也很值得去思考結果為何如此

Raspberry Pi 3: 114X

• test 18700284 us
• test_tile 2160980 us
• test_vec 358326 us
• test_gcd 6778646 us
• test_gcd_tile 1085776 us
• test_gcd_vec 164121 us

A8-5545m: 15X

• test 446682 us
• test_tile 2613739 us
• test_vec 48856 us
• test_gcd 239297 us
• test_gcd_tile 2638397 us
• test_gcd_vec 29957 us

i7-3612QM: 23.4X

• test 275879 us
• test_tile 224157 us
• test_vec 31457 us
• test_gcd 112845 us
• test_gcd_tile 92521 us
• test_gcd_vec 11808 us

20170116 UPDATED:
後續加入了 GCC 的支援, 相關數字已更新於 github 頁面上