在 ARM 平台上使用 Function Multi-Versioning (FMV) - 以使用 Android NDK 為例

Function Multi-Versioning (FMV)

過往的 CPU 發展歷程中, x86 平台由於因應各種應用需求的提出, 而陸陸續續加入了不同的指令集, 此外也可能因為針對市場做等級區隔, 支援的數量與種類也不等. 在 Linux 平台上這些 CPU 資訊可以透過 /proc/cpuinfo 取得. 而對於 CPU 不同的指令集支援去實作軟體最佳化, 是相當繁複的工作. 而這些在過去分享過的 Function Multi-Versioning (FMV) 的出現, 而大幅地簡化.

隨著時間過去, 這十年來發展迅速的 ARM 平台, 在 ARMv8 之後不斷作小幅度的翻新, 也逐漸產生類似 x86 平台上的問題. 以 A-Profile 處理架構而言從 ARMv8.0-A 一路發展到 ARMv9.4-A. 當中有著許多新硬體特性與新增指令, 有些是 core feature 而有些則是 optional, 這些也逐漸造成了 ARM CPU 平台功能性上大小不等的碎片化. 在 ARM 官網有提供了一份 "Feature names in A-profile architecture" 列表. 而這整理僅止於 ARMv9 之前. 此外更麻煩的地方是, 各大 compiler 在先前一直沒有對應支援 FMV 的功能, 如此對於通用的應用程式, 想要能夠有效使用 ARM 這些功能特性同時顧及不支援的平台, 著實是一件不簡單的事.

而事情在 2023 年中開始有了轉機, 在當時所釋出的 Clang / LLVM 16.0 開始在 ARM 平台上支援  FMV 功能. ARM 官方甚至在其系列文的 "What's new in LLVM 16" 中有特別舉例說明 FMV 的使用. 一旦能夠使用 FMV, 這對於不確定運作的目標平台的應用程式開發很有幫助. 但在當時個人第一時間找不到具有支援的 llvm toolchain, 而之後的這段時間因工作繁忙一直沒有時間再繼續做嘗試, 也就沒辦法撰寫心得. 近日回想起心中所惦記的這麼一件事, 於是特別拿了最新的 Android NDK 26d (附的是 clang 17.0.2) 來測試.

Case Study - FP16 GEMM w/ Vector Extension, Gather Auto-Vectorization

首先請安裝 Android NDK 26d 到 /opt 目錄內, 並將 /opt/android-ndk-r26d/toolchains/llvm/prebuilt/linux-x86_64/bin/ 加入 PATH 環境變數中.

由於個人在過去所分享的教材投影片中即談論過 FMV, 本文中就不再多做詳細介紹, 有需要可以參考過去的投影片, 在這裡就透過兩個實際的例子來嘗試 ARM FMV 功能, 第一個 Case 是個人常拿來作為 Compiler Vector Extension 的 Laboratory 範例的 Matrix Multiplication, 這個範例可以顯示出使用 Compiler Vector Extension 搭配 FMV 可以共伴產生廣泛支援又高效的程式碼:

typedef __fp16 fp16_16 __attribute__((ext_vector_type(16)));
#ifdef __HAVE_FUNCTION_MULTI_VERSIONING
__attribute__((target_clones("sve2", "simd")))
#endif
void gemm_vec(
    // buffer, stride
    __bf16 *a, int sa,
    __bf16 *b, int sb,
    __bf16 *c, int sc
){
    fp16_16 vb[16];
    for(int y = 0; y < 16; y++){
        vb[y] = *((fp16_16*)(b + sb*y));
    }
    for(int y = 0; y < 16; y++){
        fp16_16 vc = *((fp16_16*)(c + sc*y));
        fp16_16 va = *((fp16_16*)(a + sa*y));
        for(int x = 0; x < 16; x++)
            vc += va[x]*vb[x];
        *((fp16_16*)(c + sc*y)) = vc;
    }
}

這個範例中, 為了能產生較大的 code generation 的差異, 個人選擇使用 float16 的資料型別. 主要是因為 ARMv8 並非預設支援 fp16, 而是到了 ARMv8.2 以 optional extension 來提供 (而預設支援 SVE2 的 ARMv9 是預設支援的). 將上述這段程式令存為 arm_gemm.c 之後搭配下列指令編譯:

$ aarch64-linux-android34-clang -O3 --save-temps -c arm_gemm.c

如此當中可以觀察產生的 arm_gemm.s 檔案, 當中有 gemm_vec._Msve2 與 gemm_vec._Msimd 兩段程式, 以及用來做選擇的 gemm_vec.resolver, 可以觀察 sve 版本直接使用支援 fp16 的 fmul, 而 AdvSIMD 版本當中不斷使用 fcvt 指令來處理浮點數轉換, 行數落差也相當的大.

接著的第二個範例中是透過 auto-vectorization 來實際展示 ARM 官網上的 "Auto-vectorization examples" 內關於 Scatter and Gather 的部份, 而為了能確實透過 clang/llvm 的組合產生使用 gather 指令, 該段程式經過稍加修改後程式碼如下:

#include <stdint.h>
#ifdef __HAVE_FUNCTION_MULTI_VERSIONING
__attribute__((target_clones("sve2", "simd")))
#endif
float gather_reduce(float *restrict data, int *restrict indices, long c)
{
    float r = 0;
    #pragma clang loop vectorize(enable)
    for (long i = 0; i < c; i++) {
        r += data[indices[i]];
    }
    return r;
}

將上述這段程式令存為 arm_gather.c 之後搭配下列指令編譯:

$ aarch64-linux-android34-clang -O3 --save-temps -c arm_gather.c

同樣地, 透過觀察編譯過程中產生的 arm_gather.s 檔案, 當中有 gather_reduce._Msve2 與 gather_reduce._Msimd 的兩段程式碼, 從中可以看到 SVE 的版本的 gather_reduce function 確實有使用 L1DW 這個 SVE 指令集所專屬的 gather 指令.

ARM Feature List

最後談談 ARM 與 x86 平台在 FMV 在實務使用中不同的地方. 在 ARM 平台的 FMV target 僅止為 feature, 而 x86 FMV 的 target 則能夠是特定的 CPU arch, 特定的 CPU core 或是指定的 ISA. 那麼具體來說在 ARM 上 FMV 有哪些 feature 可以選擇呢? 可以參考ARM 這份 ACLE 文件, 當中有提供對應的 feature 列表與在指定 FMV target 的名稱, 另外如以往的 FMV target 的使用, 這些 feature 是可以透過 + 來合併選擇使用的, 像是 "sve2+memtag" 來合併使用 SVE2 與 Memory Tagging, 同樣地一如以往 "default" 表示的是基本的指令集, 產生完全不會使用到 feature 的 code.

Chisel 學習筆記 - Scala 與 Chisel 基礎語法

標題為筆記, 但這篇比較屬於心得

延續上一篇的環境建立, 這次計劃藉由 Jserv 最新的課程安排來學習 Chisel, 當然個人目標是能夠按照 Jserv 的課程規劃在期限之內完成 Lab 3, 由於個人並非 digital designer (現在這年紀也算老貓學新把戲...), 因此對於 chisel 的一切也是從基本開始, 目前計劃於一週內按照 Chisel Bootcamp 的章節完成 Chap 0 ~ 3 的所有內容. (Chap 4 部份在後續會依照進階需求選擇回來看)

週末花了點時間把 Chap 0 ~ 2 的說明研讀完, code block 的操作也跑過一次流程, 另外是在 try and error 下完成當中的所有 exercise & test. 個人覺得 Chisel Bootcamp 多少是面向俱備 Verilog / HDL 有些許基礎的人. 儘管無基礎並不會造成學習的困難, 但多少難以體會 Chisel 相較之下的優點與強項.

基本上Chap 0 ~ 2 的重點其實並不多

