2018年2月26日 星期一

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

Part 7 是記憶體格式與歸納範圍的進階技巧(非矩型, 向量化與平行化)
在實務上是用來處理複雜的計算規則, 以及對歸納的加速方式

==

Lesson 16 - RGB 影像與記憶體格式 (Part 1, Part 2)


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

using namespace Halide;
// 定義一個用以調亮 RGB 影像的 generator
class Brighten : public Halide::Generator<Brighten> {
public:
    // 定義輸入為一個三維影像, 前兩個維度為座標 x, y
    // 第三個為 color channel
   Input<Buffer<uint8_t>> input{"input", 3};
    // 為了能接收不同記憶體格式的輸入與輸出,
    // 將會以許多不同方式編譯此 generator
    // 而這是 GeneratorParam 好用的地方 (見 lesson 15)
    enum class Layout { Planar, Interleaved, Either, Specialized };
    GeneratorParam<Layout> layout{"layout",            // 預設值
            Layout::Planar,            // 自名稱到數值的對應
                {{ "planar",        Layout::Planar },
                 { "interleaved",   Layout::Interleaved },
                 { "either",        Layout::Either },
                 { "specialized",   Layout::Specialized }}};

    // 定義一個純量輸入來控制調整亮度的程度
    Input<uint8_t> offset{"offset"};
    // 宣告輸出
    Output<Buffer<uint8_t>> brighter{"brighter", 3};
    // 宣告變數
    Var x, y, c;

