The Power of Prefetching

這篇其實比"你也可以寫 SIMD 比網頁還快"早很多, 其中 I 的部分就是與這篇的  function 比較速度
====

先前在臉書上有分享過 “When Prefetching Works, When It Doesn’t, and Why”這篇 paper, 現在要用實例顯示 prefetching 帶來的效益, 今日使用的範例是 Matrix Transpose

一個簡單的 C 版本實作如下

void naive_transpose(int *src, int *dst, int w, int h){

    for(int x = 0; x < w; x++){

        for(int y = 0; y < h; y++){

             *(dst + x*h + y) = *(src + y*w + x);

        } 

    }

} 

接著請參考 http://www.randombit.net/bitbashing⋯⋯

SSE2 實作版本如下, 由於先進處理器架構有著 automatic/speculative prefetcher 所以這裡我們故意寫得比較 CPU-unfriendly 一些

void sse_transpose(int *src, int *dst, int w, int h){ 

    for(int x = 0; x < w; x+=4){

        for(int y = 0; y < h; y+=4){

            __m128i I0 = _mm_loadu_si128 ((__m128i*)(src+y*w+x)); 

            __m128i I1 = _mm_loadu_si128 ((__m128i*)(src+(y+1)*w+x));

            __m128i I2 = _mm_loadu_si128 ((__m128i*)(src+(y+2)*w+x)); 

            __m128i I3 = _mm_loadu_si128 ((__m128i*)(src+(y+3)*w+x));

            __m128i T0 = _mm_unpacklo_epi32(I0, I1); 

            __m128i T1 = _mm_unpacklo_epi32(I2, I3); 

            __m128i T2 = _mm_unpackhi_epi32(I0, I1);

            __m128i T3 = _mm_unpackhi_epi32(I2, I3); 

            I0 = _mm_unpacklo_epi64(T0, T1); 

            I1 = _mm_unpackhi_epi64(T0, T1); 

            I2 = _mm_unpacklo_epi64(T2, T3); 

            I3 = _mm_unpackhi_epi64(T2, T3); 

            _mm_storeu_si128((__m128i*)(dst+(x*h)+y), I0); 

            _mm_storeu_si128((__m128i*)(dst+((x+1)*h)+y), I1); 

            _mm_storeu_si128((__m128i*)(dst+((x+2)*h)+y), I2); 

            _mm_storeu_si128((__m128i*)(dst+((x+3)*h)+y), I3);

        } 

    }

}

接著我們來撰寫一個使用 SSE prefetch 指令的加速版本

void sse_prefetch_transpose(int *src, int *dst, int w, int h){ 

    for(int x = 0; x < w; x+=4){

        for(int y = 0; y < h; y+=4){

            #define PFDIST  8
            _mm_prefetch(src+(y+PFDIST)*w+x, _MM_HINT_T1); 

            _mm_prefetch(src+(y+PFDIST+1)*w+x, _MM_HINT_T1);

            _mm_prefetch(src+(y+PFDIST+2)*w+x, _MM_HINT_T1);

            _mm_prefetch(src+(y+PFDIST+3)*w+x, _MM_HINT_T1);

 
            __m128i I0 = _mm_loadu_si128 ((__m128i*)(src+y*w+x)); 

            __m128i I1 = _mm_loadu_si128 ((__m128i*)(src+(y+1)*w+x)); 

            __m128i I2 = _mm_loadu_si128 ((__m128i*)(src+(y+2)*w+x)); 

            __m128i I3 = _mm_loadu_si128 ((__m128i*)(src+(y+3)*w+x)); 

            __m128i T0 = _mm_unpacklo_epi32(I0, I1);

            __m128i T1 = _mm_unpacklo_epi32(I2, I3);

            __m128i T2 = _mm_unpackhi_epi32(I0, I1);

            __m128i T3 = _mm_unpackhi_epi32(I2, I3); 

            I0 = _mm_unpacklo_epi64(T0, T1); 

            I1 = _mm_unpackhi_epi64(T0, T1); 

            I2 = _mm_unpacklo_epi64(T2, T3); 

            I3 = _mm_unpackhi_epi64(T2, T3); 

            _mm_storeu_si128((__m128i*)(dst+(x*h)+y), I0); 

            _mm_storeu_si128((__m128i*)(dst+((x+1)*h)+y), I1); 

            _mm_storeu_si128((__m128i*)(dst+((x+2)*h)+y), I2); 

            _mm_storeu_si128((__m128i*)(dst+((x+3)*h)+y), I3); 

        }

    }

}

接著是測試程式(請將上述程式碼置於指定位置)