• Chisel 的特性
• Scala 的基本語法
• 建構於 Scala 之上的 Chisel DSL
  
• 電路的基本: module, combinational, control & sequential logic
• 除錯與測試方式

由於並非第一次學習建構在另一語言的 Domain Specific Language. 過往學習 Halide 的過程就是這樣的一個經驗. 然而與先前 Halide 不一樣的地方是 Halide 所依附的 C++ 個人有所基礎, 但是對於 Scala / Chisel 這樣的組合卻是完全從頭學起. 而這大概有幾個面向是需要反覆思考與刻意練習的.

第一個主要的點在於 Value 與 Typing System, 類似於 C++ 與 Halide 有著各自的 value & typing system, C++ 主要用在撰寫與控成 Halide code 的流程, Halide 體系著重在產生 target code 上. Scala 與 Chisel 也有各自的 value & typing system. Scala 也是用在撰寫 Chisel 主體與流程上, Chisel 是用以定義 hardware 使用到的 module 與產生 instance. Chap 1 ~ 2.1 強調分清楚 Scala 與 Chisel 的 value & typing system 確實是很重要的一件事.

第二個是 Chisel 的實作撰寫必需了解的 3 個階段 - circuit generation, circuit simulation 與 testing. (Halide 也有著類似的階段, 甚至都使用 Generator 這樣的名稱) Scala 與 Chisel 在使用上容易搞混的是 generation, Scala 是讓 code generation 能透過參數化來產生對應不同的 circuit (按照教學說法, 這是 Chisel 強大的地方). 因此這些都是在 compile time 決定的, 而 Chisel 的撰寫則是決定了 circuit generation 的結果. 對於一些 if / else 的思考必須區隔希望影響的是 compile time 的行為抑或者是 runtime 時的所需. (習題上來說, 觀察前後輸入的 data type 就可以快速反應)

儘管這些並不是什麼艱深難懂的心得, 然而困難的點是在熟悉與快速應用 Scala / Chisel. 重要的還是需要刻意練習來讓思維方式有著正確的對應. 若儘只是研讀一番並不足以深刻體會, 或是習得實務上的撰寫能力. Chisel Bootcamp 的作業與測試非常值得一一花時間下去練習. 藉由這些能夠一次次找到認知錯誤的地方或是盲點. 個人認為能嘗試完成 2.2 Arbiter, 2.3 Finite State Machine 與 2.5 FIR Filter Generator 這三個 exercise 應符合了基本要求. ( 2.2 中提供的 Cheatsheet 第一頁對於閱讀內容與作業相當有幫助)

OpenCL 1.2 之後(個人觀點中)重要的 feature 與 extension

最近因工作緣故, 而看到了使用 subgroup 特性的實作, 在第一時間個人難以解讀程式碼的目的與原意,  後來看了原始實作與擴充功能的說明文件才得知 subgroup 的功能特性為何. 主要因為個人對於 OpenCL 2.x 印象僅停留在 Shared Virtual Memory. 然而事實上 OpenCL 2.0 到 OpenCL 3.0 這段期間依然有數個實用的功能特性與擴充 對於 OpenCL 3.0 而言, 這些 2.0 以後加入的新 features 與 extensions 都屬於 extension. 因此可以到 The OpenCL™ Extension Specification 頁面一一查看. 以下列出個人覺得實用的 extensions:

Creating a 2D Image From A Buffer

許多時候人們很希望能透過 vload / vstore 來存取 buffer, 又希望能夠使用 GPU 的 texturing unit. 在 OpenCL 1.x 的時代, 通常是透過 ION buffer 來重複 mapping 到 Buffer 與 Image, 並且傳入對應的 kernel 中使用. 基本上這需要 GPU vendor 確認 cache coherency 的問題, 甚至提供對應的 driver 支援. 在 OpenCL 2.0 中將類似功能成為 core feature, 也就是 OpenCL 2.0 中是必要的存在. 透過 cl_khr_image2d_from_buffer 這個 extension, 可以先配置 Buffer 後再以此 Buffer 來建立 Image. 一開始提到相同目的.

Subgroups

儘管 OpenCL 中提供了 Global 與 Local work size 來對整體工作做切割, 而 Local work size 另一個重要的意義是對應的 workgroup size 是能作同步的最小單位. 在 OpenCL 2.1 中新增了這個 core feature. 在 Workgroup 中能夠再切分一個 1D 的 subgroup, 除了 subgroup 層及的 barrier 做同步外, 也提供了 3 個種類的 kernel function: reduce_op, scan_exclusive_op 與 scan_inclusive_op. 這些 operation 對於 reduction 與 integral 這兩類原本 OpenCL 不適合應付的計算提供了有效的介面, 讓 device 能夠使用硬體支援來處理這類的計算.

Mipmaps

在 OpenGL 中支援以硬體產生 Mipmaps 來做 texturing. 可以避免 aliasing. 這裡應該是同等的結構. 但除了做原本 antialiasing 用途外, 不知是否有機會利用這個 pyramid 結構來做 CV 應用.

Integer Dot Product

OpenCL 中原本只支援 float / double 型別的內積 (dot product), 然而現今 NN 應用中 8-bit integer dot product 的需求相當大. 這個 extension 有效的提供了對於 8-bit dot product 硬體加速的介面. 理論上透過 multiplication 搭配 subgroup 就可以得到相同結果, 但若支援此 extension 且俱備硬體加速, 可以得到更高的 computation throughput.

當然 OpenCL 2.x 系列還是有其他像是 SVM, device-side enqueue 與 Pipe 的功能, 有時間再來寫程式並撰文說明這些功能的實際的應用.

OpenCL 在 2023 相關的研讀資料 / 工具

儘管近年已經沒有密集使用 OpenCL

但最近同仁因工作而起的 OpenCL 相關實務討論還是令我相當感興趣與熱衷

雖然 OpenCL 在現今已不是各家密集發展的重點技術標準

然而在有實際應用場景需求, 但沒有更好的新標準, 儘管各家計算硬體開發商缺乏合作默契的情況之下(特別是 Mobile 平台), 繼續良好支援 OpenCL 似乎是實務上的一個共識了

這篇主要的目的是統整一下在 2023 的今日要著手 OpenCL 有哪些資料可以研讀

書籍

由於 OpenCL 3.0 標準是以 OpenCL 1.2 為基底, 因此過往的 OpenCL 圖書內容都還是相當實用. 至於談到 OpenCL 2.0 標準後新增的部份 (個人覺得現今 OpenCL runtime 逐漸完備下 device-side enqueue, subgroup 以及 pipe 都可以把玩一番), 而目前僅有一本參考書目可以研讀 - "Heterogeneous Computing with OpenCL 2.0", 這本書被標為 3rd edition 的原因是前面有著兩版談論 1.1 與 1.2 的兩個版本

標準與相關資訊

OpenCL 的制定組織統整了 OpenCL 相關的資源連結 - OpenCL Resource Guide

若要查找標準相關的資訊, 另外是相關的程式資源, 這裡算是參考資料翻找的入口

各 GPU vendor 提供的 OpenCL Programming Guide

由於各家還是有自己的 OpenCL runtime 與 compiler 實作, 因此在使用與實作 OpenCL 程式還是有各家平台在 Host 與 Device 中需要注意的一些細節. 這裡整理了各家目前提供的文件.

Intel

OpenCL Developer Guide for Intel® Processor Graphics

AMD 

ROCm OpenCL Programming Guide (這是 archive.org 在去年底的備份, 其實不清楚為何 AMD 再最新的文件中將 OpenCL 的部份移除)

Nvidia

NVIDIA OpenCL Best Practices Guide Version 1.0

NVIDIA OpenCL SDK Code Samples

Qualcomm

Qualcomm Snapdragon  Mobile Platform OpenCL General Programming and
Optimization (原名 Adreno OpenCL Programming Guide)

ARM