    void generate() {
        // 定義函式
        brighter(x, y, c) = input(x, y, c) + offset;

        // 排程
        brighter.vectorize(x, 16);
        // 取決於 generator 中的 'layout' 參數,
        // 將會編譯此 pipeline 來處理不同的記憶體格式
        if (layout == Layout::Planar) {
            // 此 pipeline 處理每條 scanline 為單一顏色的緊密排列.
            // 以 Lesson 10 stride 方式的描述, 在 x 的 stride 為 1
            // 這限制只允許平面的影像, 紅綠藍 channl 以下列方式擺放:

            // RRRRRRRR
            // RRRRRRRR
            // RRRRRRRR
            // RRRRRRRR
            // GGGGGGGG
            // GGGGGGGG
            // GGGGGGGG
            // GGGGGGGG
            // BBBBBBBB
            // BBBBBBBB
            // BBBBBBBB
            // BBBBBBBB

            // 這也能使用在較少使用的逐行的格式
            // 每一個 scanline 的紅綠藍依序擺放

            // RRRRRRRR
            // GGGGGGGG
            // BBBBBBBB
            // RRRRRRRR
            // GGGGGGGG
            // BBBBBBBB
            // RRRRRRRR
            // GGGGGGGG
            // BBBBBBBB
            // RRRRRRRR
            // GGGGGGGG
            // BBBBBBBB

        } else if (layout == Layout::Interleaved) {
            // 另一個常見的格式為'interleaved',
            // 對於每個 pixel 的紅綠藍數值在記憶體中依序存放:

            // RGBRGBRGBRGBRGBRGBRGBRGB
            // RGBRGBRGBRGBRGBRGBRGBRGB
            // RGBRGBRGBRGBRGBRGBRGBRGB
            // RGBRGBRGBRGBRGBRGBRGBRGB

            // 在這情形中 x 的 stride 值為 3, y 的 stride 值為影像寬的 3 倍
            // 而 c 的 stride 值為 1.
            // 能夠藉由以下來告知 Halide 輸入輸出格式為此:

            input
                .dim(0).set_stride(3) // dimension 0 (x) 的 stride 值為 3
                .dim(2).set_stride(1); // dimension 2 (c) 的 stride 值為 1

            brighter
                .dim(0).set_stride(3)
                .dim(2).set_stride(1);
            // 對於交錯格式, 可能會期望使用不同的排程
            // 藉由告知 Halide 額外的認定並檢測有 3 個 color channel 的情況,
            // 並利用這特性來讓 c loop 在最內層並且將其展開

            input.dim(2).set_bounds(0, 3); // Dimension 2 (c) starts at 0 and has extent 3.
            brighter.dim(2).set_bounds(0, 3);
            // 將 color channel loop 移置於最內層, 並且將其展開
            brighter.reorder(c, x, y).unroll(c);
            // 請注意當處理具有 alpha channel 的影像時,
            // x 的 stride 數值與 channel bound 為 4, 而非原本的 3
        } else if (layout == Layout::Either) {
            // 亦能夠移除所有的限制來編譯產生能是用在任何記憶體格式上
            // 但如此可能會很慢, 因為所有向量的載入變成了 gather (彙集)
            // 而所有向量的寫入變成了 scatter (散佈)
            input.dim(0).set_stride(Expr()); // 使用預設建構方式
                                             // 未定義的 Expr 表示沒有限制

            brighter.dim(0).set_stride(Expr());

        } else if (layout == Layout::Specialized) {
            // 藉由在執行時期告知 Halide 預期的記憶體格式, 與偵測到的 stride
            // 就能夠以良好的效率接收任何記憶體格式.
            // 首先解除當 dim(0).stride == 1 時的預設限制:

            input.dim(0).set_stride(Expr()); // 使用未定義的 Expr  來表示沒有限制

            brighter.dim(0).set_stride(Expr());
            // 建構用以偵測在執行時期是否為平面或是交錯的布林表示.
            // 這些條件必須針對所有被使用的情形之下一一檢測
            Expr input_is_planar =
                (input.dim(0).stride() == 1);
            Expr input_is_interleaved =
                (input.dim(0).stride() == 3 &&
                 input.dim(2).stride() == 1 &&
                 input.dim(2).extent() == 3);

            Expr output_is_planar =
                (brighter.dim(0).stride() == 1);
            Expr output_is_interleaved =
                (brighter.dim(0).stride() == 3 &&
                 brighter.dim(2).stride() == 1 &&
                 brighter.dim(2).extent() == 3);

            // 能夠使用 Func::specialize 來撰寫在執行時期切換的排程
            // 相關特定程式碼基於一個布林表示 (boolean Expr)
            // 該程式碼上利用著已知為 true 的 Expr.

            brighter.specialize(input_is_planar && output_is_planar);
            // 目前已對 brighter 作向量化與平行化,
            // 而此兩個列舉會繼承那些排程指示
            // 此外也能加入額外只套用在單一列舉的排程指令
            // 藉此告知 Halide 產生列舉版本的程式碼,
            // 其內容為交錯的輸入格式, 需調整次序且展開 loop

            brighter.specialize(input_is_interleaved && output_is_interleaved)
                .reorder(c, x, y).unroll(c);
            // 此外也能針對輸入為交錯, 輸出為平面而加入列舉, 或者相反的情況
            // 但列舉兩個情況已經足夠展示這個功能特性
            // 後續將會探索 Func::specialize 更創新的用法
            // 增加列舉可能會根本上增進套用條件的效能
            // 但也同時增加了需要編譯與產生的程式碼
            // 若執行檔大小是個考量, 且輸入輸出是明確的已知,
            // 應用上或許會想改用 set_stride 與 set_extent 取代
        }
    }
};
// 如同在 lesson 15 中, 註冊 generator 並且接著搭配 tools/GenGen.cpp 編譯
HALIDE_REGISTER_GENERATOR(Brighten, brighten)

// 編譯後來看看如何使用
// 這是實際上如何使用已經編譯出的 Halide pipeline 範例
// 由於並不需要依賴 libHalide, 因此無需引入 Halide.h.
//
// 取而代之地, 需要自上述程式所產生的 header 檔
#include "brighten_planar.h"
#include "brighten_interleaved.h"
#include "brighten_either.h"
#include "brighten_specialized.h"
// 將使用 Halide::Runtime::Buffer 在 pipeline 中傳遞資料
#include "HalideBuffer.h"

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

#include "clock.h"
// 提供各情況下執行時間的輔助函式. 印出自上次呼叫此函式起算的間隔時間
double tick(const char *name) {
    static double t_old = 0;
    double t_new = current_time();
    double dt = t_new - t_old;
    if (name) {
        printf("%s: %f\n", name, dt);
    }
    t_old = t_new;
    return dt;
}

