Halide学习笔记----Halide tutorial源码阅读5

Halide入门教程05

// Halide教程第五课：向量化，并行化，平铺，数据分块
// 本课展示了如何才操作函数像素索引的计算顺序，包括向量化/并行化/平铺/分块等技术

// 在linux系统中，采用如下指令编译并执行
// g++ lesson_05*.cpp -g -I ../include -L ../bin -lHalide -lpthread -ldl -o lesson_05 -std=c++11
// LD_LIBRARY_PATH=../bin ./lesson_05

#include "Halide.h"
#include <stdio.h>
#include <algorithm>
using namespace Halide;

int main(int argc, char **argv) {

Var x("x"), y("y");

// First we observe the default ordering.
{
Func gradient("gradient");
gradient(x, y) = x + y;
gradient.trace_stores();

//默认遍历像素的顺序是行优先，即内层循环沿着行方向，外层循环沿着列方向
printf("Evaluating gradient row-major\n");
Buffer<int> output = gradient.realize(4, 4);

// The equivalent C is:
printf("Equivalent C:\n");
for (int y = 0; y < 4; y++) {
for (int x = 0; x < 4; x++) {
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
printf("\n\n");

// 跟踪系统调度可以很容易理解调度系统如何工作。可以通过Halide提供的函数来打印出实际工作
// 是执行的哪种循环调度。
printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");

// Because we're using the default ordering, it should print:
// compute gradient:
//   for y:
//     for x:
//       gradient(...) = ...
}

// Reorder variables.
{
Func gradient("gradient_col_major");
gradient(x, y) = x + y;
gradient.trace_stores();

// 可以通过reorder函数来改变函数遍历的顺序，下面的语句将行方向（y）置于内层循环，而将原本的内层
// 循环调整到了外循环。也就是说y遍历比x遍历更快。是一种列优先的遍历方法
gradient.reorder(y, x);

printf("Evaluating gradient column-major\n");
Buffer<int> output = gradient.realize(4, 4);

printf("Equivalent C:\n");
for (int x = 0; x < 4; x++) {
for (int y = 0; y < 4; y++) {
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
printf("\n");

//
printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Split a variable into two.
{
Func gradient("gradient_split");
gradient(x, y) = x + y;
gradient.trace_stores();

// 原始调度中，最有效的就是split调度了，它将一个大循环，拆解成一个外部循环和一个内部循环；
// 即，将x方向的循环，拆成一个外部循环x_outer和一个内部循环x_inner
// 下面的split将x拆成x_outer,x_inner, 内循环的长度为2
Var x_outer, x_inner;
gradient.split(x, x_outer, x_inner, 2);

printf("Evaluating gradient with x split into x_outer and x_inner \n");
Buffer<int> output = gradient.realize(4, 4);

printf("Equivalent C:\n");
for (int y = 0; y < 4; y++) {
for (int x_outer = 0; x_outer < 2; x_outer++) {
for (int x_inner = 0; x_inner < 2; x_inner++) {
int x = x_outer * 2 + x_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Fuse two variables into one.
{
Func gradient("gradient_fused");
gradient(x, y) = x + y;

// 和split相反的是fuse，它将两个变量融合成一个变量，fuse的重要性并没有split高。
Var fused;
gradient.fuse(x, y, fused);

printf("Evaluating gradient with x and y fused\n");
Buffer<int> output = gradient.realize(4, 4);

printf("Equivalent C:\n");
for (int fused = 0; fused < 4*4; fused++) {
int y = fused / 4;
int x = fused % 4;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Evaluating in tiles.
// tile的中文意思是瓦片，在这里是指将图像数据拆分成和瓦片一项的小图像块
{
Func gradient("gradient_tiled");
gradient(x, y) = x + y;
gradient.trace_stores();

// 既然我们可以拆分和调整顺序，我们可以按照划分数据块的方式来进行计算。将x和y方向拆分，然后
// 调制x和y的顺序，按照小的数据块的方式来进行遍历。
// 一个小的数据块将整个图像划分成小的矩形，外层循环在tile击毙恩进行循环，遍历所有的tile。
Var x_outer, x_inner, y_outer, y_inner;
gradient.split(x, x_outer, x_inner, 4);
gradient.split(y, y_outer, y_inner, 4);
gradient.reorder(x_inner, y_inner, x_outer, y_outer);

// This pattern is common enough that there's a shorthand for it:
// gradient.tile(x, y, x_outer, y_outer, x_inner, y_inner, 4, 4);

printf("Evaluating gradient in 4x4 tiles\n");
Buffer<int> output = gradient.realize(8, 8);

printf("Equivalent C:\n");
for (int y_outer = 0; y_outer < 2; y_outer++) {
for (int x_outer = 0; x_outer < 2; x_outer++) {
for (int y_inner = 0; y_inner < 4; y_inner++) {
for (int x_inner = 0; x_inner < 4; x_inner++) {
int x = x_outer * 4 + x_inner;
int y = y_outer * 4 + y_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Evaluating in vectors.
{
Func gradient("gradient_in_vectors");
gradient(x, y) = x + y;
gradient.trace_stores();

// split能够让内层循环变量在一个划分银子内变化。这个划分因子通常是指定的，因此在编译时是一个常数
// 因此我们可以调用向量化指令来执行内部循环。在这里我们指定这个因子为4，这样就可以调用x86机器上的SSE
//指令来计算4倍宽的向量，这里充分利用了cpu的SIMD指令来加快计算
Var x_outer, x_inner;
gradient.split(x, x_outer, x_inner, 4);
gradient.vectorize(x_inner);

// 上述过程有更简单的形式
// gradient.vectorize(x, 4);
// 等价于
// gradient.split(x, x, x_inner, 4);
// gradient.vectorize(x_inner);
// 这里我们重用了x，将它当作外循环变量，稍后的调度将x当作外循环(x_outer)来进行调度

// 这次在一个8x4的矩形上执行gradient算法
printf("Evaluating gradient with x_inner vectorized \n");
Buffer<int> output = gradient.realize(8, 4);

printf("Equivalent C:\n");
for (int y = 0; y < 4; y++) {
for (int x_outer = 0; x_outer < 2; x_outer++) {
// The loop over x_inner has gone away, and has been
// replaced by a vectorized version of the
// expression. On x86 processors, Halide generates SSE
// for all of this.
int x_vec[] = {x_outer * 4 + 0,
x_outer * 4 + 1,
x_outer * 4 + 2,
x_outer * 4 + 3};
int val[] = {x_vec[0] + y,
x_vec[1] + y,
x_vec[2] + y,
x_vec[3] + y};
printf("Evaluating at <%d, %d, %d, %d>, <%d, %d, %d, %d>:"
" <%d, %d, %d, %d>\n",
x_vec[0], x_vec[1], x_vec[2], x_vec[3],
y, y, y, y,
val[0], val[1], val[2], val[3]);
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Unrolling a loop.
{
Func gradient("gradient_unroll");
gradient(x, y) = x + y;
gradient.trace_stores();

// 如果多个像素共享一些重复的（overlapping）数据，可以将循环铺平，从而共享的数据只需要载入或者
// 计算一次。它和向量化的表达方式类似。先将数据进行划分，然后将内层循环铺平。
Var x_outer, x_inner;
gradient.split(x, x_outer, x_inner, 2);
gradient.unroll(x_inner);

// The shorthand for this is:
// gradient.unroll(x, 2);

printf("Evaluating gradient unrolled by a factor of two\n");
Buffer<int> result = gradient.realize(4, 4);

printf("Equivalent C:\n");
for (int y = 0; y < 4; y++) {
for (int x_outer = 0; x_outer < 2; x_outer++) {
// Instead of a for loop over x_inner, we get two
// copies of the innermost statement.
{
int x_inner = 0;
int x = x_outer * 2 + x_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
{

e8e5
int x_inner = 1;
int x = x_outer * 2 + x_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// Splitting by factors that don't divide the extent.
{
Func gradient("gradient_split_7x2");
gradient(x, y) = x + y;
gradient.trace_stores();

// 当原来图像尺寸不能整除划分的小矩形尺寸时，最后的一行或者一列的tile在边界处会出现重复计算的现象
Var x_outer, x_inner;
gradient.split(x, x_outer, x_inner, 3);

printf("Evaluating gradient over a 7x2 box with x split by three \n");
Buffer<int> output = gradient.realize(7, 2);

printf("Equivalent C:\n");
for (int y = 0; y < 2; y++) {
for (int x_outer = 0; x_outer < 3; x_outer++) { // Now runs from 0 to 2
for (int x_inner = 0; x_inner < 3; x_inner++) {
int x = x_outer * 3;
// Before we add x_inner, make sure we don't
// evaluate points outside of the 7x2 box. We'll
// clamp x to be at most 4 (7 minus the split
// factor).
if (x > 4) x = 4;
x += x_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");

// 如果仔细查看程序的输出，你会发现有些像素点进行了不止一次计算。这是因为尺寸不能整除所致
// 由于Halide函数没有边缘效应，因此计算多次并不会产生副作用。
}

// Fusing, tiling, and parallelizing.
{
// 这里才是fuse真正发挥威力的地方。如果想要在多个维度进行并行计算，
// 可以将多个循环网fuse起来，然后在fuse后的维度下进行并行计算

Func gradient("gradient_fused_tiles");
gradient(x, y) = x + y;
gradient.trace_stores();

Var x_outer, y_outer, x_inner, y_inner, tile_index;
gradient.tile(x, y, x_outer, y_outer, x_inner, y_inner, 4, 4);
gradient.fuse(x_outer, y_outer, tile_index);
gradient.parallel(tile_index);

// 每个调度函数返回的是引用类型，因此可以按如下方式用点号连接起多次调用
// gradient
//     .tile(x, y, x_outer, y_outer, x_inner, y_inner, 2, 2)
//     .fuse(x_outer, y_outer, tile_index)
//     .parallel(tile_index);

printf("Evaluating gradient tiles in parallel\n");
Buffer<int> output = gradient.realize(8, 8);

// tile层面的调度是乱序的，但是在每一个tile内部，是行优先的计算顺序

printf("Equivalent (serial) C:\n");
// This outermost loop should be a parallel for loop, but that's hard in C.
for (int tile_index = 0; tile_index < 4; tile_index++) {
int y_outer = tile_index / 2;
int x_outer = tile_index % 2;
for (int y_inner = 0; y_inner < 4; y_inner++) {
for (int x_inner = 0; x_inner < 4; x_inner++) {
int y = y_outer * 4 + y_inner;
int x = x_outer * 4 + x_inner;
printf("Evaluating at x = %d, y = %d: %d\n", x, y, x + y);
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient.print_loop_nest();
printf("\n");
}

// 将前面演示的调度综合在一起
{
Func gradient_fast("gradient_fast");
gradient_fast(x, y) = x + y;

// tile尺寸为64x64，采用并行计算
Var x_outer, y_outer, x_inner, y_inner, tile_index;
gradient_fast
.tile(x, y, x_outer, y_outer, x_inner, y_inner, 64, 64)
.fuse(x_outer, y_outer, tile_index)
.parallel(tile_index);

// 将内部的64x64tile继续拆分成更小的tile。在x方向上采用向量化的计算（调用SIMD指令），
// 在y方向进行平铺
Var x_inner_outer, y_inner_outer, x_vectors, y_pairs;
gradient_fast
.tile(x_inner, y_inner, x_inner_outer, y_inner_outer, x_vectors, y_pairs, 4, 2)
.vectorize(x_vectors)
.unroll(y_pairs);

Buffer<int> result = gradient_fast.realize(350, 250);

printf("Checking Halide result against equivalent C...\n");
for (int tile_index = 0; tile_index < 6 * 4; tile_index++) {
int y_outer = tile_index / 4;
int x_outer = tile_index % 4;
for (int y_inner_outer = 0; y_inner_outer < 64/2; y_inner_outer++) {
for (int x_inner_outer = 0; x_inner_outer < 64/4; x_inner_outer++) {
// We're vectorized across x
int x = std::min(x_outer * 64, 350-64) + x_inner_outer*4;
int x_vec[4] = {x + 0,
x + 1,
x + 2,
x + 3};

// And we unrolled across y
int y_base = std::min(y_outer * 64, 250-64) + y_inner_outer*2;
{
// y_pairs = 0
int y = y_base + 0;
int y_vec[4] = {y, y, y, y};
int val[4] = {x_vec[0] + y_vec[0],
x_vec[1] + y_vec[1],
x_vec[2] + y_vec[2],
x_vec[3] + y_vec[3]};

// Check the result.
for (int i = 0; i < 4; i++) {
if (result(x_vec[i], y_vec[i]) != val[i]) {
printf("There was an error at %d %d!\n",
x_vec[i], y_vec[i]);
return -1;
}
}
}
{
// y_pairs = 1
int y = y_base + 1;
int y_vec[4] = {y, y, y, y};
int val[4] = {x_vec[0] + y_vec[0],
x_vec[1] + y_vec[1],
x_vec[2] + y_vec[2],
x_vec[3] + y_vec[3]};

// Check the result.
for (int i = 0; i < 4; i++) {
if (result(x_vec[i], y_vec[i]) != val[i]) {
printf("There was an error at %d %d!\n",
x_vec[i], y_vec[i]);
return -1;
}
}
}
}
}
}
printf("\n");

printf("Pseudo-code for the schedule:\n");
gradient_fast.print_loop_nest();
printf("\n");

// Note that in the Halide version, the algorithm is specified
// once at the top, separately from the optimizations, and there
// aren't that many lines of code total. Compare this to the C
// version. There's more code (and it isn't even parallelized or
// vectorized properly). More annoyingly, the statement of the
// algorithm (the result is x plus y) is buried in multiple places
// within the mess. This C code is hard to write, hard to read,
// hard to debug, and hard to optimize further. This is why Halide
// exists.
// 从前面给出的几个调度的例子可以看出，Halide对应版本的代码相对于C的代码，算法和调度分别进行分离
// 调度方便，发ingbianjinxing游动优化，而对应的C语言代码，冗长杂乱，很难写/难读/调试困难，
// 不方便进一步优化
}

printf("Success!\n");
return 0;
}

编译执行：

$ g++ lesson_05*.cpp -g -I ../include -L ../bin -lHalide -lpthread -ldl -o lesson_05 -std=c++11
$ ./lesson_05

执行结果：

由于篇幅太长，略

代码详解与效果图示：

1 算法描述，对于所有版本都是一样，和算法调度/优化高度分离，无依赖性

gradient(x, y) = x + y;

2 reorder函数，默认是x为内循环，y为外循环

gradient.reorder(y, x); // 交换x和y循环的顺序

x为内循环



y为内循环



3 split函数，将一个大的变量拆分成两个小的变量进行循环

for x in(0, 20):
//...

Func.split(x, x_outer, x_inner, 5)；

// 等价于
for x_outer in(0,4):
for x_inner in(0, 5):
// ...

4 fuse函数，将多个循环合并成一个大循环，可以方便进行多个维度的并行调度

for x in(1,4):
for y in(1,4):
//...

funcs.fuse(x, y, xy):
//等价于
for xy in (1,16):
x = xy / 4;
y = xy % 4;
// ...

5 tile将整个图像划分成多个小矩形，分别在小矩形上进行处理，每个tile内部仍然是按照行优先的调度进行计算



6 vertorize向量化计算，调用x86的SSE指令，通常是单指令多数据指令（SIMD），加快计算



7 unroll，将循环铺平，铺平并不会影响计算顺序，当循环内部有数据共享时，铺平可以将重复计算的数据只需读一次和计算一次，减少重复计算

8 原图像尺寸不能整数split后tile小矩形尺寸情形



9 tile的外层循环融合成一个整体循环，在此基础上调用多核并行

Var x_outer, y_outer, x_inner, y_inner, tile_index;
gradient.tile(x, y, x_outer, y_outer, x_inner, y_inner, 4, 4);
gradient.fuse(x_outer, y_outer, tile_index);
gradient.parallel(tile_index);

10 综合调度

// tile尺寸为64x64，采用并行计算
Var x_outer, y_outer, x_inner, y_inner, tile_index;
gradient_fast
.tile(x, y, x_outer, y_outer, x_inner, y_inner, 64, 64)
.fuse(x_outer, y_outer, tile_index)
.parallel(tile_index);

// 将内部的64x64tile继续拆分成更小的tile。在x方向上采用向量化的计算（调用SIMD指令），
// 在y方向进行平铺
Var x_inner_outer, y_inner_outer, x_vectors, y_pairs;
gradient_fast
.tile(x_inner, y_inner, x_inner_outer, y_inner_outer, x_vectors, y_pairs, 4, 2)
.vectorize(x_vectors)
.unroll(y_pairs);