Mali GPU Bifrost Valhall OpenCL Developer Guide

Mali Midgard OpenCL Developer Guide

Imagination

A quick guide to writing OpenCL kernels for PowerVR Rogue GPUs (Imagination 事實上是有 OpenCL performance optimization 相關文件, 只是比較可惜的是僅會對客戶提供)

工具

與 Optimization Guide 類似, 由於 GPU hardware counter 並沒有一致性的標準與介面, 所以這部份比較可惜的是必須仰賴各家的 Offline Compiler 與 Profiler

Intel

Intel® Graphics Performance Analyzers

Ahead-of-Time Compilation for GPU

AMD

Radeon GPU Profiler

Radeon GPU Analyzer

Nvidia

NVIDIA Visual Profiler

NVIDIA CUDA Compiler (NVCC)

Qualcomm

Snapdragon Profiler

Adreno GPU SDK

ARM

Mali Graphics Analyzer

Mali Offline Compiler

Imagination

PVRTune

PVRSHADEREDITOR

SIMD Programming Guide 簡報上線

今年定下的年度目標之一是建立比較完整的 SIMD 教學資源, 原本預計撰寫三個部份, 最後完成是有四個部份:

• Part 1 - 透過 compiler 跨硬體架構平台的 compiler vector extension 來撰寫 (architecture-independent) SIMD 程式
• Part 2 -講述一般而言 SIMD 指令有哪些類型, 開始導入 Arch-dependent Intrinsics, 並提供 NEON 與 SSE 的具體使用例子, 與 Lab2 的說明
• Part 3 - SIMD optimization 的一些基本步驟, 套用 SIMD ILP 的兩種基本面向, 並各自以 Gaussian 5x5 與 4x4 SAD 為舉例, 最後是 fuse / tiling / sliding-window / prefetch 等初步的進階技巧
• Part 4 - 講述 SIMD 的弱點也就是 SIMD 不擅長或是缺點的部份

 

  

 

Halide Introduction 2022 簡報上線

近日應公司同仁邀請分享對 Halide 做介紹
於是基於 2019 的投影片做一點整理並增加一點常被詢問的資訊

初探 clang 14.0 中的 Matrix Type

在先前 Clang 14.0 的文章中提到 Clang 14.0 新增了 Matrix Type, 出於個人的好奇, 於是就開始動手玩了一下, 這裡先說結論, Clang 14.0 中的 Matrix Type 還是開發階段, 尚不太好操作. 這裡嘗試使用 Matrix Type 的應用是 Tiled Matrix Multiplication.

比較基準的 vector type 實作

首先是先前已經以 Clang Vector Extension 所實作的 float 型別 8x8 矩陣乘法

typedef float float8 __attribute__((ext_vector_type(8)));
void gemm_vec(
    // buffer, stride
    float *a, int sa,  
    float *b, int sb,
    float *c, int sc
){
    float8 vb[8];
    for(int y = 0; y < 8; y++)
        vb[y] = *((float8*)(b + sb*y));
    }
    for(int y = 0; y < 8; y++){
        float8 vc = *((float8*)(c + sc*y));
        float8 va = *((float8*)(a + sa*y));
        for(int x = 0; x < 8; x++)
            vc += va[x]*vb[x];
        *((float8*)(c + sc*y)) = vc;
    }
}

以 matrix type 實作的過程

接著是參考 Clang Matrix Type 的說明嘗試定義 8x8 的 Matrix Type:

typedef float float8x8 __attribute__((matrix_type(8, 8)));

首先值得一提的是定義的 matrix type 所宣告的變數, 類似於 vector type 內作為 1D array, matrix type 變數能作為 2D array 的語法來存取個別數值. 這裡定義好 float8x8 型別之後, 必須宣告變數並將 2D 資料載入, 於是會發現文件中並沒有說明如何載入資料! 這部份必須參考目前的 Draft Specification, 目前 Clang 僅支援 3 個 function:

• M2 __builtin_matrix_transpose(M1 matrix)
• M __builtin_matrix_column_major_load(T *ptr, size_t row, size_t col, size_t columnStride)
• void __builtin_matrix_column_major_store(M matrix, T *ptr, size_t columnStride)

總之, 資料載入的方式有了, 那麼 column major load /store 的意義為何? 這主要是載入資料的順序, 參考 Wikipedia 關於 Row- and column-major order 頁面的圖解, 可以得知對於 column major load 的意義在於:

有點不太理解 clang 第一時間所支援的竟然不是 row major load / store. 然而使用 column major load store 的結果是所得到的 matrix 為原本所需的 transposed matrix. Anyway, 由於我們只有這個讀取方式, 所以我們先讀取原先 function 中的 a, b, c 矩陣:

float8x8 ma, mb, mc;
ma = __builtin_matrix_column_major_load(a, 8, 8, sa);
mb = __builtin_matrix_column_major_load(b, 8, 8, sb);
mc = __builtin_matrix_column_major_load(c, 8, 8, sc);

然而我們還是可以動些手腳來達成原本的 Tiled Matrix Multiplication, 首先我們原本要計算的資料為

C += A x B