int main(int argc, char **argv) {
    // 產生以交錯與平面格式儲存的影像
    // Halide::Runtime::Buffer 的預設是平面格式
    Halide::Runtime::Buffer<uint8_t> planar_input(1024, 768, 3);
    Halide::Runtime::Buffer<uint8_t> planar_output(1024, 768, 3);
    Halide::Runtime::Buffer<uint8_t> interleaved_input =
        Halide::Runtime::Buffer<uint8_t>::make_interleaved(1024, 768, 3);
    Halide::Runtime::Buffer<uint8_t> interleaved_output =
        Halide::Runtime::Buffer::make_interleaved(1024, 768, 3);
    // 提供在設置 generator 時限制的數值來檢查 stride 數值是否如預期,
    assert(planar_input.dim(0).stride() == 1);
    assert(planar_output.dim(0).stride() == 1);
    assert(interleaved_input.dim(0).stride() == 3);
    assert(interleaved_output.dim(0).stride() == 3);
    assert(interleaved_input.dim(2).stride() == 1);
    assert(interleaved_output.dim(2).stride() == 1);
    // 接著呼叫所編譯出的不同函數並一一確認執行效能

    // 啟動時鐘
    tick(NULL);
    // 在平面影像上執行平面版本, 以及在交錯影像上執行交錯版本
    // 每個版本都會執行 1000 次來取得評測結果.
    for (int i = 0; i < 1000; i++) {
        brighten_planar(planar_input, 1, planar_output);
    }
    double planar_time = tick("brighten_planar");
    // Click to show output ...

    for (int i = 0; i < 1000; i++) {
        brighten_interleaved(interleaved_input, 1, interleaved_output);
    }
    double interleaved_time = tick("brighten_interleaved");
    // Click to show output ...

   // 多數的操作上平面格式通常會叫交錯格式來得快
    assert(planar_time < interleaved_time);
    // 因為 stride 並非 generator 中所設定的, 因此下列兩個註解掉的呼叫會產生錯誤

    // brighten_planar(interleaved_input, 1, interleaved_output);
    // Error: Constraint violated: brighter.stride.0 (3) == 1 (1)

    // brighten_interleaved(planar_input, 1, planar_output);
    // Error: Constraint violated: brighter.stride.0 (1) == 3 (3)

    // 執行彈性版本並檢測效能, 應會正常運作但比上面版本都慢
    for (int i = 0; i < 1000; i++) {
        brighten_either(planar_input, 1, planar_output);
    }
    double either_planar_time = tick("brighten_either on planar images");
    assert(planar_time < either_planar_time);
    // Click to show output ...

    for (int i = 0; i < 1000; i++) {
        brighten_either(interleaved_input, 1, interleaved_output);
    }
    double either_interleaved_time = tick("brighten_either on interleaved images");
    assert(interleaved_time < either_interleaved_time);
    // Click to show output ...

    // 針對不同的格式執行列舉的版本
    // 在內部切換到等義程式碼的作法應與對不同格式編譯的版本效能一致
    for (int i = 0; i < 1000; i++) {
        brighten_specialized(planar_input, 1, planar_output);
    }
    double specialized_planar_time = tick("brighten_specialized on planar images");
    // Click to show output ...

    // if 的描述的成本應該可以被忽略, 但這個測試允許 50% 為測量誤差
    assert(specialized_planar_time < 1.5 * planar_time);

    for (int i = 0; i < 1000; i++) {
        brighten_specialized(interleaved_input, 1, interleaved_output);
    }
    double specialized_interleaved_time = tick("brighten_specialized on interleaved images");
    assert(specialized_interleaved_time < 2.0 * interleaved_time);
    // Click to show output ...

    return 0;
}

==

Lesson 17: 在非長方區域做歸納

// 這節示範如何定義使用 predicate (推斷, 使...基於)在部份的歸納區域做更新

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

using namespace Halide;

