2018年4月17日 星期二

Halide's FFT 實作與效能比較

在許多 Halide 相關論文的投影片都可以看到與 FFTW 的速度比較
宣稱以 Halide 實作的 FFT 相較為 FFTW 的 1.5X 快

像是 An Introduction to Halide (CVPR 2015
以及 Halide: Decoupling Algorithms fromSchedules for High-PerformanceImage Processing

但是一時間無法找到原始碼
而最終是透過這篇 github 上的討論
由其得知現在是在 Halide/apps/fft 目錄中

找到後就簡單多了, 接著在 Ubuntu 上準備環境
$ sudo apt-get install libfftw3-3 libfftw3-dev
接著在 Halide/apps/fft 目錄下, 首先編修 ../support/Makefile.inc
$ vi ../support/Makefile.inc
修改 HALIDE_BIN_PATH, LDFLAGS 與 LLVM_CONFIG 後如下編譯與執行
$ WITH_FFTW=1 make bench_16x16 bench_32x32 bench_48x48 bench_64x64
./bin/bench_fft 16 16 ./bin/
                       Halide                  FFTW    
    DFT type  Time (us)    MFLOP/s  Time (us)    MFLOP/s      Ratio
         c2c      0.370   27659.83      0.575   17808.48       1.55
         r2c      0.269   18998.29      0.517    9910.42       1.92
         c2r      0.271   18918.61      0.518    9888.96       1.91
./bin/bench_fft 32 32 ./bin/
                       Halide                  FFTW    
    DFT type  Time (us)    MFLOP/s  Time (us)    MFLOP/s      Ratio
         c2c      3.121   16403.89      2.372   21586.78       0.76
         r2c      1.222   20947.74      2.415   10601.24       1.98
         c2r      1.372   18661.48      2.591    9878.77       1.89
./bin/bench_fft 48 48 ./bin/
                       Halide                  FFTW    
    DFT type  Time (us)    MFLOP/s  Time (us)    MFLOP/s      Ratio
         c2c     10.113   12724.26      8.103   15879.26      0.801
         r2c      3.986   16141.35      8.057    7985.28       2.02
         c2r      4.221   15242.53      7.610    8454.44        1.8
./bin/bench_fft 64 64 ./bin/
                       Halide                  FFTW    
    DFT type  Time (us)    MFLOP/s  Time (us)    MFLOP/s      Ratio
         c2c     21.800   11273.26     11.506   21358.44      0.528
         r2c      7.588   16193.94     11.638   10558.20       1.53
         c2r      7.530   16319.53     11.560   10629.47       1.54

在 m3-7y30 平台執行結果如上, 可以觀察到常使用的 r2c/c2r 的 2D FFT 效能相當優異, 也比 2015 年的數據出色

2018年4月11日 星期三

Halide - 實務心得 2

與第一篇著重在與 C/C++ 整合與比較特別的小技巧有所不同
在第二篇中談論到的是演算應用實務上會碰到值得注意的地方
由於此系列並不是有所規劃的, 因此日後累積足夠多心得再撰寫後續文章

演算法歸納撰寫方式

Halide 對於 RDom 的處理相當一板一眼
其強調於演算與排程的脫勾
無論如何必須注意的是 Halide 並不會在演算層次上分析多餘的計算並做演算變更
因此演算描述的好壞也決定了一部份效能的結果, 所以演算上的精練也是很重要的
像是影像處理上許多的演算實則為同一 filter 套用在不同方向上 (eg: gaussian blur or box filter)
若符合這個條件時則可以考慮分為 2-stage 的作法, 這樣的好處有二:

1. 大幅減少運算
像是 5-tap filter, 若套用在 2D 上需要計算 5x5 = 25 次
若分為 2-stage 則為 5+5 = 10次, 可以說 tap 數愈高, 差異愈大

2. 可重複某一平方向結果
由於是在垂直方向套用相同 filter , 因此可以將套用一個方向的結果存下, 重複利用在不同點的計算

因此在演算的撰寫上還是需要注意演算的表示上是否有進一步優化的可能

AOT 在局部影像 ROI (Region Of Interest) 的處理

特定情況需要套用演算在特定區域上, 而非全畫面的處理
這時可以嘗試產生對應特定區域的 Buffer  物件
Halide::Runtime::Buffer 中支援 cropped 操作
對於此操作文件描述為 "Make an image that refers to a sub-range of this image along the given dimension.", 明確指出是直接參考 Buffer 的特定區域.
使用範例:
Buffer<uint8_t> buf(width, height);    //2D buffer
Buffer<uint8_t> crop_buf= buf.cropped(0, crop_x, crop_w).cropped(1, crop_y, crop_h);
如此可以將取得的 crop_buf 帶入到 AOT function 中, 即可只對特定區域做處理

Generator / Auto-Scheduler

JIT 方式儘管套用相當方便,
但因第一次呼叫都需要花時間 realize, 且並不是所有平台都支援 JIT
因此使用上對於嘗試與試驗較好;
對於固定流程之後的應用而言, 使用 AOT 較為實際.
強烈建議將實作的 Func 在獨立撰寫之後轉為 Generator 形式,
除了能夠迅速跨越不同硬體架構平台,
也能夠在完成後快速地套用 Auto Scheduler 來做優化的嘗試. (甚至如同這篇一般做不同優化方式比較)

2018年4月1日 星期日

(Halide) Generator 增進 - 非官方中文翻譯

本篇是 Halide Wiki 上的一篇 Generator Enhancement 的翻譯
由於 Halide Tutorial 中使用了新的寫法, 因此需要了解其中的改進, 儘管 Tutorial 中已經提供基本的 Generator 使用方式, 然而許多進階與注意事項在此文件提到了更多
Generator 這些改進提供了更清晰與洗煉的表示方式, 新的參數讓編譯更有彈性, 以及更強大的混合使用, 新的設計更考量了 AutoScheduler 的搭配, 值得花時間學習

====

在 2016 年 10月下旬 (https://github.com/halide/Halide/pull/1523), Halide Generators 已被增進與改善:
  • 改善 Generator 的可讀性與彈性
  • 提供機械產生的 Stubs 讓 Generator 能彼此間的使用更簡單
  • 使得與 Autoscheduler 間的整合更容易與更可靠
請注意並沒有任何的改變會讓既有存在的 Generator 無法運作(所有現在既存的 Generator 應該如其地運作); 所有現存的 Generator 在可見的未來將持續如同以往地運作.
這份文件主要在於提供改變的本質與描述如何"升級" Generator 來使用新的增進功能.

使用 Input<X> 取代 Param<X> (以及用 Input<Func> 取代 ImageParam )


Param<> 依然存在, 但 Generators 現在能夠使用一個新的類別, Input<>, 來取代. 對於純量型別, 這些實質上能夠被認為等同於 Param<>, 但因為一些後續將介紹的程式碼簡明緣故而有著不同的名字.
同樣地, ImageParam 依然存在, 但 Generators 現在能夠使用一個 Input<Func> 來取代. 這(實質上)類似一個 ImageParam, 而有著最主要的差異是背後可能不(或者可能)再有著實際的 buffer, 因此沒有事先定義的延伸範圍.
Input<Func> input{"input", Float(32), 2};
類似於下列方式, 這等義於一個背後有著能被一個的 Input<Buffer<T>> 所建構的實體 buffer 的 ImageParam:
ImageParam input{UInt(8), 2, "input"};
轉變為:
Input<Buffer<uint8_t>> input{"input", 2};
這允許程式寫作者 (相較於使用 Input<Func>) 透過 input.dim(0).extend() 以及 input.dim(1).extend() 來取得 buffer 的寬與高.

同時宣告 Input<>Param<>ImageParam (i.e.: 如果使用 Input<> 可能無法使用先前的語法) 對於 Generator 是一種錯誤. 請注意 Input<> 只能搭配 Generator 使用, 並且無法使用在其他 Halide 程式碼中; 除了在 Generator 中之外, 它並非打算來取代 Param<>


範例:
class SumColumns : Generator<SumColumns> {
  ImageParam input{Float(32), 2, "input"};

  Func build() {
    RDom r(0, input.width());
    Func f;
    Var y;
    f(y) = 0.f;
    f(y) += input(r.x, y);
    return f;
  }
};
轉變為
class SumColumns : Generator<SumColumns> {
  Input<Func> input{"input", Float(32), 2};
  Input<int32_t> width{"width"};

  Func build() {
    RDom r(0, width);
    Func f;
    Var y;
    f(y) = 0.f;
    f(y) += input(r.x, y);
    return f;
  }
};
能夠選擇性地不指定 Input<Func> 的型別與/或維度, 在這時數值直接由 Func 傳入的數值推斷. 當然, 若指定了明確的型別與維度, 將要求傳入的 Func 符合, 否則會產生編譯錯誤.
Input<Func> input{ "input", 3 };  // require 3-dimensional Func,
                                  // but leave Type unspecified  
當一個使用了 Input<ltFunc> 的 Generator 被直接編譯時 (e.g. 使用 GenGen), 而型別與(或)維度在程式碼中並沒有被指明, 此時 Input<ltFunc> 必須明確地指明. 能夠透過 GeneratorParams 搭配Input 或 Output 衍生出的名字來指定. (在上述的例子中, input 有著名為 "input.type" 的 GeneratorParam 與名為 "input.dim" 的 GeneratorParam)

明確宣告輸出


所有 Generator 的輸入可以藉由其成員來反推, 但是關於輸出的資訊在先前僅能藉由呼叫 build() 並確認回傳值來確定(可能是個 Func  或是一個 Pipeline)
而這次的改變, Generator 相對地能明確地宣告其輸出為一個成員變數, 並且提供 generator() 介面函式取代 build() (等同於 generate 沒有回傳任何數值)

範例:
class SumColumns : Generator<SumColumns> {
  Input<Func> input{"input", Float(32), 2};
  Input<int32_t> width{"width"};

  Func build() {
    RDom r(0, width);
    Func f;
    Var y;
    f(y) = 0.f;
    f(y) += input(r.x, y);
    return f;
  }
};
轉變為:
class SumColumns : Generator<SumColumns> {
  Input<Func> input{"input", Float(32), 2};
  Input<int32_t> width{"width"};

  Output<Func> sum_cols{"sum_cols", Float(32), 1};

  void generate() {
    RDom r(0, width);
    Var y;
    sum_cols(y) = 0.f;
    sum_cols(y) += input(r, y);
  }
};

當搭配  Output<Func>, 能夠選擇性地不指定 Output<Func> 的型別與/或維度; 當任何未被指明的型別必須透過 GeneratorParam 指定才能夠使用高階編譯.
請注意 Output<> 僅能搭配 Generator 使用, 並不能在其他 Halide 程式使用.
Generator 的底層實作將會驗證(在呼叫 generate() 之後)所有定義的輸出, 與定義是否符合宣告.


開發者能夠透過 Type 清單來指明一個回傳為 Tuple 的輸出:
class Tupler : Generator<Tupler> {
  Input<Func> input{"input", Int(32), 2};
  Output<Func> output{"output", {Float(32), UInt(8)}, 2};

  void generate() {
    Var x, y;
    output(x, y) = Tuple(cast(input(x, y)), cast<uint8_t>(input(x, y)));
  }
};
一個 Generator 能夠定義多個 Output (相當於悄悄地實作為一個 Pipeline):
class SumRowsAndColumns : Generator<SumRowsAndColumns> {
  Input<Func> input{"input", Float(32), 2};
  Input<int32_t> width{"width"};
  Input<int32_t> height{"height"};

  Output<Func> sum_rows{"sum_rows", Float(32), 1};
  Output<Func> sum_cols{"sum_cols", Float(32), 1};

  void generate() {
    RDom rc(0, height);
    Var x;
    sum_rows(x) = 0.f;
    sum_rows(x) += input(x, rc);

    RDom rr(0, width);
    Var y;
    sum_cols(y) = 0.f;
    sum_cols(y) += input(rr, y);
  }
};
並且允許指定輸出為任意的純量型別(除了 Handle 型別外); 這僅只是在零維度 Func 上的語法糖果, 但相當的方便, 特別是用在多個輸出上:
class Sum : Generator<Sum> {
  Input<Func> input{"input", Float(32), 2};
  Input<int32_t> width{"width"};
  Input<int32_t> height{"height"};

  Output<Func> sum_rows{"sum_rows", Float(32), 1};
  Output<Func> sum_cols{"sum_cols", Float(32), 1};
  Output<float> sum{"sum"};

  void generate() {
    RDom rc(0, height);
    Var x;
    sum_rows(x) = 0.f;
    sum_rows(x) += input(x, rc);

    RDom rr(0, width);
    Var y;
    sum_cols(y) = 0.f;
    sum_cols(y) += input(rr, y);

    RDom r(0, width, 0, height);
    sum() = 0.f;
    sum() += input(r.x, r.y);
  }
}; 

請注意同時定義 build() 與 generate() 對 Halide 是一個錯誤. 

輸入與輸出陣列


藉由使用一個型別陣列來作為型別參數, 亦能夠透過使用新的語法來宣告一個輸入或輸出的陣列:
// Takes exactly 3 images and outputs exactly 3 sums.
class SumRowsAndColumns : Generator<SumRowsAndColumns> {
  Input<Func[3]> inputs{"inputs", Float(32), 2};
  Input<int32_t[2]> extents{"extents"};
  Output<Func[3]> sums{"sums", Float(32), 1};
  void generate() {
    assert(inputs.size() == sums.size());
    // assume all inputs are same extent
    Expr width = extent[0];
    Expr height = extent[1];
    for (size_t i = 0; i < inputs.size(); ++i) {
      RDom r(0, width, 0, height);
      sums[i]() = 0.f;
      sums[i]() += inputs[i](r.x, r.y);
     }
  }
}; 
陣列大小能夠不指定, 但有著注意事項:
  • 對於 ahead-of-time 編譯, 輸入必須有著明確的大小, 編譯時藉由一個  GeneratorParam 帶入 (e.g., pyramid.size=3)
  • 對於透過一個 Stub 的 JIT 編譯, 輸入陣列大小將自傳入的向量推斷.
  • 對於 ahead-of-time 編譯, 輸出可能在編譯時藉由一個 GeneratorParam 指定了具體的大小 (e.g., pyramid.size=3), 或是能夠透過 resize 的方式指定大小.
class Pyramid : public Generator<Pyramid> {
public:
    GeneratorParam<int32_t> levels{"levels", 10};
    Input<Func> input{ "input", Float(32), 2 };
    Output<Func[]> pyramid{ "pyramid", Float(32), 2 };
    void generate() {
        pyramid.resize(levels);
        pyramid[0](x, y) = input(x, y);
        for (int i = 1; i < pyramid.size(); i++) {
            pyramid[i](x, y) = (pyramid[i-1](2*x, 2*y) +
                               pyramid[i-1](2*x+1, 2*y) +
                               pyramid[i-1](2*x, 2*y+1) +
                               pyramid[i-1](2*x+1, 2*y+1))/4;
        }
    }
};
一個未指定大小的 Input/Output 陣列必須在高階編譯時有著明確的大小; 現在有著基於命名的 GeneratorParam 用來設定此數值, (在上面的範例為"pyramid.size").
請注意 Input 與 Output 陣列都支援來自 std::vector<> 有限的介面:
  • operator[]
  • size()
  • begin()
  • end()
  • resize() (Output only)
    
    

自 Build() 中分開排程

一個 Generator 現在能夠將既存的 build() 分為兩個介面函式:
void generate() { ... }
void schedule() { ... }

如此一個 Generator 必須將所有對於中介 Func 排程的程式碼移到 schedule() . 請注意這表示可排程的 Func, Var, 等等將需要宣告為 Generator 的成員變數. (因為
Output<> 需要宣告為成員變數, 這相當簡單, 但對於中介需要排程 Func 可能需要相當變動.)

範例:
class Example : Generator<Example> {
  Output<Func> output{"output", Float(32), 2};

  void generate() {
    Var x, y;

    Func intermediate;
    intermediate(x, y) = SomeExpr(x, y);

    output(x, y) = intermediate(x, y);

    intermediate.compute_at(output, y);
  }
};
轉變為:
class Example : Generator<Example> {
  Output<Func> output{"output", Float(32), 2};

  void generate() {
    intermediate(x, y) = SomeExpr(x, y);
    output(x, y) = intermediate(x, y);
  }

  void schedule() {
    intermediate.compute_at(output, y);
  }

  Func intermediate;
  Var x, y;
};
請注意輸出的 Func 並不需要一個對於 compute_at()store_at() 的排程指令: 要不就是潛在地 compute_root() (當直接編譯為濾波器), 或是透過呼叫者來排程 (如後面將看到的, 當被用作為子項目元件).

即便中介的 Halide 程式碼並沒有必要的排程 (e.g. 全為 inline), 依然應該提供一個空的schedule() 來使其明確與簡潔.

範例:
class ExampleInline : Generator<ExampleInline> {
  Output<Func> output{"output", Float(32), 2};

  void generate() {
    Var x, y;
    output(x, y) = SomeExpr(x, y);
  }
};
轉變為
class ExampleInline : Generator<ExampleInline> {
  Output<Func> output{"output", Float(32), 2};

  void generate() {
    output(x, y) = SomeExpr(x, y);
  }

  void schedule() {
    // empty
  }

  Var x, y;
};

當需要時將 GeneratorParam 轉為 ScheduleParam


GeneratorParam 現在透過新的 ScheduleParam 型別來改進. 所有被 schedule() 所使用的 generator 參數必須宣告為 ScheduleParam 而非 GeneratorParam. 這有著兩個目的:
  • 允許在任意的 Generator 之間陳列與傳遞排程資訊的方式 (之後將會看到).
  • 使得哪些 GeneratorParams 是作為排程用途, 這有助於未來的 Autoscheduler 運作.
請注意有些一般性的 GeneratorParam 約定早已如同 ScheduleParam 動作(特別是 vectorizeparallelize); 這僅只是型式化先前的慣例.

GeneratorParamScheduleParam 將共存在單一個 namespace (i.e., 將 GeneratorParamScheduleParam 以相同命名宣告是個錯誤).

當一個 GeneratorParam 能用在 Generator 中任何地方使用 (無論 generate() 或者schedule() ), 一個 ScheduleParam 僅能夠限縮在 schedule() 中存取. (在未來, 這將可能會是編譯時期的限制)

請注意當 GeneratorParam 必須能夠對字串雙向序列化 (如同 GeneratorParams), 一些 ScheduleParam 數值並非可序列化, 因為可能參照了執行時期的 Halide 結構 (特別是, LoopLevel, 並不能在一般情況下可靠地以命名指定). 嘗試從 GenGen 設定如此的 ScheduleParam  將產生編譯時期錯誤.

範例:
class Example : Generator<Example> {
  GeneratorParam<int32_t> iters{"iters", 10};
  GeneratorParam<bool> vectorize{"vectorize", true};

  Func generate() {
    Var x, y;
    vector<Func> intermediates;
    for (int i = 0; i < iters; ++i) {
      Func g;
      g(x, y) = (i == 0) ? SomeExpr(x, y) : SomeExpr2(g(x, y));
      intermediates.push_back(g);
    }
    Func f;
    f(x, y) = intermediates.back()(x, y);

    // Schedule
    for (auto fi : intermediates) {
      fi.compute_at(f, y);
      if (vectorize) fi.vectorize(x, natural_vector_size());
    }
    return f;
  }
};
轉變為:
class Example : Generator<Example> {
  GeneratorParam<int32_t> iters{"iters", 10};
  ScheduleParam<bool> vectorize{"vectorize", true};

  Output<Func> output{"output", Float(32), 2};

  void generate() {
    for (int i = 0; i < iters; ++i) {
      Func g;
      g(x, y) = (i == 0) ? SomeExpr(x, y) : SomeExpr2(g(x, y));
      intermediates.push_back(g);
    }
    output(x, y) = intermediates.back()(x, y);
  }

  void schedule() {
    for (auto fi : intermediates) {
      fi.compute_at(output, y);
      if (vectorize) fi.vectorize(x, natural_vector_size());
    }
  }

  Var x, y;
  vector<Func> intermediates;
};
請注意 ScheduleParam 也能夠有著其他有趣的數值, 特別是 LoopLevel:
class Example : Generator<Example> {
  // Specify a LoopLevel at which we want intermediate Func(s)
  // to be computed and/or stored.
  ScheduleParam<LoopLevel> intermediate_compute_level{"level", "undefined"};
  ScheduleParam<LoopLevel> intermediate_store_level{"level", "root"};
  Output<Func> output{"output", Float(32), 2};

  void generate() {
    intermediate(x, y) = SomeExpr(x, y);
    output(x, y) = intermediate(x, y);
  }

  void schedule() {
    intermediate
      // If intermediate_compute_level is undefined,
      // default to computing at output's rows
      .compute_at(intermediate_compute_level.defined() ?
                  intermediate_compute_level :
                  LoopLevel(output, y))
      .store_at(intermediate_store_level);
  }

  Func intermediate;
  Var x, y;
};
請注意 ScheduleParam<LoopLevel> 能夠預設為 "root", "inline", 或 "undefined"; 所有其他的數值 (e.g. Func-and-Var) 必須在實際程式碼中指定. (LoopLevel(Func, Var) 明確地無法藉由像是 "func.var" 的命名來指定; 雖然 Halide 內部使用這慣例, 但目前無法保證在任意集合的 Generator 使用唯一的 Func 命名)

請注意使用一個未定義的 LoopLevel 作為排程是一個錯誤

修訂的 RegisterGenerator 語法


在以往, 必須註冊藉由在廣域下明確的實體化一個 RegisterGenerator 來註冊 Generator:
Halide::RegisterGenerator<MyGen> register_my_gen{"my_gen"};
這依然有效, 但導入了更簡單的註冊巨集:
HALIDE_REGISTER_GENERATOR(MyGen, my_gen) // no semicolon at end
若想要針對特定 Generator 產生一個 Stub, 必須使用這新的註冊巨集, 並加入該資訊到宣告之中:
// We must forward-declare the name we want for the stub, 
// inside the proper namespace(s). None of the namespace(s) 
// may be anonymous (if they are, failures will occur at Halide
// compilation time).
namespace SomeNamespace { class MyGenStub; }
HALIDE_REGISTER_GENERATOR(MyGen, "my_gen", SomeNamespace::MyGenStub)
若對於第3個參數的 stub 合格命名沒有被適當地宣告, 將產生編譯錯誤. 完全合格的命名必須至少有著一個 namespace (也就是, 一個全域命名是不會被接受的).

Generator Stubs


先以使用範例做開始, 接著再回頭解釋過程. 若有著一個想重複使用的 RGB-to-YCbCr 元件:
class RgbToYCbCr : public Generator<RgbToYCbCr> {
  Input<Func> input{"input", Float(32), 3};
  Output<Func> output{"output", Float(32), 3};
  void generate() { ... conversion code here ... }
  void schedule() { ... scheduling code here ... }
};
RegisterGenerator register_me{"rgb_to_ycbcr"};
GenGen 現在能夠產生一個類Func 的包覆一個 generator 的 stub 類別, (通常) 以一個副檔名 為 ".stub.h" 檔案發出. 看起來像是:
/path/to/rgb_to_rcbcr.stub.h:

  // MACHINE-GENERATED
  class RgbToYCbCr : public GeneratorStub {
    struct Inputs { 
       // All the Input<>s declared in the Generator are listed here,
       // as either Func or Expr
       Func input;
    };
    struct GeneratorParams { ... };
    struct ScheduleParams { ... };

    // ctor, with required inputs, and (optional) GeneratorParams.
    RgbToYCbCr(GeneratorContext* context,
               const Inputs& inputs,
               const GeneratorParams& = {}) { ... }

    // Output(s)
    Func output;

    // Overloads for first output
    operator Func() const { return output; }
    Expr operator()(Expr x, Expr y, Expr z) const  { return output(x, y, z); }
    Expr operator()(std::vector<Expr> args) const  { return output(args); }
    Expr operator()(std::vector<Var> args) const  { return output(args); }

    void schedule(const ScheduleParams &params = {});
  };
請注意這只是一個 "header-only" 類別; 所有的函式都是 inline 方式(或是 template-multilinked, etc) 所以並沒有對應要合入的 .cpp . 並且注意到有一個 "by-value", 基於內部處理的類別, 如同多數其他在 Halide 中的類別(e.g. Func, Expr, etc).

接著像是如此般地使用:
#include "/path/to/rgb_to_rcbcr.stub.h"

class AwesomeFilter : public Generator<AwesomeFilter> {
 public:
  Input<Func> input{"input", Float(32), 3};
  Output<Func> output{"output", Float(32), 3};

  void generate() {
    // Snap image into buckets while still in RGB.
    quantized(x, y, c) = Quantize(input(x, y, c));

    // Convert to YCbCr.
    rgb_to_ycbcr = RgbToYCbCr(this, {quantized});

    // Do something awesome with it. Note that rgb_to_ycbcr autoconverts to a Func.
    output(x, y, c) = SomethingAwesome(rgb_to_ycbcr(x, y, c));
  }
  void schedule() {
    // explicitly schedule the intermediate Funcs we used
    // (including any reusable Generators).
    quantized.
      .vectorize(x, natural_vector_size<float>())
      .compute_at(rgb_to_ycbcr, y);
    rgb_to_ycbcr
      .vectorize(x, natural_vector_size<float>())
      .compute_at(output, y);

    // *Also* call the schedule method for all reusable Generators we used,
    // so that they can schedule their own intermediate results as needed.
    // (Note that we may have to pass them appropriate values for ScheduleParam,
    // which vary from Generator to Generator; since RgbToYCbCr has none,
    // we don't need to pass any.)
    rgb_to_ycbcr.schedule();
 }

 private:
  Var x, y, c;
  Func quantized;
  RgbToYCbCr rgb_to_ycbcr;

  Expr Quantize(Expr e) { ... }
  Expr SomethingAwesome(Expr e) { ... }
};
值得一提的是所有對子元件的輸入必須在子元件建構時明確地提供 (以其 ctor 參數的方式); 呼叫者必須負責提供這些資訊 (並沒有自呼叫者到子元件自動傳遞的概念)

RgbToYCbCr 有著輸入或輸出陣列會如何? 像是:
class RgbToYCbCrMulti : public Generator<RgbToYCbCrMulti> {
  Input<Func[3]> inputs{"inputs", Float(32), 3};
  Input<float> coefficients{"coefficients", 1.f};
  Output<Func[3]> outputs{"outputs", Float(32), 3};
  ...
};
在這例子中, 所產生的 RgbToYCbCrMulti 類別對於輸入需要 vector-of-Func (或 vector-of-Expr) , 並且提供 vector-of-Func 作為輸出成員:
class RgbToYCbCrMulti : public GeneratorStub {
    struct Inputs { 
       std::vector<Func> inputs;
       std::vector<Expr> coefficients;
    };
    RgbToYCbCr(GeneratorContext* context,
               const Inputs& inputs,
               const GeneratorParams& = {}}) { ... }

    ...

    std::vector<Func> outputs;
};
RgbToYCbCr 有著多個輸出會如何? 像是:
class RgbToYCbCrMulti : public Generator<RgbToYCbCrMulti> {
  Input<Func> input{"input", Float(32), 3};
  Output<Func> output{"output", Float(32), 3};
  Output<Func> mask{"mask", UInt(8), 2};
  Output<float> score{"score"};
  ...
};
在這例子中, 所產生的  RgbToYCbCrMulti 類別將所有的輸出搭配對應宣告在 Generator 中的命名作為結構成員:
struct RgbToYCbCrMulti {
    ...
    Func output;
    Func mask;
    Func score;
};
請注意純量輸出為了一致性依然以(零維度)函式方式表現. (並且注意到 "output" 並不是一個特別的命名; 只是發生在 Generator 第一個輸出的命名上)

亦請注意第一個輸出總是以 "是" 與 "有著"  的關係表示: RgbToYCbCrMulti 複載了必要的運算元, 因此使用其 "output" 就如同使用以一個 Func, 在此例中:
struct RgbToYCbCrMulti {
    ...
    Func output;

    operator Func() const { return output; }
    Expr operator()(Expr x, Expr y, Expr z) const  { return output(x, y, z); }
    Expr operator()(std::vector<Expr> args) const  { return output(args); }
    Expr operator()(std::vector<Var> args) const  { return output(args); }
    ...
};
這(的確)是多餘的, 但這是考量過的: 這讓大部分(單一輸出)的情況更為便利,  然而對多輸出的情況沒有影響.

consumer 可能如此般地使用:
#include "/path/to/rgb_to_rcbcr_multi.stub.h"

class AwesomeFilter : public Generator<AwesomeFilter> {
  ...
  void generate() {
    rgb_to_ycbcr_multi = RgbToYCbCrMulti(this, {input});
    output(x, y, c) = SomethingAwesome(rgb_to_ycbcr_multi.output(x, y, c),
                                       rgb_to_ycbcr_multi.mask(x, y),
                                       rgb_to_ycbcr_multi.score());
  }
  void schedule() {
    rgb_to_ycbcr_multi.output
      .vectorize(x, natural_vector_size<float>())
      .compute_at(output, y);
    rgb_to_ycbcr_multi.mask
      .vectorize(x, natural_vector_size<float>())
      .compute_at(output, y);
    rgb_to_ycbcr_multi.score
      .compute_root();
    // Don't forget to call the schedule() function.
    rgb_to_ycbcr_multi.schedule();
  }
};
RgbToYCbCr 之中有著想要設定來做 code generation 的 GeneratorParam 會如何? 在這情況, 當呼叫建構子時將傳入一個數值到一個 generator_params 的選擇性欄位.
class RgbToYCbCr : public Generator<RgbToYCbCr> {
  GeneratorParam<Type> input_type{"input_type", UInt(8)};
  GeneratorParam<bool> fast_but_less_accurate{"fast_but_less_accurate", false};
  ...
};
這可能會產生一個不一樣定義的 GeneratorParams, 其有著每個 GeneratorParam 欄位並以適當預設值作初始:
struct GeneratorParams {
  Halide::Type input_type{UInt(8)};
  bool fast_but_less_accurate{false};
};
接著可能手動地填入:
class AwesomeFilter : public Generator<AwesomeFilter> {
  void generate() {
    ...
    GeneratorParams generator_params;
    generator_params.input_type = Float(32);
    generator_params.fast_but_less_accurate = true;
    rgb_to_ycbcr = RgbToYCbCr(this, input, generator_params);
    ...
  }
}
或者, 在 C++ 編譯時期得知其型別時, 能更精練地使用一個範本化建構子(templated constructor):
class AwesomeFilter : public Generator<AwesomeFilter> {
  void generate() {
    ...
    rgb_to_ycbcr = RgbToYCbCr::make<float, true>(this, input);
    ...
  }
}
RgbToYCbCr 中有著 ScheduleParam 會如何?
class RgbToYCbCr : public Generator<RgbToYCbCr> {
  ScheduleParam<LoopLevel> level{"level"};
  ScheduleParam<bool> vectorize{"vectorize"};

  void generate() {
    intermediate(x, y) = SomeExpr(x, y);
    output(x, y) = intermediate(x, y);
  }

  void schedule() {
    intermediate.compute_at(level);
    if (vectorize) intermediate.vectorize(x, natural_vector_width<float>());
  }

  Var x, y;
  Func intermediate;
};
在這情形中, 產生的 stub 程式碼對於 ScheduleParams 將有著一個不一樣的宣告:
struct ScheduleParams {
  LoopLevel level{"undefined"};
  bool vectorize{false};
};

並且可能如此地呼叫:
class AwesomeFilter : public Generator<AwesomeFilter> {
  ...
  void schedule() {
    rgb_to_ycbcr
      .vectorize(x, natural_vector_size<float>())
      .compute_at(output, y);

    rgb_to_ycbcr.schedule({
      // We want any intermediate products also at compute_at(output, y)
      LoopLevel(output, y),
      // vectorization: yes please
      true
    });
  }
  ...
}


2018年3月15日 星期四

Halide - 實務心得 1

近日著手開始以 Halide 撰寫 AOT 程式 發現有一些特別的東西沒有寫在教學文件上

 

 1. buffer pointer of Halide::Runtime::Buffer

宣告一個 Buffer object 並不困難像是:
Halide::Runtime::Buffer buf(640, 480);
而在教學文件中很常看到 Buffer 的使用 
但是範例上常常是直接載入 jpg/png 圖檔, 
或是透過下列的方式來讀取或寫入 (x, y) 位置的數值: (見 Lesson 10 part 2)
buf(x, y) = val

而對於 C/C++ 而言是可以取得 buffer pointer 的
取得方式也很簡單, 透過 Buffer::data() 這個成員函式的呼叫
以上述例子就是:
unsigned char* buf_ptr = (unsigned char*)(buf.data());
使用由 data() 取得的 pointer 就必須注意 layout, 預設的是 Planar
若需要使用 interleaving data 則要透過 set_stride  來設定 (見 Lesson 16)

2. Buffer allocation

儘管一般開發 Halide 程式使用的裝置離不開 CPU/GPU
但是 Halide 俱備介面, 只要實作 code generation 與 runtime interface 就可以接上
以 Qualcomm HVX 來說, DSP backend 的 Buffer 使用與 CPU 有相當不同
透過的是 Buffer::device_malloc() 來配置特定計算裝置的記憶體
並以此避免 memcpy 的呼叫, 來達到 zero-copy 的目的
而 device_malloc 需要帶入 Device API 的 interface
像是在 Qualcomm 官方 Halide 文件中提到 Qualcomm Hexagon DSP buffer 配置的呼叫為:
buf.device_malloc(halide_hexagon_device_interface());
如此即可配置針對 Hexagon DSP 所使用的記憶體空間, 並且無需 data copy

3. Prefetch schedule directive

計算單元與記憶體間資料的傳輸方式一直是改善效能的重點
個人經驗是撰寫 SIMD code 透過 prefetch 能改善不少效能
其他像是 GPU 俱備 local/shared memory,
有些 DSP 也有著自己的 scratchpad memory
透過重疊(overlapping) 執行時間與傳輸時間
這部分都牽涉到 2D 的 Prefetch 或 DMA 的呼叫使用

在 2015 年剛接觸 Halide 時, 與當時熟悉 Halide 同事討論
當時給的答覆是並無法對應 prefetch 這樣的行為
而日前第一時間在 Qualcomm Halide SDK 中的範例看到 prefetch 排程指令時:
(像是 Halide SDK 中  HALIDE_Tools/2.0/Halide/hexagon_benchmarks/gaussian5x5_generator.cpp )
output.hexagon().prefetch(input, y, 2) ... ;
個人一度以為這只是 Qualcomm 為了 Hexagon DSP 而客製的 API
然而近日發現 Halide 中的 Prefetch.h說明文件
又挖了 History 後才發現在 2016 年中 Halide 已加入了對於 prefetch 行為的支援
這也讓 Halide 能夠更彈性的支援不同處理器的運作特性
而類似於 data 在 Boundary 的處理 Prefetch 也有 PrefetchBoundStrategy

4. AOT target

在 Lesson 15 並沒有提及 Generator Enhancements
在後續 Halide 對於 Generator 分為 generate() 與 schedule()
顧名思義, generate() 是用來定義 pipeline 的, 而 schedule() 是處理排程的
若希望有效地支援不同平台, 就需要知道 code generation target 的特性
這裡主要透過二個 API 與二個 public member
* Halide::GeneratorContext::get_target() - 這可以得到 Halide::Target 的物件, 透過這物件可以獲得相關資訊, 最直接的就是 Target::to_string() 可以得到指令傳入的 target 字串
* Halide::Target::has_feature() - 可以確認是否支援某些特性像是 SSE/AVX, ARMv7s or HVX
* Halide::Target.arch - X86, ARM, Hexagon ...
* Halide::Target.os - Linux, Windows ...
如此可以透過條件判別來針對平台特性呼叫在 schedule() 中呼叫對應的排程指令


簡短列出這些日子看到的有趣東西, 一方面也做個紀錄


2018年3月5日 星期一

Halide Tutorial 非官方中文翻譯 - Part 8

這是最後四節的中譯部分
官方 Tutorial 的翻譯也到此告一段落
若有不恰當, 錯譯或是更好的表示方式還請不吝指教
以下為前7個部分的聯結

Halide Tutorial 非官方中譯 - Part 1
Halide Tutorial 非官方中譯 - Part 2
Halide Tutorial 非官方中譯 - Part 3
Halide Tutorial 非官方中譯 - Part 4
Halide Tutorial 非官方中譯 - Part 5
Halide Tutorial 非官方中譯 - Part 6
Halide Tutorial 非官方中譯 - Part 7

後續希望能針對應用這些技巧提供心得與 case study
==

Lesson 19: 包裝函式 (Wrapper Funcs)

// 這節示範如何使用 Func::in 與 ImageParam::in 對一個函式在不同地方以不同方式排程
// 並且階段地自 Func 或 ImageParam 載入

#include "Halide.h"
#include  <stdio.h>

using namespace Halide;

int main(int argc, char **argv) {
    // 宣告後續會使用到的變數
    Var x("x"), y("y"), xo("xo"), yo("yo"), xi("xi"), yi("yi"); 
    // 這節將使用 Func::in 與 ImageParam::in 指令來包裝一個 Func 或是 ImageParam
    {
        // 考量一個簡單二階 pipeline
        Func f("f_local"), g("g_local");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y) + 3;

        f.compute_root();

        // 這將產生下列巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     g(x, y) = 2 * f(x, y) + 3

        // 使用 Func::in 能夠在介於 f 與 g 之間插入一個單獨排程的函式
        Func f_in_g = f.in(g);
        f_in_g.compute_root(); 
        // 同樣地, 能夠如下地串聯起來:
        // f.in(g).compute_root();
        // 這將產生下列巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     f_in_g(x, y) = f(x, y)
        // for y:
        //   for x:
        //     g(x, y) = 2 * f_in_g(x, y) + 3

        g.realize(5, 5); 
        // 下圖為視覺化結果
 
        // 排程指令 f.in(g) 會以包裝後的函式取代在 'g' 中 'f' 的呼叫.
        // 實際上它以下列方式改寫了上面原本的 pipeline
        {
            Func f_in_g("f_in_g"), f("f"), g("g");
            f(x, y) = x + y;
            f_in_g(x, y) = f(x, y);
            g(x, y) = 2 * f_in_g(x, y) + 3;

            f.compute_root();
            f_in_g.compute_root();
            g.compute_root();
        }
        // 在隔離方式上, 像轉換的計算似乎沒多大意義, 但這能被應用在多樣的排程技巧上
    }

    { 
        // 在上面的排程, 僅有 'g' 所產生的 'f' 呼叫被取代
        // 其他對 f 產生的呼叫依然直接地使用 'f'
        // 若希望廣域地以一個包裝函式取代所有 'f' 的呼叫, 那麼就使用 f.in()
        // 考慮下列有兩個 stage 使用 f 的三階 pipeline:
        Func f("f_global"), g("g_global"), h("h_global");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y);
        h(x, y) = 3 + g(x, y) - f(x, y);
        f.compute_root();
        g.compute_root();
        h.compute_root(); 
        // 將所有在 'g' 與 'h' 中的 'f' 呼叫以一個包裝函式取代:
        f.in().compute_root();

        // 等義的巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     f_in(x, y) = f(x, y)
        // for y:
        //   for x:
        //     g(x, y) = 2 * f_in(x, y)
        // for y:
        //   for x:
        //     h(x, y) = 3 + g(x, y) - f_in(x, y)

        h.realize(5, 5); 
        // 下圖為視覺化結果
 
    }

    { 
        // 亦能夠給予 g 與 h 屬於各自的 f 包裝函式
        // 這次在各自巢狀 loop 內部排程, 這是無法以單一包裝函式所達到的

        Func f("f_unique"), g("g_unique"), h("h_unique");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y);
        h(x, y) = 3 + g(x, y) - f(x, y);

        f.compute_root();
        g.compute_root();
        h.compute_root();

        f.in(g).compute_at(g, y);
        f.in(h).compute_at(h, y);

        // 這產生了下列巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     f_in_g(x, y) = f(x, y)
        //   for x:
        //     g(x, y) = 2 * f_in_g(x, y)
        // for y:
        //   for x:
        //     f_in_h(x, y) = f(x, y)
        //   for x:
        //     h(x, y) = 3 + g(x, y) - f_in_h(x, y)

        h.realize(5, 5);
        // 下圖為視覺化結果

     }

    {
        // 至今這可能像是許多無意義的記憶體複製
        // 而 Func::in 能與其他排程指令合併來達到不同的用途
        // 首先確認對於不同的函式建立了個別的實現並且以不同方式排程
 
        // 以幾乎相同的 pipeline 開始
        Func f("f_sched"), g("g_sched"), h("h_sched");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y);
        // h 使用 f 的一個較遠區域
        h(x, y) = 3 + g(x, y) - f(x + 93, y - 87);

        // 這次以 inline 方式使用 f
        // f.compute_root();
        g.compute_root();
        h.compute_root();

        f.in(g).compute_at(g, y);
        f.in(h).compute_at(h, y);

        // g 與 h 現在透過個別的包裝函式來呼叫 f. The wrappers are
        // 排程是透過包裝函式, 若無則 f 是以 inline 方式置於兩個包裝函式中.
        // 這兩者將獨立計算呼叫其所需的 f 區域.
        // 當以 compute_root 排程 f 函式, 將會計算 g 與 h 所需範圍的區域方塊,
        // 而當中可能多數都是不被需要的資料.
 
        // 因此亦能夠各自不同地排程
        // 為了排程, 包裝函式自原始函式繼承了純變數, 
        // 因此能夠使用在定義 f 時使用一樣的 x 與 y:
        f.in(g).vectorize(x, 4);
        f.in(h).split(x, xo, xi, 2).reorder(xo, xi);

        // 請注意乎要 f.in(g) 在第二次呼叫時回傳的是在第一次呼叫已被建構的包裝函式,
        // 並不會建構一個新的.

        h.realize(8, 8);
        // 下圖為視覺化結果


        // 請注意由於 f 是以 inline 方式被兩個包裝函式使用,
        // 因此是包裝函式執行 f 的計算工作, 而非只是從已經計算好的資料載入
    }

    {
        // Func::in 對於透過 Func 階段載入可能在堆疊或GPU 共享記憶體的中介 buffer 相當有用
 
        // 考量一個轉置部份 compute_root 計算後結果的 pipeline:

        Func f("f_transpose"), g("g_transpose");
        f(x, y) = sin(((x + y) * sqrt(y)) / 10);
        f.compute_root();

        g(x, y) = f(y, x);

        // 執行策略上要載入一個 f 的 4x4 tile 到暫存器中, 於暫存器中轉置並寫到 g 的 4x4 tile.
        // 這裡使用 Func::in 來表示:
 
        Func f_tile = f.in(g);

        // 如此有了一個三階 pipeline:
        // f -> f_tile -> g 
        // f_tile 會載入多個 f 向量, 寫入暫存器中並轉置, g 將這些資料寫回主記憶體.
        g.tile(x, y, xo, yo, xi, yi, 4, 4)
            .vectorize(xi)
            .unroll(yi);
 
        // 這裡將在 g 的 tile 計算 f_transpose,
        // 並且使用 Func::reorder_storage 來說明 f_transpose 應以 column-major 方式存放
        // 因此對於 g 的載入這些資料能以向量載入的形式達成
       f_tile.compute_at(g, xo)
            .reorder_storage(y, x)
            .vectorize(x)
            .unroll(y); 
        // 必須要小心確認 f_transpose 只以常數索引方式存取
        // 全部於各自 compute_at 層次存在的 loop 的 unrolling/vectorization
        // 有著這樣的特性. 僅被常數索引存取的記憶體配置將會被提升到暫存器中.
 
        g.realize(16, 16);
        // 下圖為視覺化結果

    }

    {
        // ImageParam::in 如同 Func::in 的行為, 能夠以類似的方式做階段載入
        // 不再以轉置為例, 而用 ImageParam::in 來階段載入輸入影像至 GPU 共享記憶體,
        // 如同明確管理的 cache 一般有效地使用 shared/local memory.
        ImageParam img(Int(32), 2);

        // 計算些微模糊化的輸入影像
        Func blur("blur");
        blur(x, y) = (img(x - 1, y - 1) + img(x, y - 1) + img(x + 1, y - 1) +
                      img(x - 1, y    ) + img(x, y    ) + img(x + 1, y    ) +
                      img(x - 1, y + 1) + img(x, y + 1) + img(x + 1, y + 1));

        blur.compute_root().gpu_tile(x, y, xo, yo, xi, yi, 8, 8); 
        // 以 ImageParam::in 建構的包裝函式 有著以 _0, _1 .. 等等所命名的純變數
        // 對 blur 每個 tile 的變數排程, 並且將 _0 與 _1 對應到 gpu 執行緒
        img.in(blur).compute_at(blur, xo).gpu_threads(_0, _1);
         // Without Func::in, computing an 8x8 tile of blur would do
        // 8*8*9 loads to global memory. With Func::in, the wrapper
        // does 10*10 loads to global memory up front, and then blur
        // does 8*8*9 loads to shared/local memory.
 
        // 在不使用 Func::in 情況下, blur 計算 8x8 tile 需要 8*8*9 次 global memory 載入
        // 而使用 With Func::in, 包裝函式於開始會作 10*10 次 global memory 載入
        // 接著 blur 函式會自 shared/local memory 作 8*8*9 次載入

        // 如同 Lesson 12 , 選擇適當的 GPU API
        Target target = get_host_target();
        if (target.os == Target::OSX) {
            target.set_feature(Target::Metal);
        } else {
            target.set_feature(Target::OpenCL);
        }

        // 建立一個有趣的輸入影響
        Buffer input(258, 258);
        input.set_min(-1, -1);
        for (int y = input.top(); y <= input.bottom(); y++) {
            for (int x = input.left(); x <= input.right(); x++) {
                input(x, y) = x * 17 + y % 4;
            }
        }

        img.set(input);
        blur.compile_jit(target);
        Buffer<int> out = blur.realize(256, 256);

        // 確認結果一致
        for (int y = out.top(); y <= out.bottom(); y++) {
            for (int x = out.left(); x <= out.right(); x++) {
                int val = out(x, y);
                int expected = (input(x - 1, y - 1) + input(x, y - 1) + input(x + 1, y - 1) +
                                input(x - 1, y    ) + input(x, y    ) + input(x + 1, y    ) +
                                input(x - 1, y + 1) + input(x, y + 1) + input(x + 1, y + 1));
                if (val != expected) {
                    printf("out(%d, %d) = %d instead of %d\n",
                           x, y, val, expected);
                    return -1;
                }
            }
        }
    }

    {
        // Func::in 也能夠用來以相同巢狀 loop 將多階合併為單一函式
        // 考量下列計算每個 pixel 數值的並由左至右翻轉每個 scanline 的 pipeline
        Func f("f_group"), g("g_group"), h("h_group");

        // 初始化 f
        f(x, y) = sin(x - y);
        RDom r(1, 7);

        // 由左至右翻轉
        f(r, y) = (f(r, y) + f(r - 1, y)) / 2;

        // 由右至左翻轉
        f(7 - r, y) = (f(7 - r, y) + f(8 - r, y)) / 2;
        // 接著嘗試複雜的存取模式: 一個循環的 45 度旋轉
        g(x, y) = f((x + y) % 8, (x - y) % 8);

        // 由於以複雜的方式存取, 因此 f 應以 compute_root 排程
        // 但這也表示 f 所有的階層將在巢狀 loop 中計算:

        // for y:
        //   for x:
        //     f(x, y) = sin(x - y)
        // for y:
        //   for r:
        //     f(r, y) = (f(r, y) + f(r - 1, y)) / 2
        // for y:
        //   for r:
        //     f(7 - r, y) = (f(7 - r, y) + f(8 - r, y)) / 2
        // for y:
        //   for x:
        //     g(x, y) = f((x + y) % 8, (x - y) % 8);
        // 若以共同的 loop y 來對 f 作排程能夠得到更佳的區域特性(locality)
        // 能夠以包裝函式計算 f 的 scanline  來達成:
        f.in(g).compute_root();
        f.compute_at(f.in(g), y);
        // 對於一個函式有更新階 (update stage)時, 
        // f 有著預設在 consumer 最內層 loop 計算的排程, 這裡也就是 f.in(g)
        // 因此產生如下有較佳區域性的巢狀 loop:
 
        // for y:
        //   for x:
        //     f(x, y) = sin(x - y)
        //   for r:
        //     f(r, y) = (f(r, y) + f(r - 1, y)) / 2
        //   for r:
        //     f(7 - r, y) = (f(7 - r, y) + f(8 - r, y)) / 2
        //   for x:
        //     f_in_g(x, y) = f(x, y)
        // for y:
        //   for x:
        //     g(x, y) = f_in_g((x + y) % 8, (x - y) % 8);
        // 能額外對 f 的初始作向量化, 並且將 pixel 數值自 f 帶給其包裝函式:
        f.vectorize(x, 4);
        f.in(g).vectorize(x, 4);

        g.realize(8, 8);
        // 下圖為視覺化結果

    }

    printf("Success!\n");

    return 0;
}