然而這當中的 A, B, C, 透過 __builtin_matrix_column_major_load 所得到的已經是 transposed matrix (這裡我們用 A', B', C' 代表), 因此我們可以透過預先透過 __builtin_matrix_transpose 來轉回 A, B, C. 正確計算後能再轉回原先的 column major matrix 來存回正確的數值. 

ma = __builtin_matrix_transpose(ma);
mb = __builtin_matrix_transpose(mb);
mc = __builtin_matrix_transpose(mc);
mc += ma*mb;
mc = __builtin_matrix_transpose(mc); //transposed again
__builtin_matrix_column_major_store(mc, c, sc);

這裡使用了四次 __builtin_matrix_transpose, 僅僅是為了達到計算上正確的 C=AxB, 由於在回存時必須保持 transposed 狀態, 所以如果計算時能產生對應 C' 的資料即可, 於是我們可以這麼計算:

C'+=B'xA'

於是我們整個實作就可以簡化如下:


void gemm_vec(
    float *a, int sa,
    float *b, int sb,
    float *c, int sc
){
    float8x8 ma, mb, mc;
    ma = __builtin_matrix_column_major_load(a, 8, 8, sa);
    mb = __builtin_matrix_column_major_load(b, 8, 8, sb);
    mc = __builtin_matrix_column_major_load(c, 8, 8, sc);
    mc += mb * ma;
    __builtin_matrix_column_major_store(mc, c, sc);
}

最後是編譯時要注意, clang 必須加入 -fenable-matrix 這個參數

以這種實作搭配 AVX2 在 Intel Core i5-8350U 測試, 對於 naive, vector type 與 matrix type 來測試用於兩個 512x512 matrix multiplication 的效能

因此 single thread 來執行 naive, vector type, matrix type 三種實作的時間分別為 129.9ms, 26.0ms, 27.0ms. 因此以 matrix type 可以得到十分接近的 vector extension 的性能, 但是實作上相當簡潔易懂.

到此許多人應該會關心 clang 未來是否支援 row-major 的形式? 個人認為是會的, 首先 Intel AMX 指令集所使用的即是 row major 的 load/store. 再者是在 2020 年 LLVM 開發者大會中 "Matrix Support in LLVM and Clang" 的演講的倒數第二張投影片 "Remaining Work" 明確標示了 Row-major support.

使用後個人想提的幾個 idea:

定義 matrix 的同時, 應同時定義出 row 與 column 的 vector 型別, 目前不知如何取得定義語法, 這裡暫時以 ".row_vec" 與 ".col_vec" 來表示. 也就是說當我們定義了:

typedef float float8x8 __attribute__((matrix_type(8, 8)));

那麼可以用此定義一併定義出對應用在 row 與 col 的 vector type (或是分開定義但是 compiler 偵測 vector length):

float8x8.row_vec vrow, vrow2;
float8x8.col_vec vcol, vcol2;

如此可以合併使用 vector extension 與 matrix type 來處理 row major 或是 column major 的數值設定:

float8x8 mat;
...
vrow = *((float8x8.row_vec*)ptr_row);
vcol = *((float8x8.col_vec*)ptr_col);
mat.row[1] = vrow;        // set a row
mat.col[2] = vcol;        // set a column

另外可以做 row vector 與  column vector 與矩陣相乘運算:

float8x8 mat;
...
vrow = *((float8x8.row_vec*)ptr_row);
vcol = *((float8x8.col_vec*)ptr_col);
vrow2 = vcol*mat;
vcol2 = mat*vrow;

 

Halide 實務心得 5

在實作影像處理演算的過程或是實務的功能整合中, 都圍繞著 Halide::Runtime:: Buffer 的使用. 有時難以避免需要初始化預設值或是從其他 Buffer 物件複製資料. 這些可以透過 Buffer::fill() 或是 Buffer::copy_from 來做基本的填值處理.

對於一些相當 routine 但是卻又需要計算的部份,  通常會透過 loop 來走過每個座標, 甚至每個一個 value. 像是做簡單的 gamma correction 可能就會做這樣的動作:

Buffer<uint8_t, 2> src = Buffer<uint8_t>(w, h);
...
Buffer<uint8_t, 2> dst = Buffer<uint8_t>(w, h);
for(int y = 0; y < dst.height(); y++){
    for(int x = 0; x < dst.width(); x++){
        for(int c = 0; c < dst.channel(); c++){
            dst(x, y, c) = gamma(src(x, y, c));
        }
    }
}

或是可能是格式轉換 (像是 RGB <-> YUV 或是這裡的 Demosaic):

Buffer<uint8_t, 2> src = Buffer<uint8_t>(w, h);
...
Buffer<uint8_t, 2> dst = Buffer<uint8_t>(w, h);
for(int y = 0; y < dst.height(); y++){
    for(int x = 0; x < dst.width(); x++){
        uint8_t r, g, b;
        demosaic(src(x, y), r, g, b);
        dst(x, y, 0) = r;
        dst(x, y, 1) = g;
        dst(x, y, 2) = b;
    }
}

一旦數量很多, 或是 dimension 很大一一撰寫 loop 除了費時費力外, 也讓 code 不易維護與可讀性不佳. 對於這種問題, Halide 透過了 C++ lambda 的方式能夠簡短地處理, Halide Buffer 提供了兩個在其官方教學沒有提到的兩個實用的方式: for_each_value 與 for_each_element :

對於 Gamma 的例子可以簡短地改寫為:

Buffer<uint8_t, 2> dst = Buffer<uint8_t>(w, h);
dst.for_each_value([&](uint8_t &dval, uint8_t &sval){
    dval = gamma(sval);
}, src);

而對於像 demosaic 的例子則可以改寫為:

Buffer<uint8_t, 2> dst = Buffer<uint8_t>(w, h);
dst.for_each_element([&](int x, int y){
        uint8_t r, g, b;
        demosaic(src(x, y), r, g, b);
        dst(x, y, 0) = r;
        dst(x, y, 1) = g;
        dst(x, y, 2) = b;
});

 這樣的作法也適合套用在一些特殊的初始化.

Halide 實務心得 4

這陣子因為同事在弄 Hexagon DSP 且需要實作 Laplacian Pyramid 來完成所需功能, 由於這是 Halide 當初推出時 paper 內說明的範例, 所以搜尋了一下檔名, 建議他參考 Halide git repo 內的 Laplacian 範例.

今日與演算開發者開會討論時, 這位同事建議演算開發者考慮去使用 Halide language, 並且以當中的 downsample 為範例解說給演算開發的同事聽. 這時我注意到了 downsample 中有著先前沒看過, 另外也沒在任何 Halide 文件 / 教學文件提到的東西:

using Halide::_;

最有趣的是透過使用 "_" 所實作的 downsample Funciton:

Func downsample(Func f)
{
    using Halide::_;
    Func downx, downy;
    downx(x, y, _) = (
        f(2 * x - 1, y, _) +
        3.0f * (f(2 * x, y, _) + f(2 * x + 1, y, _)) +
        f(2 * x + 2, y, _)
    ) / 8.0f;
    downy(x, y, _) = (
        downx(x, 2 * y - 1, _) +
        3.0f * (downx(x, 2 * y, _) + downx(x, 2 * y + 1, _)) +
        downx(x, 2 * y + 2, _)
    ) / 8.0f;
    return downy;
}

很明顯地, "Halide::_" 是用來處理 function 內 argument / dimension 不定數量的問題. 然而詳細的使用應該有正式的解說, 經過搜尋後發現在 Halide 文件中的 Var 當中, 這個功能稱為 Implicit Variable Constructor. 文件中的解釋為: "Implicit variables are injected automatically into a function call if the number of arguments to the function are fewer than its dimensionality and a placeholder ("_") appears in its argument list. Defining a function to equal an expression containing implicit variables similarly appends those implicit variables, in the same order, to the left-hand-side of the definition where the placeholder ('_') appears." 這也就是說在 function 呼叫若使用了 "_" 在參數 list 中 , 一旦參數的數目少於函數所需參數, 將會自動對應遞補足. 透過使用 implicit variable "_" 來定義一個函數, 這表示函數有著以 "_" 出現所代表的參數 list. 

透過 implicit variable constructor 如此可以無論傳進的 Func f 有幾個參數, 像是對於影像處理的函數而言, 最前面兩個參數能固定被使用作為座標的 (x, y), 而對應的維度在處理時將直接補足對應. 因此無論純 2D 黑白影像可以套用外, 3D (有3個參數, 分別為 Width, Height, Channels) 的 RGB / RGBA / YUV444 等等影像都可以直接套用這個 downsample 來處理, 更重要的是參數更多的更高維度資料都可以因此而透過單一實作來處理. 另外應可以搭配 RDom / RVar 來使用, 像是上述程式碼中的 downx / downy 的範圍, 搭配 RDom 來使用, 應該能實作出更為簡潔的 Func 表示方式.

Halide 實務心得 3

在 Halide 的使用上會有錯覺地認為 Halide::Runtime::Buffer 的使用必須與 libHalide.so or libHalide.a linking 才可以. 但其實 Halide::Runtime::Buffer 是可以單獨使用的, 只需要 header files, 基本上一般的程式只要確認有:

#include <HalideBuffer.h>

並且確認編譯時有在 include path 即可, 單獨的使用上並不需要任何 libHalide, 相較使用 raw pointer 或是 STL 元件, 這樣的方式有很多的好處:

1. 俱備 reference count 特性
   
2. 提供類似函數的存取介面
   
3. 能夠彈性的配置 planar 或 interleaved 資料
4. 能夠帶入或取得 raw pointer
5. 可搭配 halide_image_io.h 來載入/另存 jpeg or png 圖檔
   

俱備 reference count 特性

現今的 C++ 元件多半具有這樣的特性, 除了使用上的簡潔彈性外, 好處是可以減少程式中充斥自行管理 buffer/object 管理上的配置與釋放相關的程式碼, 除了有效避免 memory leak / double free 之外, 也能降低因 pointer 操作的錯誤發生. 而這樣的特性也讓程式中對於 Buffer 的指定與傳遞更為便利.

提供類似函數的存取介面

Halide::Runtime::Buffer 並不是單純地提供 Halide 所使用的 Buffer object 存在, 每個所配置產生的 Buffer object 是能夠做存取操作的, 與 STL 不同的地方是, Halide::Runtime::Buffer 使用的是類似函式呼叫的形式來存取資料, 假設我們如下宣告了一個 Buffer.

Buffer<uint8_t> rgb_planes(1920, 1080, 3);

概念上,  這配置了 Width - 1920, Height - 1080, Channel - 3 的 Buffer, 也就是 sizeof(uint8_t) x 1920x1080x3 大小的記憶體空間, 每個 channel 有 1920x1080 , 而要存取該特定的位置, 可以使用下列方式:

// write to a pixel of a channel

rgb_plane(100, 100, 0) = 128

// read a pixel of a channel

uint_8 pix_val = rgb_plane(200, 200, 1)

可以減少計算 buffer offset 的錯誤, 程式碼也更為簡潔, 而因為採用了類似函數的形式, 程式的邏輯過程會更接近 functional 的感覺. 另外是也方便將查表轉為計算實作的驗證流程.

能夠彈性的配置 planar 或 interleaved 資料

許多的圖形資料並非單存以基本的 planar 形式存在, 難免會需要操作 interleaved RGB or RGBA 的資料, 以 Halide 參數順序代表資料格式與 loop 巢狀結構的概念, 那麼要宣告先前 rgb_plane 的 interleaved 版本, 必須這麼宣告:

Buffer<uint8_t> rgb_planes(3, 1920, 1080);

這的確是可行的, 但操作 planar 與 interleave 的程式流程就會使用不同的 index 次序, 而為了解決這樣的方式 Halide::Runtime::Buffer 提供了 make_interleaved 的介面來建構這樣的 Buffer object. 而最大的優點為, 這能夠保有與上一點相同次序的存取方式, 而對應底層不同的 memory layout.

能夠帶入或取得 raw pointer

為了銜接一些其他的實作或操作, 像是使用 OpenCV library 或是以 SIMD 加速的程式, 在 Buffer 的宣告上是能夠傳入外部的 pointer. 相對地需要 Buffer 對應的 buffer pointer 是能夠透過 Halide::Runtime::Buffer::data() 這個呼叫取得. 如此就能銜接/整合/驗證以不同方式實作的資料處理功能.

可搭配 halide_image_io.h 來載入/另存 jpeg or png 圖檔

在 Halide tutorial 中使用的 load_image 與 save_image 對於實務上是很便利的工具,  這也能在 include halide_image_io.h 之後個別來使用, 而需要注意的是這會需要 link libpng 與 libjpeg.

Raspberry Pi Zero W Project Part 3 - Software Interrupt & Handler example

For ARM architecture, the instruction for software interrupt is SWI or SVC (ARMv7).

Using the instruction is very simple:

SVC #imm

When an ARM processor executes the instruction, it will switch to SVC/SWI mode and jump to the SWI or SVC hander in exception vector table specified by Vector Base Address Register (VBAR).

Similar to previous parts, please get from repo https://github.com/champyen/rpiz_bare_metal.git 

please checkout the commit: 637086c2

$ git checkout 637086c2

It is not difficult to use SWI/SVC instruction. But it requires experience to apply it well. Using SWI/SVC instrution for system calls, we need to:

• design system call table ( SWI number , system call pair)
  
• implement SWI/SVC handler, get the system call number
  
• implement system calls with corresponding number
  

For get the SWI (or system call) number, it is describe in ARM's "SWI Handler's" document. The SWI number is encoded as lowest 24 bit in SWI/SVC instruction. Therefore we could get SWI/SVC number by fetching the instruction itself and clearing the MSb 8bits with BIC instruction right after entering SWI hander:

ldr     r0, [lr,#-4]
bic     r0, r0, #0xff000000

After getting SWI number, we could pass it to system_handler as an argument saved in R0. The system_hander is very similar to ISR, but we don't need to check the hardware status to know which hardware interrupts CPU. We just need to handler the SWI or system call by the SWI number. Of course, as IRQ has Interrupt latency SWI/SVC instruction consumes more cycles than normal instructions for most CPUs.

For Application processor it is easy to understand the meaning of mode switching and SWI handling. In fact, Cortex-R, and Cortex-M processor also has SWI/SVC instruction. For most Embedded / RTOS developers, they think it is useless or trivial in such processor.

For some embedded or RTOS, mutual-exclusive applications can be loaded on demand. This can be done by well-designed linker-script (e.g.: different LAs with same). For such Embedded / RTOS, SWI / SVC instruction is useful to maintain an API layer between Applications and Kernel. An intuitive / naive way is to setup and maintain table of function pointers on the kernel side (Of course we still need to setup API index for using the function as SWI number). For development it has an drawback: it is hard to maintain or adjust or expand the table. With SWI/SVC instruction, no table forwarding is needed. Besides it provides flexibility to group system calls and reserve numbers for future needs. (Since for table-based method, it required to reserve a huge table).

Raspberry Pi Zero W Project Part 2 - Interrupt Service Routine (ISR) example

After implementing naive function and printf, let's add ISR to it.

please checkout the commit: 91abc0d3

$ git checkout 91abc0d3

There are many changes from previous commit:

head.S

it becomes more complicated. As we know, its context starts at 0x8000 address.

start:
    ldr     pc, reset_target        /* 0x00 mode: svc */
    ldr     pc, undefined_target    /* 0x04 mode: ? */
    ldr     pc, swi_target          /* 0x08 mode: svc */
    ldr     pc, prefetch_target     /* 0x0c mode: abort */
    ldr     pc, abort_target        /* 0x10 mode: abort */
    ldr     pc, unused_target       /* 0x14 unused */
    ldr     pc, irq_target          /* 0x18 mode: irq */
    ldr     pc, fiq_target          /* 0x1c mode: fiq */

reset_target:           .word   reset_entry
undefined_target:       .word   undefined_entry
swi_target:             .word   syscall_entry
prefetch_target:        .word   prefetch_entry
abort_target:           .word   abort_entry
unused_target:          .word   unused_entry
irq_target:             .word   irq_entry
fiq_target:             .word   fiq_entry

After loading the binary, it will jump to a routine named "reset_entry". Before we trace the reset_entry. The 8 "ldr pc, XXXXXX" instructions are so called Exception Vector Table. (FIQ is a special irq mode, it has a advantage - the implementation can be start at the location, the jump is not necessary. Therefore one jump delay is saved. ) It is used to handle system exceptions. Each has corresponding privileged mode to it. Besides, each mode has dedicated LR and SP registers - this means OS / firmware implementation should take care of stack space arrangement for the mode:

In fact, for reset_entry here, its major work is setting stack for each mode:

reset_entry:
    /* set VBAR to 0x8000 */
    mov r0, #0x8000
    mcr p15, 0, r0, c12, c0, 0

    /* (PSR_FIQ_MODE|PSR_FIQ_DIS|PSR_IRQ_DIS) */
    mov r0,#0xD1
    msr cpsr_c,r0
    ldr sp, stack_fiq_top

... other 4 modes ...

    /* (PSR_SVC_MODE|PSR_FIQ_DIS|PSR_IRQ_DIS) */
    mov r0,#0xD3
    msr cpsr_c,r0
    ldr sp, stack_svc_top

    cpsie i
    bl  bare_metal_start

In addition to stack assignment and jump to bare_metal_start , there are two key points here:

1. setup Vector Base Address Register (VBAR) - From ARMv6, the exception vector can be placed other than 0x00000000 and 0xFF000000. This is achieved by setting VBAR, please refer to "3.2.43 c12, Secure or Non-secure Vector Base Address Register" in ARM1176JZF-S TRM.
   
2. enable interrupt - cpsie instruction
   

And we have to trace isr_entry:

irq_entry:
    stmfd  sp!, {r0-r12, lr}
    add     lr, pc, #4
    bl      isr_entry
    ldmfd   sp!, {r0-r12, lr}

    subs    pc, lr, #4

For ISR, it is not surprised to backup and restore all (non-dedicated) registers. The most interesting things are - LR register setting and the instruction to leave IRQ mode. For LR setting it is easy to figure out, the target return address is the 'bl isr_entry' not 'add lr, pc, #4". That's the main reason to save "pc+4" to LR. And for leaving each mode, please refer to "2.12.2 Exception entry and exit summary" of ARM1176JZF-S TRM.

isr.c

There 3 functions in the source file: timer_enable, timer_check and isr_enty. Here we use "System Timer" in BCM2835, please refer to Chap 7 and Chap 12 of "BCM2835 Peripheral specification". Besides the IRQ number of System Timer is not listed in the document, please refer to the link of  "errata and some additional information" on the page.

The timer_enable will enable System Timer 1 or 3 by index and timer_check is used to clear IRQ state and update next timeout interrupt. Therefore the isr_enty just check status and call timer_check for clear the IRQ.

bare_metal.c

For demonstrate IRQ and main thread's progress, a busy loop with counter is added. The loop will print out a number when specified condition is met. And Timer is enabled before the loop, You can see the timer tick with ISR and the main thread keeps counting.

void bare_metal_start(void)
{
    int base = 0;
    asm volatile (
        "mov %0, sp\n\t" : "=r" (base)
    );

    printf("\n\n%s:%x: Hello World! %s %s %d\n\n", __func__, base, __DATE__, __TIME__, __LINE__);
    printf("enter busy loop\n");

    timer_enable(1);

    volatile int i = 0;
    while(1){
        if((i++ & 0x00FFFFFF) == 0)
            printf("%d\n", i);
    }
}

 

 

 

 

 

Linux Developer 都應該學習的 Container 工具 - Apptainer / Singularity

現今對於 Daily Linux Developer / User 面對不同程式/開發版本環境感到很頭疼, 常常疲於

• 執行舊版程式需要安裝舊版本 Library, 設定 RPATH / LD_LIBRARY_PATH 
  
• 開發需求建立不同的版本 SDK 開發/執行環境, 在較舊系統上需要新軟體, 或是新系統上需要舊軟體工具(基於驗證或是綁定版本等等需求)
• 短暫需要的工具, 安裝了一堆套件一一全部移除相當麻煩
• 嘗鮮希望測試實驗性版本軟體

面對這些林林總總開發上的需求, 但是 ... 卻一定不希望弄亂/髒平日使用的主要系統環境

而在接觸到 Singularity 之前

個人認為 python 的 Virtual Environment (VENV) 是很酷的東西

因為可以很快的在同路徑與環境下直接切換不同版本的開發環境

而儘管 docker / VM 也能夠提供多種版本的開發環境

但是 docker / VM 也著一些相對不便的缺點: 

• 需要 root 才會有比較方便的使用體驗
• 需要登入 / 切換使用者環境, 不是開發者既有的 Home
• isolation 對 service 是優點, 但對開發完全成為一個重大缺點
  
• 目錄資料需要複製 / 轉移 / 同步
  

本文談到的 rootless 與環境功能相信 docker 都能提供, 但是流程相對煩雜非常非常多

相關工具也並非官方提供, 會有版本造成無法運作的問題

而透過使用 Singularity 可以很輕鬆達成完整環境類似 VENV 的效果

Singularity  的好處很多, 對於開發者而言:

• 是 container 技術, 不像 VM 有效能減損問題
  
• 基本的設定與使用非常簡單
  
• 能直接在既有的 Home 目錄下切換環境, 無同步上的困擾
• 可使用 docker hub 的資源
• 簡單的 image / sandbox 概念設計
• 透過 fakeroot 就可以建立/修改自己需要的 image 環境, 無需 root
• 符合 OCI 標準,  能提供 Docker 那類 service 導向的 instance 環境
  

Singularity 因為無需 root 權限的特性廣受 HPC 相關應用服務的歡迎與採用

但是 Singularity 的文件數量並不像 Docker 那樣豐富, 但還是有相當不錯的介紹

個人推薦瑞士University of Bern 中 Science IT 單位提供的 Singularity Container 指引

當然還是推薦官方網站與官方文件中的 User Guide 與 Admin Guide

本文單純是以 Developer 角度的推薦說明文

包含安裝過程與詳細的操作與進階功能請參考上述文件連結

建立基本 SIF image, 基本操作

安裝 Singularity 後, 當然立刻會想要抓個環境來試試看

$ singularity build ubuntu_18.04.sif docker://ubuntu:18.04

接著你會看到類似下列的畫面

大小約 26MB 左右, 可想而知環境是相當陽春的, 但接著嘗試用看看吧

$ singularity shell ubuntu_18.04.sif

從截圖中可以看到主系統是 Ubuntu 20.10 但 SIF 環境是 18.04, 此外還可以比較一些工具

完全是不同的世界! 更棒的是第一張截圖顯示了還是在同一個路徑下

這時候好奇嘗試 apt-get update,  或是 sudo apt-get update 會發現完全無法執行

那麼該怎麼安裝想要的套件呢?

這時候你需要建立 sandbox (也就是 image 的 rootfs)

再透過調整好的 sandbox 來建立 sif image

建立 Sandbox

在 University of Bern 的投影片中絕大多數需要 sudo 的操作都可以透過 --fakeroot 參數取代

像是建立 sandbox 原本的指令會是:

$ sudo singularity build --sandbox ubuntu_18.04 docker://ubuntu:18.04

但是透過 fakeroot 特性, 現在只需要:

 $ singularity build --fakeroot --sandbox ubuntu_18.04 docker://ubuntu:18.04

如此你就建立了可以修改的 sandbox, sandbox 的功能與 SIF image 相同

但最大的不同是它可以 writable! 所以可以用來修改與建立想要的環境

接著讓我們以安裝 gcc 與 clang 為例

首先透過 fakeroot 以 writable 方式開啟 ubuntu_18.04 環境, 並 apt-get update

$ singularity shell --fakeroot --writable ubuntu_18.04

接著來以 apt-get install gcc clang 指令試試看安裝 gcc & clang

要安裝的是 gcc-7 與 clang-6.0 完全就是  ubuntu 18.04 的套件版本

透過 sandbox 建立 SIF image 檔案

安裝完後當然可以直接使用 gcc / clang, 但這裡我們先離開 sandbox 環境

接著要把修改好的 sandbox 封裝成方便轉移與可供他人使用的 SIF image 檔案

將 sandbox 轉為 sif 的指令為:

$ singularity build --fakeroot ubuntu_18.04_dev.sif ubuntu_18.04

是的, build 的來源不僅只能是 docker hub, 還能由客製的 sandbox 與 singularity hub

在此觀察了 sif 檔案大小, 發現大小增加到 210MB 上下

這個新建立的 ubuntu_18.04_dev.sif 就可以提供給任何人來使用

編譯測試

這裡建立了一個 test/test.c 作為切換的展示

而 test/test.c 內容為:

#include <stdio.h>

int main(void)
{
    printf("This is a singularity test.\n");
}

這個測試是展示在相同的目錄下切換不同系統環境

清除 sandbox 目錄

由於 sandbox 是透過 fakeroot 身份建立的, 因此可能會有無法用自身帳號移除部份檔案/目錄的情況

這時候請執行 fakeroot 指令轉換為 (fake) root

雖然顯示為 root, 但是事實上並不是真的 root 也無法執行實際系統上需要 root 的指令 (這裡嘗試執行 apt-get update)

這時用此方式就可以順利移除產生的 sandbox

如何? Singularity 是不是簡單好吃不黏牙?

UPDATE - 由於 Sylabs 的公司化與透過 Singularity 與營利, 然而又並沒有改名, 因此原本社群版本在 2021 年 11 月時更名為 Apptainer

關於 clang/gcc vector extension 的 unaligned data load / store

從 ARM 平台開始學習著手寫 SIMD intrinsics 有個相當大的好處 - 在 ARM 平台上 vector load/store 在任何 memory address 都可以使用 vld1q_u32 or vst1q_u32 這類的 intrinsics (另外像是 vld3q_u8 / vst3q_u8 對於 interleaving RGB 非常好用, 在 x86 上就很頭大, 但這是題外話了), 也就是說 ARM 的 vector load / store 可以不用考慮 vector width alignment. 

然而 ARM NEON 的 vector load / store 指令設計是少數, 在多數的 SIMD ISA 平台都不是如此, 像是 Intel AVX / SSE 就必須顧慮到 alignment ( 請參考 Intel Intrinsics Guide 網站 ):
以 SSE 而言:

• _mm_load_si128 - 這是用來針對 aligned address
• _mm_loadu_si128 - 這是使用來處理 unaligned address

對於 AVX :

• _mm256_load_si256 - 這是用來針對 aligned address
• _mm256_loadu_si256 - 這是使用來處理 unaligned address

在 unaligned address 使用 aligned instruction 就會發生 crash / error, 因此不可不慎, 而這在硬體上會反映在內部 pipeline 狀態到 memory request 的頻率/數量, 對於效能有一定程度的影響. 在多數 DSP 上由於使用 banked SRAM, 因此有可能有更嚴格的 alignment 限制. 

在 clang / gcc 中使用 vector extension 是可以較為簡單撰寫 SIMD-friendly code 的方法之一, 然而無論在 clang 的 vector extension 說明 或是 gcc 的 vector extension 說明 都採用 aligned vector 的形式, 所以通常在撰寫 data loading 時就必須注意 pointer alignment:

// clang vector type
typedef float float4 __attribute__((ext_vector_type(4)));
...
float4 *vdata = *((float4*)(data_ptr + offset));

對於像是 CV / Image Processing 通常必須注意 padding, boundary 與 data_ptr 的 offset, 最簡單的方式當然是透過 memalign 之類的方式配置 buffer, 然後將 layout 設計好符合 SIMD width. 然而並非所有演算都可以如此輕鬆, 一些像是 connected conponents 之類的演算, 能夠有一定程度的 SIMD, 但是又並不保證起點由 aligned address 處理, 這時通常處理方式有二:

1. 使用 unaligned read + indications/select
   
2. 轉為 aligned-based loops

兩者都有人使用, 這裡要談的不是 algorithm 層級的事情, 且 1. 的部份會簡單一些, 那麼問題就在於該如何讓 compiler 產生 unaligned load / store. 對於 gcc / clang vector extension 而言這在於 type alignment 屬性 (詳情請參考 StackOverflow 對此問題的回覆) , 將上述的 float4 型別宣告改為:

typedef float float4 __attribute__((ext_vector_type(4), aligned(1)));

如此 float4 型別的 alignment 會是 1-byte, 如此 compiler 會轉為產生 unaligned load / store 指令, 便能在 c code 當中於任意位置做資料的存取與處理.

以 halide 實作 block matching 演算法 - block matching with Halide language

由於 multi-frame 與 multi-view 的影像處理需求的提升,  image alignment 的需求日益提升, 而相關的演算也很多, 像是基於 pyramid-based, 9-steps, hexagon or diamond search. 而大小的部份也有著 match / apply 不同大小的 non-local means, 然而這些演算的基本還是源自於特定 search window 的 block matching. 本文重點在於以 Halide language 描述 matching algorithm, 至於 scheduling 的優化不在此文重點. (或是擇日分享)

Step 0 - 輸入格式 (input format - single channel as example)
首先假設要處理的是 image 的某個 channel, multi-channel 的影像則僅是套用到多個 frame 或是以 channel 輸入, 同常有一張是 base image 這裡以 in 代表, 另一張為 reference image 以 ref 代表.

Step 1 - 邊界處理, 用來避免 search block 範圍超過邊界 (boundary handling, to avoid error of part of target block is out of frame)

Func input = BoundaryConditions::repeat_edge(in);
Func refer = BoundaryConditions::repeat_edge(ref);

Step 2 - 用以代表 block 內所有 input 與 reference 的 pixel 的座標 (setup the RDom that represent all pixels of corresponding input & refer block)

// the position in motion vector output table

Var bid_x, bid_y;

// variable for motion vector

Var mv_x, mv_y;

RDom rblk(0, BLOCK_SIZE, 0, BLOCK_SIZE)
// position in input frame
Expr blk_x = (bid_x * BLOCK_SIZE) + rblk.x;

Expr blk_y = (bid_y * BLOCK_SIZE) + rblk.y;
// position in reference frame

Expr ref_x = blk_x + mv_x;
Expr ref_y = blk_y + mv_y;

Step 3 - 建立比對標準, 這裡以 SAD 為例 ( matching criteria calculation (Sum of Absolute Difference(SAD) as example))

Expr in_val = i16(input(blk_x, blk_y);

Expr ref_val = i16(refer(ref_x, ref_y));

Func sad;

sad(mv_x, mv_y, bid_x, bid_y) = sum(abs(in_val - ref_val));

Step 4 - 搜尋與輸出, 這裡是透過 RDom + argmin 搭配來做比對搜尋並取得結果 (search and output, here we use RDom + argmin for matching and finding)

// setup search window
RDom search_win(-range/2, range, -range/2, range);

// argmin will find the one in search_window which get minimal sad.

Tuple me = argmin(sad(search_win.x, search_win.y, bid_x, bid_y));

// output the motion estimation result to MV table

Var ch;

out_mv(bid_x, bid_y, ch) = me[ch];

比較不直覺的地方在於兩個 RDom 套用的層次不同, 第一個 RDom 用來計算 matching criteria, 二個 RDom 是用來作為 motion estimation 尋找 mv, 總之這是以 search window 的方式作窮舉的 block matching 定義, 然而可以延伸與改進用在一開始所說的各種 matching 演算法上

撰寫各 x86 CPU 平台適用的 compiler vector extension code

在套用 gcc / clang 的 -ftree-vectorize or compiler vector extension 時來產生適當的 SIMD optimization 時,  同時會產生一個問題 - 當編譯時套用的 architecture 參數(e.g.: -mavx2) 所使用的指令集與運作時的 CPU 不一致時可能會出現 SIGILL (illegal instruction) 的問題. 過往可能會採用 CPU feature detection + dispatching code 的撰寫, 大致上邏輯如下:

// AVX-512 version

void __foo_avx512(void){  ... }

// AVX version

void __foo_avx2(void){ ... }

// SSE4 version

void __foo_sse4(void){  ... }

// C version

void __foo(void){  ... }

void foo(void)

{

    if(_avx512_available())

        __foo_avx512();

    else if(_avx2_available())

        __foo_avx2();

    else if(_sse4_available())

        __foo_sse4();

    else

        __foo();

}

這是過去長久以來很平常的作法, 而最近研究在接觸 Clear Linux, 研究了一下 Intel 採用的方式, 發現了有趣的事(這些 Intel 都有在 2016 年的簡報 - Improving Linux Performance with GCC latest technologies 中說明), 主要是 gcc 4.8 中已新增了 Function Multi-Versioning (FMV)的功能, 能透過函數屬性更方便地做到原本的目的, 以上述的例子來說就會變成:

// AVX-512 version

__attribute__ ((target ("avx512f")))

void foo(void){  ... }

// AVX2 version

__attribute__ ((target ("avx2")))

void foo(void){ ... }

// SSE4 version

__attribute__ ((target ("sse4.2")))

void foo(void){  ... }

// C version

__attribute__ ((target ("default")))

void foo(void){  ... }

使用 FMV 功能的好處首先所有的 function 都是相同的名稱, 僅是屬性中的 target 不同, 而GCC 會自動產生偵測 CPU 功能並從中選擇對應 foo 的 code, 這點可以省去維護上面 dispatching code 的問題(若有多個函數,  增刪都會是很冗長的工作). 使用 __attribute__ ((target ("TARGET_NAME"))) 的語法方式目前已被 gcc / icc / clang 所採用.

然而並不是所有人都有心力以 SIMD ISA 去實作版本, 這時複製多份相同的 C code 來做不同 SIMD ISA 的 auto-vectorization 似乎不是很實際的方式, 因此 gcc 進一步提供了有趣的屬性功能 - target_clones:

// C version

__attribute__ ((target_clones ("avx512f", "avx2", "sse4.2", "default")))

void foo(void){  ... }

而 GCC FMV 功能亦能夠與 compiler vector extension 合併使用, 因此以 compiler vector 實作後函數可以使用上述 target_clones 的方式, 如此 compiler 會自動針對不同的 SIMD ISA 生出多版本的, 以個人常用來教學的 tiled 8x8 gemm 來說就會像是:

#define TSIZE 8
#if defined (__clang__)
typedef float vfloat __attribute__((ext_vector_type(TSIZE)));
#else
typedef float vfloat __attribute__ ((vector_size (TSIZE*4)));
#endif

__attribute__ ((target_clones ("avx512f", "avx2", "avx", "sse4.2", "default")))
void gemm_vec(float *a, int sa, float *b, int sb, float *c, int sc)
{
    vfloat vb[TSIZE];
        for(int y = 0; y < TSIZE; y++){
            vb[y] = *((vfloat*)(b + sb*y));
        }
        for(int y = 0; y < TSIZE; y++){
            vfloat vc = *((vfloat*)(c + sc*y));
            vfloat va = *((vfloat*)(a + sa*y));
                for(int x = 0; x < TSIZE; x++){
            vc += va[x] * vb[x];
                }
                *((vfloat*)(c + sc*y)) = vc;
        }
}

接著以 gcc 編譯:

$ gcc -O3 -shared mm.c -o libgemm.so

透過 objdump 觀察

$ objdump -t libgemm.so

可以看到中間有如下的輸出:

0000000000001720 l     F .text    0000000000000c96 gemm_vec.default.4
00000000000023c0 l     F .text    00000000000005df gemm_vec.avx512f.0
00000000000029a0 l     F .text    00000000000005df gemm_vec.avx2.1
0000000000002f80 l     F .text    000000000000060c gemm_vec.avx.2
0000000000003590 l     F .text    0000000000000c95 gemm_vec.sse4_2.3

此外值得一提的是在 2019 年初 Clear Linux 計劃中 Intel 釋出了一個 FMV patch generator, 透過 make-fmv-patch 這個工具能針對  C/C++ code 自動產生對應的 FMV 屬性 patch, 套用後重新編譯即可產生一體適用的 SIMD optimized libraries / executables

clang & gcc vector extension 探究心得

這文章換整了把一個以前納悶許久的功能做了研究與釐清後的心得 - vector selection in GCC / Clang

在 GCC / Clang 的 vector extension 中 "?:" operation 的使用, 基本上搜尋 "gcc vector extension" or "clang vector extension" 你一定會看到文件對於 ?: 這操作寫著是支援的 甚至 clang 的文件 這麼列著.

 

 

然而按照常用的 CL 使用方式去寫 code 會發現沒有辦法編譯過.
對於 OpenCL 中的使用, 請參考 Programming with OpenCL C 一書的 Vector Operators - Conditional Operator

對於 gcc/clang 支援的 ?: 是指 a ? b : c 這樣的計算中
的 a 必須是 scalar type, b, c 可以是 scalar or vector ... 因此如此是無法達到最好用的 vector element selection 功能, 也就是當我們撰寫 c = a > b ? a : b; 時於 vector 上等同於:

u32x8 c = {
    a[0] > b[0] ? a[0] : b[0],
    a[1] > b[1] ? a[1] : b[1],
    a[2] > b[2] ? a[2] : b[2],
    a[3] > b[3] ? a[3] : b[3],
    a[4] > b[4] ? a[4] : b[4],
    a[5] > b[5] ? a[5] : b[5],
    a[6] > b[6] ? a[6] : b[6],
    a[7] > b[7] ? a[7] : b[7],
};

以這個例子來說當宣告:

u32x8 a, b, c;
...
c = a > b ? a : b;

對於 x86 AVX ISA 而言你會希望 compiler 能使用下列 intrinsics / instruction:
_mm256_max_epu32 / vpmaxud

但是很不幸的是在 gcc / clang 這樣的 code 無法被編譯, 在不考慮直接使用 AVX intrincs 的情況下, 搜尋後會得到兩個可能的建議

1. 使用 loop 讓 compiler 優化, 也就是:**

for(int i = 0; i < 8; i++)
    c[i] = a[i] > b[i] ? a[i] : b[i];

如此你會得到如下圖的使用 loop 結果, 雖然有使用 AVX register, 但是計算上完全沒有好處

 

 

2. 這個結果比較需要技巧才會找到, 也就是 bitwise operation:**

u32x8 d = (u32x8)(a > b);
u32x8 c = (d & a) | (~d & b);

如此你會得到圖三的結果, 以 code generation 的結果這可能是相對較好的方式

 

然而這完全比不上直接使用 avx intrinsics 的下圖的結果:

u32x8 c = _mm256_max_epu32(a, b);

g++/clang++

然而還是有最後希望 - GCC/g++, 事實上 gcc 文件中如下圖有提到 C++ 模式支援 element-wise

由於個人印象中並沒有這樣的方式, 查了一下 是 GCC 5 / Clang 10 開始支援這樣的模式 (只能說在最需要使用的時候 android gcc, Google 只用到了 4.9 ... 之後的版本 vector extension 沒特別去用, 加上逐漸重心轉到使用 clang), 如此可以支援一開始期望的語法, 而且 code generation是很漂亮的(g++:下圖1/clang++:下圖2, g++ 的比較好), 非 clang 不可時透過 g++/clang++ 編譯產生 .s 再整入 clang project 是唯一可以考慮的作法 (否則就是 code 可讀性與移植性變差 + 努力查 intrinsics table 了)

1. 

 

2.

 

除此之後後續做了些實驗額外得知了一些限制/資訊

gcc & clang compatibility

在先前 gcc / clang 各自有著自己的 vector extension 與宣告方式, 個人也針對兩者做了不同宣告, 然而有趣的地方是實驗過程中, 曾經忘了切換, 然而可以順利編譯過, 這也表示 gcc / clang 互相吃了對方的 vector extension 的宣告方式, 因此個人一開始以為在各自 compiler 上兩個 vector extension 應該是互通的. 這點在昨日得知 clang 在實作時會考量儘量與 gcc 一致, 猜想 gcc 在一些特性上也對 clang 做了類似的考量.

OpenCL-like vector swizzle

這是被詢問的問題, 因此注意到了差異, 基本上只有 clang 提供了接近 OpenCL 的 vector swizzle 方式:

u32x4 a, b, c;
...
c.s02 = a.s02; // or c.xz = a.xz;
c.s13 = b.s02; // or c.yw = b.xz;

當然還可以針對特定 vector lane 作 value assign.
比較麻煩的是 initial value, clang 並沒有提供如同 OpenCL 的彈性, 而且必須是以 { ... } 的方式一一擺放, 而不如 OpenCL 以 ( ... ) 且能夾雜大小不同的 vector 作為數個數值. 在 clang 上要使用這個特性, 必須使用 clang 專屬的 vector extension 宣告方式.

vector as array

最早是 clang 開始提供的, 而且不只是 vector extension, 對於 SIMD intrinsics 的 data type clang 都支援. 這個功能相當實用, 在 android 7 時期導入 clang 可以說幫了大忙, 這點對於 C Model 與 SIMD 版本實作的轉換與驗證簡化不少.
使用上類似:

u32x4 a, b, c, vout;
...
vout = fa*a + fb*b + c;
for(int i = 0; i < 4; i++){
    assert(out[ofs+i] == vout[i]);
}

而後續 gcc 也在 vector 支援了這個功能, 然而必須使用 gcc 自己的 vector extension 宣告方式