int main(int argc, char **argv) {
    // 在 lesson 9, 已經學習過如何使用 RDom 去定義在 Halide 更新定義中使用的歸納區域
    // 然而, 以 RDom 定義的範圍通常是長方形, 而更新發生在這長方區域中的每個點
    // 在一些情況中, 可能會想在一些非長方區域做運算
    // 例如: 圓形, 能藉由使用 RDom::where 指令來達成
    {
         // 以此純定義開始:
        Func circle("circle");
        Var x("x"), y("y");
        circle(x, y) = x + y;
        // 這裡要在中心點為(3, 3) 半徑為 3 的圓範圍內更新數值為原有的平方
        // 為了達成此目的, 首先使用 RDom 定義能涵蓋此圓的最小方形
        RDom r(0, 7, 0, 7);
       // 這個方形並不需要最小的. 事實上這方形只要覆蓋要更新的圓形區域,
        // 可以是任意大小. 然而這個方形愈緊密, 產生的 loop 範圍也愈緊密
        // Halide 會自動地儘可能緊密 loop 的範圍, 但一般來說定義最小的方形會較好

        // 接著使用 RDom::where 來定義此方形上的推斷條件
        // 如此只會在推斷條件為真的區域上更新
        // 這表示在這圓形的區域
        r.where((r.x - 3)*(r.x - 3) + (r.y - 3)*(r.y - 3) <= 10);

        // 在定義推斷條件後, 接著定義更新動作
        circle(r.x, r.y) *= 2;

        Buffer<int> halide_result = circle.realize(7, 7);       
        // 下圖為視覺化結果



        // 等義的 C code:
        int c_result[7][7];
        for (int y = 0; y < 7; y++) {
            for (int x = 0; x < 7; x++) {
                c_result[y][x] = x + y;
            }
        }
        for (int r_y = 0; r_y < 7; r_y++) {
            for (int r_x = 0; r_x < 7; r_x++) {
                // Update is only performed if the predicate evaluates to true.
                if ((r_x - 3)*(r_x - 3) + (r_y - 3)*(r_y - 3) <= 10) {
                    c_result[r_y][r_x] *= 2;
                }
            }
        }

        // 檢查結果是否一致:
        for (int y = 0; y < 7; y++) {
            for (int x = 0; x < 7; x++) {
                if (halide_result(x, y) != c_result[y][x]) {
                    printf("halide_result(%d, %d) = %d instead of %d\n",
                           x, y, halide_result(x, y), c_result[y][x]);
                    return -1;
                }
            }
        }
    }

    {
        // 也能夠在 RDom 上定義多個推斷條件. 在此嘗試在三角形區域中做更新
        // 為了達到這目的, 定義三個推斷條件, 各對應為三角形的一個邊
        Func triangle("triangle");
        Var x("x"), y("y");
        triangle(x, y) = x + y;
        // 首先定義能涵蓋三角形範圍的最小方形
        RDom r(0, 8, 0, 10);
         // 接著藉由呼較多次 RDom::where 將三個推斷條件加入 RDom
        r.where(r.x + r.y > 5);
        r.where(3*r.y - 2*r.x < 15);
        r.where(4*r.x - r.y < 20);

        // 也能如下地能夠將多個推斷條件包裹成一個:
        // r.where((r.x + r.y > 5) && (3*r.y - 2*r.x < 15) && (4*r.x - r.y < 20));

        // 定義更新步驟.
        triangle(r.x, r.y) *= 2;

        Buffer<int> halide_result = triangle.realize(10, 10);

         // 下圖為視覺化結果

        // 等義 C code:
        int c_result[10][10];
        for (int y = 0; y < 10; y++) {
            for (int x = 0; x < 10; x++) {
                c_result[y][x] = x + y;
            }
        }
        for (int r_y = 0; r_y < 10; r_y++) {
            for (int r_x = 0; r_x < 8; r_x++) {
                // 更新僅會在推斷條件為 true 時執行
                if ((r_x + r_y > 5) && (3*r_y - 2*r_x < 15) && (4*r_x - r_y < 20)) {
                    c_result[r_y][r_x] *= 2;
                }
            }
        }

        // 確認結果一致:
        for (int y = 0; y < 10; y++) {
            for (int x = 0; x < 10; x++) {
                if (halide_result(x, y) != c_result[y][x]) {
                    printf("halide_result(%d, %d) = %d instead of %d\n",
                           x, y, halide_result(x, y), c_result[y][x]);
                    return -1;
                }
            }
        }
    }

    {
        // 推斷條件的使用並不只侷限於歸納範圍的變數 (r.x, r.y, ...)
        // 亦能夠參考使用其他更新定義上的自由變數
        // 甚至是呼叫使用其他 Funcs 或是遞迴呼叫一樣的 Func 像是:
        Func f("f"), g("g");
        Var x("x"), y("y");
        f(x, y) = 2 * x + y;
        g(x, y) = x + y;

        // 這個歸納範圍的推斷條件取決於 'f' 的初始值.
        RDom r1(0, 5, 0, 5);
        r1.where(f(r1.x, r1.y) >= 4);
        r1.where(f(r1.x, r1.y) <= 7);
        f(r1.x, r1.y) /= 10;

        f.compute_root();

        // 這個例子使用其他 Func 呼叫.
        RDom r2(1, 3, 1, 3);
        r2.where(f(r2.x, r2.y) < 1);
        g(r2.x, r2.y) += 17;

        Buffer<int> halide_result_g = g.realize(5, 5);
        // 下圖為視覺化結果

        // 'f' 部份等義的 C code:
        int c_result_f[5][5];
        for (int y = 0; y < 5; y++) {
            for (int x = 0; x < 5; x++) {
                c_result_f[y][x] = 2 * x + y;
            }
        }
        for (int r1_y = 0; r1_y < 5; r1_y++) {
            for (int r1_x = 0; r1_x < 5; r1_x++) {
                // 當推斷條件為 true 時才會執行更新
                if ((c_result_f[r1_y][r1_x] >= 4) && (c_result_f[r1_y][r1_x] <= 7)) {
                    c_result_f[r1_y][r1_x] /= 10;
                }
            }
        }

        // 'g' 的等義 C code:
        int c_result_g[5][5];
        for (int y = 0; y < 5; y++) {
            for (int x = 0; x < 5; x++) {
                c_result_g[y][x] = x + y;
            }
        }
        for (int r2_y = 1; r2_y < 4; r2_y++) {
            for (int r1_x = 1; r1_x < 4; r1_x++) {
                // 當推斷條件為 true 時才會執行更新
                if (c_result_f[r2_y][r1_x] < 1) {
                    c_result_g[r2_y][r1_x] += 17;
                }
            }
        }

        // 確認結果是否一致:
        for (int y = 0; y < 5; y++) {
            for (int x = 0; x < 5; x++) {
                if (halide_result_g(x, y) != c_result_g[y][x]) {
                    printf("halide_result_g(%d, %d) = %d instead of %d\n",
                           x, y, halide_result_g(x, y), c_result_g[y][x]);
                    return -1;
                }
            }
        }
    }

    printf("Success!\n");

    return 0;
}