==

Lesson 20: 複製函式(Cloning Funcs)


// 這節示範如何使用 Func::clone_in 來建構一個函式的複製

#include "Halide.h"
#include <stdio.h>

using namespace Halide;

int main(int argc, char **argv) {
    // 首先宣告後續會用到的變數
    Var x("x"), y("y"), xo("xo"), yo("yo"), xi("xi"), yi("yi");

    // 這結是關於使用 Func::clone_in 指令來複製一個 Func 函式
    {
        // 考量一個簡單二階的 pipeline:
        Func f("f_single"), g("g_single"), h("h_single");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y) + 3;
        h(x, y) = f(x, y) + g(x, y) + 10;

        f.compute_root();
        g.compute_root();
        h.compute_root();

        // 這產生下列巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     g(x, y) = 2 * f(x, y) + 3
        // for y:
        //   for x:
        //     h(x, y) = f(x, y) + g(x, y) + 10

        // 使用 Func::clone_in, 能夠以獨立排程的 'f' 複製 取代 'g' 中的 'f' 呼叫:
        Func f_clone_in_g = f.clone_in(g);
        f_clone_in_g.compute_root();

        // 同樣地, 能夠以串聯方式來排程:
        // f.clone_in(g).compute_root();

        // 這產生下列巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     f_clone_in_g(x, y) = x + y
        // for y:
        //   for x:
        //     g(x, y) = 2 * f_clone_in_g(x, y) + 3
        // for y:
        //   for x:
        //     h(x, y) = f(x, y) + g(x, y) + 10

        h.realize(5, 5);

        // 排程指令 f.clone_in(g) 以一個 'f' 的複製取代 'g' 中所有 'f' 的呼叫, 並回傳此複製
        // 本質上來說, 它以下列方式改寫上方原有的 pipeline:
        {
            Func f_clone_in_g("f_clone_in_g"), f("f"), g("g"), h("h");
            f(x, y) = x + y;
            f_clone_in_g(x, y) = x + y;
            g(x, y) = 2 * f_clone_in_g(x, y) + 3;
            h(x, y) = f(x, y) + g(x, y) + 10;

            f.compute_root();
            f_clone_in_g.compute_root();
            g.compute_root();
            h.compute_root();
        }
    }

    {
        // 在上面的排程只有 g 所產生的 'f' 的呼叫被取代
        // 其他對於 'f' 的呼叫依然直接呼叫 'f' 函式. (這裡也就是. 'h' 依然呼叫 'f' 而非複製)
        // 若期望以一個複製取代 'g' 與 'h' 對於 'f' 的呼叫, 只要呼叫 f.clone_in({g, h}).

        // 考量一個三階 pipeline, 其中兩階參考使用到 f:
        Func f("f_group"), g("g_group"), h("h_group"), out("out_group");
        f(x, y) = x + y;
        g(x, y) = 2 * f(x, y);
        h(x, y) = f(x, y) + 10;
        out(x, y) = f(x, y) + g(x, y) + h(x, y);

        f.compute_root();
        g.compute_root();
        h.compute_root();
        out.compute_root();

        // 以一個複製來取代 'g' 與 'h' 中所有對 'f' 的呼叫:
        f.clone_in({g, h}).compute_root();

        // 等義的巢狀 loop:
        // for y:
        //   for x:
        //     f(x, y) = x + y
        // for y:
        //   for x:
        //     f_clone(x, y) = x + y
        // for y:
        //   for x:
        //     g(x, y) = 2 * f_clone(x, y)
        // for y:
        //   for x:
        //     h(x, y) = f_clone(x, y) + 10
        // for y:
        //   for x:
        //     out(x, y) = f(x, y) + g(x, y) + h(x, y)

        out.realize(5, 5);
    }

    {
        // 一個應用 Func::clone_in() 的情況是當兩個 consumers 使用相當的 producer
        // 然而所需的範圍有非常的區隔. 考慮以下情況的範例:
        Func f("f"), g("g"), h("h");
        f(x) = x;
        g(x) = 2 * f(0);
        h(x) = f(99) + 10;

        // 將 'f' 以 compute_root 方式排程.
        f.compute_root();
        // 由於 'g' 與 'f' 都參考使用了 'f',  需要的 'f' 在 x 維度的範圍是 [0, 99]
        // 等義的巢狀 loop 為:
        // for x = 0 to 99
        //   f(x) = x
        // for x:
        //   g(x) = 2 * f(0)
        // for x:
        //   h(x) = f(99) + 10

        // 若 'f' 的計算成本相當昂貴, 那麼 'g' 與 'h' 最好能有著自己的一份 'f',
        // 來避免不必要的 'f' 計算.
        // 透過下列來對 'g' 與 'h' 建構區隔的 'f' 複製:
        f.clone_in(g).compute_root(); 
        // 等義的巢狀 loop 為:
        // f(0) = x
        // f_clone(99) = x
        // for x:
        //   g(x) = 2 * f_clone(0)
        // for x:
        //   h(x) = f(99) + 10
    }

    printf("Success!\n");

    return 0;
}
 