#include  

#include 

#include  

#include 
//for using x86-SSE2 intrinsics you have to include this

#include 
 

// PUT naive_transpose, sse_transpose, sse_prefetch_transpose HERE 

// ...
 

#define TEST_W 4096

#define TEST_H 4096
int main(void){

    //verify the result of 4x4 matrix 

    { 

        int testin[16] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }; 

        int testout[16];
        for(int y = 0; y < 4; y++){

            for(int x = 0; x < 4; x++)

                printf(" %2d", testin[y*4+x]);

            printf("\n");

        } 

        printf("\n");

        sse_transpose(testin, testout, 4, 4); 

        for(int y = 0; y < 4; y++){ 

            for(int x = 0; x < 4; x++) 

                printf(" %2d", testout[y*4+x]); 

            printf("\n");

        }

    }
 

    {

        struct timeval stime, etime; 

        int *src = (int*)malloc(sizeof(int)*TEST_W*TEST_H);

        int *out0 = (int*)malloc(sizeof(int)*TEST_W*TEST_H); 

        int *out1 = (int*)malloc(sizeof(int)*TEST_W*TEST_H); 

        int *out2 = (int*)malloc(sizeof(int)*TEST_W*TEST_H);

        srand(time(NULL));

        for(int y = 0; y < TEST_H; y++){ 

            for(int x = 0; x < TEST_W; x++)

                *(src + y*TEST_W + x) = rand();

        }
        gettimeofday(&stime, NULL); 

        sse_prefetch_transpose(src, out0, TEST_W, TEST_H); 

        gettimeofday(&etime, NULL);

        printf("sse prefetch: %ld us\n", (etime.tv_sec - stime.tv_sec)*1000000 + (etime.tv_usec - stime.tv_usec));
 

        gettimeofday(&stime, NULL); 

        sse_transpose(src, out1, TEST_W, TEST_H);

        gettimeofday(&etime, NULL); 

        printf("sse: %ld us\n", (etime.tv_sec - stime.tv_sec)*1000000 + (etime.tv_usec - stime.tv_usec));
 

        gettimeofday(&stime, NULL); 

        naive_transpose(src, out2, TEST_W, TEST_H);

        gettimeofday(&etime, NULL); 

        printf("naive: %ld us\n", (etime.tv_sec - stime.tv_sec)*1000000 + (etime.tv_usec - stime.tv_usec));

        free(src);

        free(out0); 

        free(out1); 

        free(out2); 

     }
     return 0; 

}

存成 test.c, 並且按照下列方式編譯

gcc -msse2 -o test test.c

於個人的 Core i7 3612QM 執行得到下列結果

  0  1  2  3 

  4  5  6  7

  8  9 10 11 

 12 13 14 15
 

  0  4  8 12 

  1  5  9 13 

  2  6 10 14

  3  7 11 15
 

sse prefetch:  57782 us

sse: 117184 us 

naive: 228870 us

儘管現今 CPU 有著強大的 auto-prefetcher, 但是並不是所有的演算都有著規律或是線型的記憶體存取模式, 如此透過優化 prefetching 依然是能夠有相當的效能增進.

你也可以寫 SIMD 比寫網頁還快 - II

第一篇並沒有介紹 clang 中的 OpenCL kernel 一些基礎 型別宣告部分, 已經在上一篇出現過了, 所以這裡就從使用開始. 首先是 vector 的初始化

typedef int int4 __attribute__((ext_vector_type(4)));

typedef int int8 __attribute__((ext_vector_type(8)));

int a = 3;

int4 va = {1, 2, 3, 4}; 

int4 vb = 3;

接著是個別的 access, 下列的方式很適合做資料格式上的操作 在等號的兩邊可以用不同的表示法, 沒有需要一致:

va.x = a; 

// the below two lines are equivalent. 

vb.xyw = va.xzw; 

vb.s013 = va.s023; 

// repeat is ok 

va = vb.xyyz; 

// high and low part 

va.hi = va.lo;

va.lo = vb.hi;

// even and odd 

va.odd = vb.even;

Vector 簡單的計算如下, 可以參照 clang - vectors and extended vectors 的列表:

// scalar multiply 

va = a * vb;

// element-wise multiply

va = va * vb;

va = va | vb;

va = -va;

不同型別的 vector 轉換:

typedef float float4 __attribute__((ext_vector_type(4))); 

// format convertion 

float4 vfa = __builtin_convertvector(va, float4); 

// data reinterpret

float4 vfa = (float4)va;

接著來實戰一番吧, 下列為以 bayer format 平均 GRBG 4個像素, 轉為1/4大小灰階圖的範例

