本文是halide编程指南的连载，已同步至公众号
第 17章: 非矩形域上的约化
// 本课程演示如何使用谓词定义迭代简化域子集的更新（难懂）.
// linux:
// g++ lesson_17*.cpp -g -I <path/to/Halide.h> -L <path/to/libHalide.so> -lHalide -lpthread -ldl -o lesson_17 -std=c++11
// LD_LIBRARY_PATH=<path/to/libHalide.so> ./lesson_17
// On os x:
// g++ lesson_17*.cpp -g -I <path/to/Halide.h> -L <path/to/libHalide.so> -lHalide -o lesson_17 -std=c++11
// DYLD_LIBRARY_PATH=<path/to/libHalide.dylib> ./lesson_17
// If you have the entire Halide source tree, you can also build it by// running:// make tutorial_lesson_17_predicated_rdom/
#include "Halide.h"
#include <stdio.h>
using namespace Halide;
int main(int argc, char **argv) {
// 在第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);
// 边界框不必是最小的。事实上，这个框可以是任意大小的，只要它覆盖了我们想要更新的区域。但是，边界框越紧，生成的循环边界就越紧。Halide会在可能的情况下自动收紧循环边界，但一般来说，最好定义一个最小的边界框.
// 然后，我们使用RDom：：where在该边界框上定义谓词，这样仅当给定谓词的计算结果为true时（即在循环区域内）才执行更新.
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);
// 请看下面的一个可视化做了什么.
图1
// 等效的C代码:
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++) {
// 仅当谓词的计算结果为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);
// 请看下面的一个可视化.
图2
// 等效的C代码:
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;
}
}
}
}
{
// 谓词不仅限于RDom的变量（r.x，r.y，…）。它还可以引用更新定义中的自由变量，甚至可以调用其他Func，或者对同一Func进行递归调用。例如:
Func f("f"), g("g");
Var x("x"), y("y");
f(x, y) = 2 * x + y;
g(x, y) = x + y;
// 此RDom的谓词依赖于“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);
// 可视化图
图3
// f的等效c代码:
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代码:
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;
}
图1
图2
图3
第18章：使用rfactor分解关联约化
// 本课程演示如何使用调度指令“rfactor”并行化或向量化关联约简.
// linux 这样运行:
// g++ lesson_18*.cpp -g -I <path/to/Halide.h> -L <path/to/libHalide.so> -lHalide -lpthread -ldl -o lesson_18 -std=c++11
// LD_LIBRARY_PATH=<path/to/libHalide.so> ./lesson_18
// os x:
// g++ lesson_18*.cpp -g -I <path/to/Halide.h> -L <path/to/libHalide.so> -lHalide -o lesson_18 -std=c++11
// DYLD_LIBRARY_PATH=<path/to/libHalide.dylib> ./lesson_18
// 在原码树:
// make tutorial_lesson_18_parallel_associative_reductions.
#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<uint8_t> input(8, 8, "input");
for (int y = 0; y < 8; ++y) {
for (int x = 0; x < 8; ++x) {
input(x, y) = (rand() % 256);
}
}
{
// 正如前面在第9课中提到的，并行化作为归约域一部分的变量是很棘手的，因为这些变量之间可能存在数据依赖关系.
// 考虑第9课中的直方图示例:
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);
// 可视化图.
图1
// 我们可以对直方图存储桶的初始化进行矢量化，但是由于更新定义中的r.x和r.y之间存在数据依赖关系（即更新是指在上一次迭代中计算的值），因此如果不引入竞争条件，我们就不能并行化或矢量化r.x或r.y。以下代码将产生错误:
// histogram.update().parallel(r.y);
}
{
//但是，请注意直方图操作（这是一种和归约）是关联的。加速关联约化的一个常见技巧是将约化域分割成更小的切片，计算每个切片上的部分结果，然后合并结果。由于每个切片的计算是独立的，我们可以在切片上并行化.
// 回到直方图示例，我们通过定义一个独立计算每一行直方图的中间函数，将简化域分割成行:
Func intermediate("intm_par_manual");
intermediate(i, y) = 0;
RDom rx(0, input.width());
intermediate(input(rx, y) / 32, y) += 1;
// 然后我们定义第二个阶段，对这些部分结果求和:
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);
// 看一个可视化的图
图2.
}
{
// 这种关联约简的手动因式分解可能很繁琐，而且容易出现错误。尽管手动处理直方图相当容易，但它会很快变得复杂，特别是当RDom可能有一个谓词（RDom：：where）时，或者当函数缩减为多维元组时.
// 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将<RVar，Var>对的列表作为输入，其中包含要“并行化”的缩减变量（RVar）列表。在生成的intermediate Func中，所有对该归约变量的引用都替换为对“纯”变量（Vars）的引用。由于通过构造，VAR是无竞争条件的，所以中间缩减现在可以在这些维度上并行化。不在列表中的所有缩减变量都将从原始函数中删除，并“提升”到中间函数.
// 为了生成与手动分解版本相同的代码，我们执行以下操作:
Func intermediate = histogram.update().rfactor({{r.y, y}});
// 我们将{r.y，y}作为参数传递给rfactor，以使直方图在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<int> halide_result = histogram.realize(8);
// 请参见下面的可视化内容.
图3
// 等效c:
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;
}
}
/* 并行*/
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考虑关联约简，我们可以单独使用调度探索各种因子分解策略。给定相同的序列直方图代码:
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);
// 现在我们正在计算每列的独立直方图，我们可以对列进行矢量化.
intermediate.compute_root().update().vectorize(u, 8);
// 请注意，由于对内部维度进行矢量化会改变值添加到最终直方图桶计算的顺序，因此此技巧仅在关联约简是关联*and*交换的情况下有效。rfactor将尝试证明这些属性成立，如果不能成立，它将抛出一个错误.
// 向量化初始化.
intermediate.vectorize(x, 8);
histogram.vectorize(x, 8);
Buffer<int> halide_result = histogram.realize(8);
// 请参见下面的可视化内容.
图4
// 等效c:
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++) {
/* 向量化*/
for (int u_i = 0; u_i < 8; u_i++) {
c_intm[u * 8 + 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;
}
}
}
{
//我们还可以一次在多个维度上切分一个缩减域。这一次，我们将计算域的分块上的部分直方图.
Func histogram("hist_rfactor_tile");
histogram(x) = 0;
RDom r(0, input.width(), 0, input.height());
histogram(input(r.x, r.y) / 32) += 1;
// 我们先把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来生成一个中间函数，独立地计算每个图块的直方图.
Func intermediate = histogram.update().rfactor({{rx_outer, u}, {ry_outer, v}});
// 我们现在可以在瓦片上平行化中间层.
intermediate.compute_root().update().parallel(u).parallel(v);
// 我们还将最外层的平铺索引重新排序，以提供经典的平铺遍历.
intermediate.update().reorder(rx_inner, ry_inner, u, v);
// 向量化初始化.
intermediate.vectorize(x, 8);
histogram.vectorize(x, 8);
Buffer<int> halide_result = histogram.realize(8);
// 请参见下面的可视化内容.
图5
// 等效c:
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;
}
}
}
/* 并行*/
for (int v = 0; v < input.height() / 2; v++) {
/* 并行*/
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;
}
图1
图2
图3
图4
163
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