==

Lesson 21: 自動排程器(Auto-Scheduler) (Part 1, Part 2)


// 至今都是手工撰寫 Halide 排程, 但是有可能要求 Halide 建議一個合理的排程.
// 這稱為自動排程.
// 這節示範如何使用自動排程器來產生一個可剪貼套用並能後續改善的 CPU 排程
#include "Halide.h"
#include <stdio.h>

using namespace Halide;

// 定義一個自動排成器
class AutoScheduled : public Halide::Generator {
public:
    Input<Buffer<float>>  input{"input", 3};
    Input<float>          factor{"factor"};

    Output<Buffer<float>> output1{"output1", 2};
    Output<Buffer<float>> output2{"output2", 2};

    Expr sum3x3(Func f, Var x, Var y) {
        return f(x-1, y-1) + f(x-1, y) + f(x-1, y+1) +
               f(x, y-1)   + f(x, y)   + f(x, y+1) +
               f(x+1, y-1) + f(x+1, y) + f(x+1, y+1);
    }

    void generate() {
        // 這裡使用 Harris corner detection 作為演算
        Func in_b = BoundaryConditions::repeat_edge(input);

        gray(x, y) = 0.299f * in_b(x, y, 0) + 0.587f * in_b(x, y, 1) + 0.114f * in_b(x, y, 2);

        Iy(x, y) = gray(x-1, y-1)*(-1.0f/12) + gray(x-1, y+1)*(1.0f/12) +
                   gray(x, y-1)*(-2.0f/12) + gray(x, y+1)*(2.0f/12) +
                   gray(x+1, y-1)*(-1.0f/12) + gray(x+1, y+1)*(1.0f/12);

        Ix(x, y) = gray(x-1, y-1)*(-1.0f/12) + gray(x+1, y-1)*(1.0f/12) +
                   gray(x-1, y)*(-2.0f/12) + gray(x+1, y)*(2.0f/12) +
                   gray(x-1, y+1)*(-1.0f/12) + gray(x+1, y+1)*(1.0f/12);

        Ixx(x, y) = Ix(x, y) * Ix(x, y);
        Iyy(x, y) = Iy(x, y) * Iy(x, y);
        Ixy(x, y) = Ix(x, y) * Iy(x, y);
        Sxx(x, y) = sum3x3(Ixx, x, y);
        Syy(x, y) = sum3x3(Iyy, x, y);
        Sxy(x, y) = sum3x3(Ixy, x, y);
        det(x, y) = Sxx(x, y) * Syy(x, y) - Sxy(x, y) * Sxy(x, y);
        trace(x, y) = Sxx(x, y) + Syy(x, y);
        harris(x, y) = det(x, y) - 0.04f * trace(x, y) * trace(x, y);
        output1(x, y) = harris(x + 2, y + 2);
        output2(x, y) = factor * harris(x + 2, y + 2);
    }