void bayer_convert(char *src, char *gray, int w, int h) 

{

    for(int y = 0; y < h; y+=2){ 

        for(int x = 0; x < w; x+=16){

            //bayer raw data fetch

            ushort8 pix00 = __builtin_convertvector(*((uchar8*)(src+x)), ushort8);

            ushort8 pix01 = __builtin_convertvector(*((uchar8*)(src+x+8)), ushort8); 

            ushort8 pix10 = __builtin_convertvector(*((uchar8*)(src+x+w)), ushort8);

            ushort8 pix11 = __builtin_convertvector(*((uchar8*)(src+x+w+8)), ushort8);

            ushort8 pix_g0;

            pix_g0.lo = pix00.even;

            pix_g0.hi = pix01.even; 

            ushort8 pix_r;

            pix_r.lo = pix00.odd;

            pix_r.hi = pix01.odd;

            ushort8 pix_b;

            pix_b.lo = pix10.even;

            pix_b.hi = pix11.even;

            ushort8 pix_g1;

            pix_g1.lo = pix10.odd;

            pix_g1.hi = pix11.odd;

            // average! so simple

            ushort8 pix_out = (pix_g0 + pix_r + pix_b + pix_g1) ;

            // write out 

            *((uchar8*)(gray+(x/2))) = __builtin_convertvector(pix_out, uchar8); 

        }

        src += w*2; 

        gray += w;

    }

}

如何? 比起看 intrinsics function 來得清楚與簡單吧

如果還想了解 OpenCL vector 的操作可以參考這篇文章

有興趣記得編譯時加入 --save-temps 來看使用的指令!!! 下一篇會介紹進階的使用方式

你也可以寫 SIMD 比寫網頁還快 - I

看到 Jim Huang 要開”從無到有自幹 boot loader 比寫網頁還快” 除了幫他打廣告外, 也趕快搭這個熱潮, 證明人人可以寫 SIMD 比網頁還快 去年寫了系列 MMX/SSE 去處理 Matrix Transpose 但是看到 SSE/AVX intrinsics 很多人大概就眼睛花了 若加上實際應用還要從 Intel 網頁查找適當的 intrinsics 會覺得難以下手 搭配適當的工具, 你也可以是 SIMD 高手! 這裡最重要的工具是 - clang clang 如同 gcc 一般提供了許多 extension 其中很重要的一點是 Vectors and Extended Vectors 這個 extension 請特別注意 link 中的 OpenCL vector, 底下的表格表示此方式提供相當全面的 operator 而 clang 中使用 OpenCL vector 這即是這裡所推荐的方式 (還可以學會 OpenCL kernel 的部份撰寫!) 首先來以 clang 方式定義 OpenCL vector 的型別

typedef int int8 __attribute__((ext_vector_type(8)));

UPDATED: 對於 GCC 支援相近語法, 其宣告為 (其中 vector_size 以 byte 計)

typedef int int8 __attribute__ ((vector_size (32)));

這裡定義了長度為 8 的 int vector type 接著我們來實作先前時做過的 matrix transpose 功能

void transpose8x8(int *src, int *dst, int w, int h){
    for(int x = 0; x < w; x+=8){
        for(int y = 0; y < h; y+=8){
            int8 row0 = *((int8*)(src+y*w+x));
            int8 row1 = *((int8*)(src+(y+1)*w+x));
            int8 row2 = *((int8*)(src+(y+2)*w+x));
            int8 row3 = *((int8*)(src+(y+3)*w+x));
            int8 row4 = *((int8*)(src+(y+4)*w+x));
            int8 row5 = *((int8*)(src+(y+5)*w+x));
            int8 row6 = *((int8*)(src+(y+6)*w+x));
            int8 row7 = *((int8*)(src+(y+7)*w+x));
            int8 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
            tmp0 = (int8){row0.s0, row1.s0, row2.s0, row3.s0, row4.s0, row5.s0, row6.s0, row7.s0};
            tmp1 = (int8){row0.s1, row1.s1, row2.s1, row3.s1, row4.s1, row5.s1, row6.s1, row7.s1};
            tmp2 = (int8){row0.s2, row1.s2, row2.s2, row3.s2, row4.s2, row5.s2, row6.s2, row7.s2};
            tmp3 = (int8){row0.s3, row1.s3, row2.s3, row3.s3, row4.s3, row5.s3, row6.s3, row7.s3};
            tmp4 = (int8){row0.s4, row1.s4, row2.s4, row3.s4, row4.s4, row5.s4, row6.s4, row7.s4};
            tmp5 = (int8){row0.s5, row1.s5, row2.s5, row3.s5, row4.s5, row5.s5, row6.s5, row7.s5};
            tmp6 = (int8){row0.s6, row1.s6, row2.s6, row3.s6, row4.s6, row5.s6, row6.s6, row7.s6};
            tmp7 = (int8){row0.s7, row1.s7, row2.s7, row3.s7, row4.s7, row5.s7, row6.s7, row7.s7};
            *((int8*)(dst+(x*h)+y)) = tmp0;
            *((int8*)(dst+((x+1)*h)+y)) = tmp1;
            *((int8*)(dst+((x+2)*h)+y)) = tmp2;
            *((int8*)(dst+((x+3)*h)+y)) = tmp3;
            *((int8*)(dst+((x+4)*h)+y)) = tmp4;
            *((int8*)(dst+((x+5)*h)+y)) = tmp5;
            *((int8*)(dst+((x+6)*h)+y)) = tmp6;
            *((int8*)(dst+((x+7)*h)+y)) = tmp7;
        }
    }
}