==

Lesson 18: Factoring an associative reduction using rfactor


// 這節示範如何使用排程指令 'rfactor' 來對結合歸納作平行化或向量化

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

using namespace Halide;

int main(int argc, char **argv) {
    // 宣告後續會使用到的變數
    Var x("x"), y("y"), i("i"), u("u"), v("v");

    // 以亂數建立輸入
    Buffer input(8, 8, "input");
    for (int y = 0; y < 8; ++y) {
        for (int x = 0; x < 8; ++x) {
            input(x, y) = (rand() % 256);
        }
    }

    {
        // 如同先前在 Lesson 9 提到的, 對歸納範圍部份的變數平行化是難以處理的,
        // 因為資料相依性可能橫跨多個變數

        // 考量 Lesson 9 中的 histogram 範例:
        Func histogram("hist_serial");
        histogram(i) = 0;
        RDom r(0, input.width(), 0, input.height());
        histogram(input(r.x, r.y) / 32) += 1;

        histogram.vectorize(i, 8);
        histogram.realize(8);
        // 下圖為視覺化結果
       
        // 能夠對初始化 histogram bucket 動作做向量化,
        // 由於在更新定義中對於有著對 x, y 方向的資料相依性
        // (換句話說更新參考了先前更新所計算出的數值),
        // 因此無在不產生 race condition 下對 r.x 或 r.y 作平行化或是向量化
        // 下列程式碼將產生一個錯誤:
        // histogram.update().parallel(r.y);
    }

    {
        // 然而, 請注意 histogram 是個俱備結合特性的操作 (是一種總和歸納)
        // 一個通常用來加速結合性歸納的技巧是分成多個部份,
        // 各別計算各自的部份結果, 最後將這接部份結果整併.
        // 由於每個分割的部份的計算是獨立的, 因此能夠對其做平行化
        // 回到 histogram 的例子, 藉由定義中介函數獨立計算每列的 histogram
        // 以此將歸納範圍以數列的方式分割:
        Func intermediate("intm_par_manual");
        intermediate(i, y) = 0;
        RDom rx(0, input.width());
        intermediate(input(rx, y) / 32, y) += 1;

        // 接著定義用以累加局部結果的第二個 stage:
        Func histogram("merge_par_manual");
        histogram(i) = 0;
        RDom ry(0, input.height());
        histogram(i) += intermediate(i, ry);
        // 由於中介資訊不再有著 y 方向上的資料相依性
        // 因此能夠對 y 做平行化:
        intermediate.compute_root().update().parallel(y);

        // 依然能夠對初始化作向量化
        intermediate.vectorize(i, 8);
        histogram.vectorize(i, 8);

        histogram.realize(8);

        // 下圖為視覺化結果

    }

    {
        // 對於結合歸納的手動分解是枯燥且容易出錯的.
        // 儘管對於 histogram 而言手動處理相當簡單,
        // 但若 RDom 有著推斷條件(RDom::where), 或當函數歸納到多維 tuple 上,
        // 情況將會相當快速地變得複雜
        // Halide 提供了一個透過排程指令 'rfactor' 的方式來處理此類的分解.
        // rfactor 將具結合性更新定義分為兩部份
        // 首先是一個用以計算部份歸納區域結果的中介
        // 接著以整併局部結果的新定義來取代既有的更新定義

        // 使用 rfactor, 一點都無需改變演算法:
        Func histogram("hist_rfactor_par");
        histogram(x) = 0;
        RDom r(0, input.width(), 0, input.height());
        histogram(input(r.x, r.y) / 32) += 1;
        // 分解的動作透過 rfactor 而被移至排程中.
        // rfactor 的輸入為包含可平行化的歸納變數(RVars)的成對清單.


        // 在產生的中介函式 Func, 所有參考這些歸納變數的部份將以"純"變數取代.
        // 由於藉由建構, 變數不會造成 race condition,
        // 因此中介歸納在這些維度上為可平行化.
        // 不在此清單的所有歸納變數將自原定義移除, 並移到中介函式中

        // 下列程式能產生同等於手動分解版本:
        Func intermediate = histogram.update().rfactor({{r.y, y}});
        // 將 {r.y, y} 作為 rfactor 參數傳入,
        // 如此使得 histogram 類似於手動分解版本在 y 維度上可平行處理
        intermediate.compute_root().update().parallel(y);
        // 在這個情況中, 僅對一個變數維度做分割,
        // 能夠將大括號省略, 直接以下列方式撰寫 rfactor
        // Func intermediate = histogram.update().rfactor(r.y, y);
        // 如同先前一般, 對初始作向量化
        intermediate.vectorize(x, 8);
        histogram.vectorize(x, 8);
        // 務必注意到 rfactor(一般來說, 或是歸納的分解) 只對於結合性歸納有效.
        // 結合性歸納有著無論計算如何分群結果都一致的良好特性
        // 若 rfactor 無法證明歸納的結合特性, 將會丟出錯誤.

        Buffer halide_result = histogram.realize(8);

        // 下圖為視覺化結果

        // 等義的 C code:
        int c_intm[8][8];
        for (int y = 0; y < input.height(); y++) {
            for (int x = 0; x < 8; x++) {
                c_intm[y][x] = 0;
            }
        }
        /* parallel */ for (int y = 0; y < input.height(); y++) {
            for (int r_x = 0; r_x < input.width(); r_x++) {
                c_intm[y][input(r_x, y) / 32] += 1;
            }
        }

        int c_result[8];
        for (int x = 0; x < 8; x++) {
            c_result[x] = 0;
        }
        for (int x = 0; x < 8; x++) {
            for (int r_y = 0; r_y < input.height(); r_y++) {
                c_result[x] += c_intm[r_y][x];
            }
        }

        // 確認結果一致:
        for (int x = 0; x < 8; x++) {
            if (c_result[x] != halide_result(x)) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    {
        // 現在能夠透過排程指令 'rfactor' 分解結合性歸納
        // 使用此排程能探究多種的分解策略
        // 以同樣的 histogram 為例:
        Func histogram("hist_rfactor_vec");
        histogram(x) = 0;
        RDom r(0, input.width(), 0, input.height());
        histogram(input(r.x, r.y) / 32) += 1;

        // 相較先前使用 r.y, 這次改以 r.x 套用在 rfactor 上來對直行做分解
        Func intermediate = histogram.update().rfactor(r.x, u);
        // 現在對每行去獨立計算 histogram
        // 如此能在直行上作向量化
        intermediate.compute_root().update().vectorize(u, 8);
        // 請注意因為對內層維度向量化改變了數值累加到最後 histogram 的計算次序,
        // 因此這個技巧僅適用在俱備結合性(associative)與交換性(commutative)的歸納上
        // rfactor 將嘗試證明具有這些特性, 若無法則會丟出錯誤

        // 對初始做向量化
        intermediate.vectorize(x, 8);
        histogram.vectorize(x, 8);

        Buffer<int> halide_result = histogram.realize(8);

        // 下圖為視覺化結果


        // 等義的 C code:
        int c_intm[8][8];
        for (int u = 0; u < input.width(); u++) {
            for (int x = 0; x < 8; x++) {
                c_intm[u][x] = 0;
            }
        }
        for (int r_y = 0; r_y < input.height(); r_y++) {
            for (int u = 0; u < input.width() / 8; u++) {
                /* vectorize */ for (int u_i = 0; u_i < 8; u_i++) {
                    c_intm[u*4 + u_i][input(u*8 + u_i, r_y) / 32] += 1;
                }
            }
        }

        int c_result[8];
        for (int x = 0; x < 8; x++) {
            c_result[x] = 0;
        }
        for (int x = 0; x < 8; x++) {
            for (int r_x = 0; r_x < input.width(); r_x++) {
                c_result[x] += c_intm[r_x][x];
            }
        }

        // 檢查結果是否一致:
        for (int x = 0; x < 8; x++) {
            if (c_result[x] != halide_result(x)) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    {
        // 也能夠一次對多維度的歸納範圍做分割
        // 這次將以 tile 方式對範圍做分割並計算部分 histogram
        Func histogram("hist_rfactor_tile");
        histogram(x) = 0;
        RDom r(0, input.width(), 0, input.height());
        histogram(input(r.x, r.y) / 32) += 1;

        // 首先以 4 為單位同時分割 r.x 與 r.y .
        RVar rx_outer("rx_outer"), rx_inner("rx_inner");
        RVar ry_outer("ry_outer"), ry_inner("ry_inner");
        histogram.update()
            .split(r.x, rx_outer, rx_inner, 4)
            .split(r.y, ry_outer, ry_inner, 4);

        // 接著呼叫 rfactor 來產生中介函式來獨立計算每個 tile 的 histogram
        Func intermediate = histogram.update().rfactor({{rx_outer, u}, {ry_outer, v}});

        // 現在對每個 tile 的中介做平行化
        intermediate.compute_root().update().parallel(u).parallel(v);

        // 亦能夠調整最外層 tile 索引的次序來做經典的 tile 尋訪
        intermediate.update().reorder(rx_inner, ry_inner, u, v);

        // 對初始做向量化
        intermediate.vectorize(x, 8);
        histogram.vectorize(x, 8);

        Buffer halide_result = histogram.realize(8);

        // 下圖為視覺化結果

        // 等義的 Code :
        int c_intm[4][4][8];
        for (int v = 0; v < input.height() / 2; v++) {
            for (int u = 0; u < input.width() / 2; u++) {
                for (int x = 0; x < 8; x++) {
                    c_intm[v][u][x] = 0;
                }
            }
        }
        /* parallel */ for (int v = 0; v < input.height() / 2; v++) {
            /* parallel */ for (int u = 0; u < input.width() / 2; u++) {
                for (int ry_inner = 0; ry_inner < 2; ry_inner++) {
                    for (int rx_inner = 0; rx_inner < 2; rx_inner++) {
                        c_intm[v][u][input(u*2 + rx_inner, v*2 + ry_inner) / 32] += 1;
                    }
                }
            }
        }

        int c_result[8];
        for (int x = 0; x < 8; x++) {
            c_result[x] = 0;
        }
        for (int x = 0; x < 8; x++) {
            for (int ry_outer = 0; ry_outer < input.height() / 2; ry_outer++) {
                for (int rx_outer = 0; rx_outer < input.width() / 2; rx_outer++) {
                    c_result[x] += c_intm[ry_outer][rx_outer][x];
                }
            }
        }

        // 確認結果是否一致:
        for (int x = 0; x < 8; x++) {
            if (c_result[x] != halide_result(x)) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    printf("Success!\n");

    return 0;
}

沒有留言:

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

標題為筆記, 但這篇比較屬於心得 延續 上一篇 的環境建立, 這次計劃藉由 Jserv 最新的 課程安排 來學習 Chisel, 當然個人目標是能夠按照 Jserv 的課程規劃在 期限之內 完成 Lab 3, 由於個人並非 digital designer (現在這年紀也算老貓學...