    void schedule() {
        if (auto_schedule) {
            // 自動排程器需要依序評估所有 input/output 大小以及參數數值,
            // 來比較不同的選擇方案來決定一個好的排程.

            // 使用 set_bounds_estimate() 來提供輸入影像的每個維度評估 (最小以及延伸數值),
            // set_bounds_estimate() 帶入對應維度的 (最小, 延伸數值) 作為參數
            input.dim(0).set_bounds_estimate(0, 1024);
            input.dim(1).set_bounds_estimate(0, 1024);
            input.dim(2).set_bounds_estimate(0, 3);

            // 使用 set_estimate 來評估參數數值(parameter values).
            factor.set_estimate(2.0f);

            // 使用 estimate() 來評估 pipeline 輸出的每個維度(最小與延伸數值)
            // estimate() 帶入 (dim_name, 最小, 延伸數值) 作為參數
            output1.estimate(x, 0, 1024)
                   .estimate(y, 0, 1024);

            output2.estimate(x, 0, 1024)
                   .estimate(y, 0, 1024);

            // 技術上來說, 評估數值可以是任意值, 但愈接近應用數值形式, 產生的排程愈好.

            // 為了自動排程 pipeline, 不需要做其他事情:
            // 每個 Generator 潛在地有著一個名為 "auto_schedule" 的 GeneratorParam ;
            // 若將其設為 true, Halide 將對 pipeline 所有的輸出自動地呼叫 auto_schedule().

            // 每個 Generator 潛在地有著一個名為 "machine_params" 的 GeneratorParams ;
            // 這能允許指定計算架構的特對於自動排成器,  通常在 Makefile 中指定
            // 若沒有指定, 自動排成器將使用通用 CPU 的預設架構參數

            // (這段屬於架構專有解釋, 略譯).
            // Let's see some arbitrary but plausible values for the machine parameters.
            //
            //      const int kParallelism = 32;
            //      const int kLastLevelCacheSize = 16 * 1024 * 1024;
            //      const int kBalance = 40;
            //      MachineParams machine_params(kParallelism, kLastLevelCacheSize, kBalance);
            //
            // The arguments to MachineParams are the maximum level of parallelism
            // available, the size of the last-level cache (in KB), and the ratio
            // between the cost of a miss at the last level cache and the cost
            // of arithmetic on the target architecture, in that order.

            // 請注意當使用自動排程器時,  不能對 pipeline 套用任何排程, 否則將換產生錯誤
            // 目前自動排程器無法處理局部排程的 pipeline.

            // 當 HL_DEBUG_CODEGEN 設為 3 或更大數值, 排程將會透過 stdout 印出
            // 一個更好的方式是加入 "schedule" 到給 Generator 的 -e 旗標
            // (在 CMake 與 Bazel 這是透過使用 "extra_outputs" 旗標.)

            // 產生的排程為能夠直接剪貼小幅修改而使用的 Halide C++ 原始碼檔案
            // 這能作為起始的排程, 並再反覆地改善.
            // 請注意目前的自動排程器僅能夠產生 CPU 排程,
            // 且搭配 tiling, 簡單的向量化與平行化
            // 並不會處理 line buffer, 儲存次序調整或是歸納處理

            // 這個範例對於上述宣告的評估與架構參數, 自動排程器傳回下列排程:
            //
            // Var x_i("x_i");
            // Var x_i_vi("x_i_vi");
            // Var x_i_vo("x_i_vo");
            // Var x_o("x_o");
            // Var x_vi("x_vi");
            // Var x_vo("x_vo");
            // Var y_i("y_i");
            // Var y_o("y_o");
            //
            // Func f0 = pipeline.get_func(3);
            // Func f1 = pipeline.get_func(7);
            // Func f11 = pipeline.get_func(14);
            // Func f2 = pipeline.get_func(4);
            // Func output1 = pipeline.get_func(15);
            // Func output2 = pipeline.get_func(16);
            //
            // {
            //     Var x = f0.args()[0];
            //     f0
            //         .compute_at(f11, x_o)
            //         .split(x, x_vo, x_vi, 8)
            //         .vectorize(x_vi);
            // }
            // {
            //     Var x = f1.args()[0];
            //     f1
            //         .compute_at(f11, x_o)
            //         .split(x, x_vo, x_vi, 8)
            //         .vectorize(x_vi);
            // }
            // {
            //     Var x = f11.args()[0];
            //     Var y = f11.args()[1];
            //     f11
            //         .compute_root()
            //         .split(x, x_o, x_i, 256)
            //         .split(y, y_o, y_i, 128)
            //         .reorder(x_i, y_i, x_o, y_o)
            //         .split(x_i, x_i_vo, x_i_vi, 8)
            //         .vectorize(x_i_vi)
            //         .parallel(y_o)
            //         .parallel(x_o);
            // }
            // {
            //     Var x = f2.args()[0];
            //     f2
            //         .compute_at(f11, x_o)
            //         .split(x, x_vo, x_vi, 8)
            //         .vectorize(x_vi);
            // }
            // {
            //     Var x = output1.args()[0];
            //     Var y = output1.args()[1];
            //     output1
            //         .compute_root()
            //         .split(x, x_vo, x_vi, 8)
            //         .vectorize(x_vi)
            //         .parallel(y);
            // }
            // {
            //     Var x = output2.args()[0];
            //     Var y = output2.args()[1];
            //     output2
            //         .compute_root()
            //         .split(x, x_vo, x_vi, 8)
            //         .vectorize(x_vi)
            //         .parallel(y);
            // }
        } else {
            // 這裡能夠宣告先前手寫排程或是剪貼自動排程器所產生的排程.
            // 為了比較自動排程與基本的排程差異, 這裡使用了很陽春的排程
            gray.compute_root();
            Iy.compute_root();
            Ix.compute_root();
        }
    }
private:
    Var x{"x"}, y{"y"}, c{"c"};
    Func gray, Iy, Ix, Ixx, Iyy, Ixy, Sxx, Syy, Sxy, det, trace, harris;
};