將上列內容與先前的範例結合, 最後就是編譯了

clang -O3 -mavx -mfma -o test test.c

 以下為在 Intel(R) Core i7-3612QM 的執行結果

clang vector: 35759 us
sse: 41246 us
naive: 146310 us

這樣是不是很快呢? 甚至比之前直接用 SSE intrinsics 寫的還快!!!(灑花) 更棒的是這樣的 SIMD 撰寫, 透過 clang 是可以在各種硬體平台的 可以自己試試看用 NDK 然後跑在自己的手機看看喔!!

The era of OpenCL 2.0 is comming

Although OpenCL 2.0 spec was finalized & released in Dec 2013. But it takes almost one year for vendors to implement 2.0 runtime. In the recent months, Intel and AMD have released their SDK (Intel SDK, AMD SDK) for OpenCL 2.0.

AMD provides a series of articles about the new features of OpenCL 2.0 on the AMD Developer Blog:
OpenCL 2.0 - Shared Virtual Memory
OpenCL 2.0 - Pipes
OpenCL 2.0 - Device Enqueue and Workgroup Built-in Functions
OpenCL 2.0 - Generic Address Space and Program Scope Variables
OpenCL 2.0 - Image Enhancement

Besides these, AMD also provides a source code package of OpenCL 2.0 Samples. It is worth to take a look. You can get the key points of OpenCL 2.0 from the topics of sample package and articles. IMHO, SVM, Pipes and Device Enqueue are the most important three features among the new features of OpenCL 2.0.

==
儘管 OpenCL 2.0 規格在2013年十二月就發佈了. 但是各家廠商花了近一年實作 OpenCL 2.0 runtime. 最近幾個月Intel 與 AMD 相繼釋出各自的 SDK (Intel SDK, AMD SDK) for OpenCL 2.0.

AMD 在其 Developer Blog 還發佈了一系列關於 OpenCL 2.0 新特性的文章:
OpenCL 2.0 - Shared Virtual Memory
OpenCL 2.0 - Pipes
OpenCL 2.0 - Device Enqueue and Workgroup Built-in Functions
OpenCL 2.0 - Generic Address Space and Program Scope Variables
OpenCL 2.0 - Image Enhancement

此外, AMD 也提供了 OpenCL 2.0 Samples 的原始碼, 值得一看. OpenCL 2.0 的關鍵點能透過這些原始碼與文章的標題觀察得到. 其中 SVM, Pipe 與 Device Enqueue 是個人認為 OpenCL 2.0 新特性之中最重要的3個.

Gaming 或許是 Android 5.0 Lollipop 其中一項重點

Android 5.0 在近日正式釋出了
而在 10/17 也已釋出了 SDK for Android 5.0
對於 Android 5.0 網路上相關介紹文件不少
由於先前並沒有花心思看 Android L Preview 的介紹
所以這幾日花了點時間消化相關新的系統特性

這裡想指出 Google 對於 Android 平台在 Gaming 的野心
首先是 Android 5.0 除了 OpenGL ES 3.1 的支援外提出的 Android Extension Pack(AEP)
The extension pack supports:
* Guaranteed fragment shader support for shader storage buffers, images, and atomics (Fragment shader support is optional in OpenGL ES 3.1.)
* Tessellation and geometry shaders * ASTC (LDR) texture compression format
* Per-sample interpolation and shading
* Different blend modes for each color attachment in a frame buffer