// 如同在 lesson 15, 註冊 generator 後搭配 tools/GenGen.cpp 來一同編譯
HALIDE_REGISTER_GENERATOR(AutoScheduled, auto_schedule_gen)

// 編譯之後來看如何使用

// 這份 code 實際使用編譯出的 pipeline, 因此不會使用到 libHalide,
// 因此無需引入 Halide.h.
//
// 相對地, 這使用了上面的程式所產生的 header files
#include "auto_schedule_false.h"
#include "auto_schedule_true.h"

// 使用 Halide::Runtime::Buffer 類別來傳遞 pipeline 的資料輸入輸出
#include "HalideBuffer.h"
#include "halide_benchmark.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

int main(int argc, char **argv) {
    // 宣告與初始輸入影像
    Halide::Runtime::Buffer<float> input(1024, 1024, 3);

    for (int c = 0; c < input.channels(); ++c) {
        for (int y = 0; y < input.height(); ++y) {
            for (int x = 0; x < input.width(); ++x) {
                input(x, y, c) = rand();
            }
        }
    }

    Halide::Runtime::Buffer<float> output1(1024, 1024);
    Halide::Runtime::Buffer<float> output2(1024, 1024);
    // 為了評測而對每個版本 (沒有自動排程版本以及自動排程版本)執行數次
    double auto_schedule_off = Halide::Tools::benchmark(2, 5, [&]() {
        auto_schedule_false(input, 2.0f, output1, output2);
    });
    printf("Manual schedule: %gms\n", auto_schedule_off * 1e3);

    double auto_schedule_on = Halide::Tools::benchmark(2, 5, [&]() {
        auto_schedule_true(input, 2.0f, output1, output2);
    });
    printf("Auto schedule: %gms\n", auto_schedule_on * 1e3);

    // auto_schedule_on 應該會比較快, 因為 auto_schedule_off 版本的排程很簡陋,
    assert(auto_schedule_on < auto_schedule_off);

    return 0;
}

Halide's FFT 實作與效能比較

在許多 Halide 相關論文的投影片都可以看到與 FFTW 的速度比較 宣稱以 Halide 實作的 FFT 相較為 FFTW 的 1.5X 快 像是 An Introduction to Halide (CVPR 2015 以及 Halide: Decoupling...