其中前兩點是最重要的,
請搜尋 Tessellation 稍微了解其目的
另外必須知道 Tessellation 是 OpenGL 4.0/DirectX 11 後的特性
Tessellation & Geometry Shader 是近年 GPU 提升模型動態細緻化的重大演進
對於 AEP 詳細的規範可以參考 Khronos 的 spec

而另外一個點在於此次 Google 裝置除了 Nexus 6/9 外
還有個不太受到注目的 Nexus Player,
多數的 Nexus Player 分析著重於 Android TV 與具有語音控制/搜尋的遙控器上, 這點無可厚非
但 Ｎexus Player 有許多有趣的特性(像是它是第一個官方 x86 Nexus 裝置)
其一便是 Nexus Player 提供了選購的無線藍牙 GamePad
目前已知 Nexus Player 中的 Intel Atom 所使用的 GPU 是 PowerVR 6000 系列
這系列的 GPU 俱備了支援 AEP 的規格
此兩點顯示 Google 有意透過 Nexus Player 做為 Android Gaming Console 的示範平台的野心
儘管 Android Extension Pack 對於 Android 並非必要項目,
但在智慧裝置利潤逐漸下滑的今日, 平台差異化顯得重要
AEP這樣明確規範的規格, 是很容易成為各家IC生產商軍備競賽的比較條件之一...
況且, 若無市場戰略目標, 維持目前 Android 圖形需求, 跟隨 OpenGL ES 標準也已足夠 (若要與 iOS Metal 競爭那又是另外一個面向的事情了)
並無制定這套延伸 OpenGL ES 擴充標準的必要

OpenCL Kernel Optimization Tips 簡報上線

這段日子工作與 OpenCL 相關
因此碰觸到了一些 Mobile 與 Desktop 的 OpenCL Runtime
面對 OpenCL 程式對不同的硬體平台做優化的過程
有許多的層面必須要考量, 因此將這些想法集結為這份 slides
日後希望有機會能update加入 example 與 case study.

OpenCL Kernel Optimization Tips from Champ Yen

OpenCL Utility/Framework 專案 - CLScript

開了自己的第一個 github project
https://github.com/champyen/clscript

這一兩年由於 GPGPU 應用的興起,
心中一直有個關於 video codec 的想法想要付諸實現
近日開始著手動工, 儘管有些地方的設計還沒頭緒
但是基本的問題還是需要解決
CLScript 這個專案因此這些基本問題而產生, 目前百廢待興, 文件與 code 還需要作進一步改進, 但是初版已可以玩一玩
OpenCL 其中一個讓 programmer 感到麻煩的地方在於
OpenCL Host Code 的撰寫很冗長,
且對於一個 function call 需要相當的前置作業
撰寫 kernel code 對應的 Host Program 是無趣的一件事
CLScript 這是個 utility/framework 層面的專案
其目的就在於簡化 OpenCL 的使用流程, 解決上述問題
讓開發者能更專注於 OpenCL 中核心的 kernel code 的開發
C/C++ framework 使用上, 建構的概念是將 .cl file 視為類似於 C library
在 test/test.cpp 中可以看到參考用法
using namespace clscript;
...
CSRuntime runtime(CS_DEV_GPU);
CSLib lib(&runtime, clFileName);
CSBuffer buf(&runtime, 4096*sizeof(cl_int));
CSWorkSize gWS(4096);
lib.exec("test", gWS, NullWorkSize, 0, &buf);

CSLib 的 exec 界面會將傳進的不定數目參數作設定與傳遞
對 C/C++ 使用者而言, 就像是呼叫一個特定的 function
另外一個 exec 界面是為了 csutil CLI tool 而設的
csutil 為直接測試 CL kernel code 的工具, 可以不用寫 Host Program 而直接測試kernel code 程式
csutil 有四大參數, 目前預設使用 GPU device
1. CL 檔案
2. 在 1. 檔案中, 所要使用的 function 名稱
3. WorkSize
     各 dimension 以 ',' 區隔, global 與 local 以 ':' 區隔, local worksize部份並非必要
4. function 所使用的參數
     * float - f:1.0
     * double - d:1.0
     * int - i:1
     * buffer - b:sz=4096,if=/dev/zero,of=out.bin
    這是最麻煩的, 設定參數又有三個sz, if 與 of(皆非必要), sz 為 buffer 大小, 未設置預設為 global worksizes 相乘, if 為用來初始化 buffer 的資料, 會讀入與 buffer size 相同大小的內容, 若無必要初始化無須指定, of為 buffer 結果的輸出, 若為 temp buffer 無須指定

以 source 中附上的 test.cl 就可以這樣使用
./csutl/csutil ../test.cl test 4096 i:2 b:sz=16384,if=/dev/zero,of=out.bin