当前位置: 首页 > 工具软件 > Halide > 使用案例 >

Halide示例学习五

路昆杰
2023-12-01

函数重定义操作

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    {
        // 这里有着var到expr的转换,是一种原始定义
        Func f;
        f(x, y) = x + y;
        // 这里是单点定义
        f(3, 7) = 42;

        // 递归定义:表示在同一位置递归调用本身值f(x,y) = g(f(x,y), var)
        f(x, y) = f(x, y) + 17;

        // 特指某一行进行更新
        f(x, 3) = f(x, 0) * f(x, 10);

        // 特指某一列进行更新
        f(0, y) = f(0, y) / f(3, y);

        // 更新定义的规则是:左右两边的var参数必须是(x,其他),或者是(其他,y),这跟之前的f定义是一致的。
        // 以下部分是合法的
        f(x, 17) = x + 8; // y >= 18
        f(0, y) = y * 8; // x >= 4
        f(x, x + 1) = x + 8; // y >= x + 1
        f(y / 2, y) = f(0, y) * 17; // x >= y / 2

        // 下面的部分是非法的

        // f(x, 0) = f(x + 1, 0);
        // First argument to f on the right-hand-side must be 'x', not 'x + 1'.

        // f(y, y + 1) = y + 8;
        // Second argument to f on the left-hand-side must be 'y', not 'y + 1'.

        // f(y, x) = y - x;
        // Arguments to f on the left-hand-side are in the wrong places.

        // f(3, 4) = x + y;
        // Free variables appear on the right-hand-side but not the left-hand-side.

        // 这里需要根据前面的重定义函数来评估x,y的尺寸关系
        f.realize({100, 101});

        // 每一次的重定义,都会产生多次更新动作,那么在realize的时候将会体现每次的重定义
        Func g("g");
        g(x, y) = x + y;    // Pure definition
        g(2, 1) = 42;       // First update definition
        g(x, 0) = g(x, 1);  // Second update definition

        g.trace_loads();//跟踪读数据
        g.trace_stores();//跟踪写数据

        g.realize({4, 4});
        // 下述是等效C代码
        int result[4][4];
        // Pure definition
        for (int y = 0; y < 4; y++) {
            for (int x = 0; x < 4; x++) {
                result[y][x] = x + y;
            }
        }
        // First update definition
        result[1][2] = 42;
        // Second update definition
        for (int x = 0; x < 4; x++) {
            result[0][x] = result[1][x];
        }
    }

    return 0;
}
# 结果如下所示
Store g.0(0, 0) = 0
Store g.0(1, 0) = 1
Store g.0(2, 0) = 2
Store g.0(3, 0) = 3
Store g.0(0, 1) = 1
Store g.0(1, 1) = 2
Store g.0(2, 1) = 3
Store g.0(3, 1) = 4
Store g.0(0, 2) = 2
Store g.0(1, 2) = 3
Store g.0(2, 2) = 4
Store g.0(3, 2) = 5
Store g.0(0, 3) = 3
Store g.0(1, 3) = 4
Store g.0(2, 3) = 5
Store g.0(3, 3) = 6
Store g.0(2, 1) = 42
Load g.0(0, 1) = 1
Store g.0(0, 0) = 1
Load g.0(1, 1) = 2
Store g.0(1, 0) = 2
Load g.0(2, 1) = 42
Store g.0(2, 0) = 42
Load g.0(3, 1) = 4
Store g.0(3, 0) = 4

RDom下标使用

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 把更新放置到内部循环中
    {
        // Starting with this pure definition:
        Func f;
        f(x, y) = (x + y) / 100.0f;

        // Say we want an update that squares the first fifty rows. We
        // could do this by adding 50 update definitions:

        // f(x, 0) = f(x, 0) * f(x, 0);
        // f(x, 1) = f(x, 1) * f(x, 1);
        // f(x, 2) = f(x, 2) * f(x, 2);
        // ...
        // f(x, 49) = f(x, 49) * f(x, 49);

        // Or equivalently using a compile-time loop in our C++:
        // for (int i = 0; i < 50; i++) {
        //   f(x, i) = f(x, i) * f(x, i);
        // }

        // 使用RDom变量,我们可以把限定范围的更新组织成循环
        RDom r(0, 50);
        f(x, r) = f(x, r) * f(x, r);
        Buffer<float> halide_result = f.realize({100, 100});

        // See figures/lesson_09_update_rdom.mp4 for a visualization.

        // The equivalent C is:
        float c_result[100][100];
        for (int y = 0; y < 100; y++) {
            for (int x = 0; x < 100; x++) {
                c_result[y][x] = (x + y) / 100.0f;
            }
        }
        for (int x = 0; x < 100; x++) {
            for (int r = 0; r < 50; r++) {
                // RDom这个限定轴,被内嵌为内部循环
                c_result[r][x] = c_result[r][x] * c_result[r][x];
            }
        }

        // Check the results match:
        for (int y = 0; y < 100; y++) {
            for (int x = 0; x < 100; x++) {
                if (fabs(halide_result(x, y) - c_result[y][x]) > 0.01f) {
                    printf("halide_result(%d, %d) = %f instead of %f\n",
                           x, y, halide_result(x, y), c_result[y][x]);
                    return -1;
                }
            }
        }
    }


    return 0;
}

使用RDom计算直方图

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 现在使用重定义来实现直方图
    {
        Func histogram("histogram");
        // 直方图置零
        histogram(x) = 0;
        // 定义一个多维的归化维度
        RDom r(0, input.width(), 0, input.height());
        // 统计input的直方图(0,255)灰度值之间的个数
        histogram(input(r.x, r.y)) += 1;
        // 只保留[0,255]之间的值
        Buffer<int> halide_result = histogram.realize({256});
        // 等效C代码如下所示
        int c_result[256];
        for (int x = 0; x < 256; x++) {
            c_result[x] = 0;
        }
        for (int r_y = 0; r_y < input.height(); r_y++) {
            for (int r_x = 0; r_x < input.width(); r_x++) {
                c_result[input(r_x, r_y)] += 1;
            }
        }

        // Check the answers agree:
        for (int x = 0; x < 256; 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;
            }
        }
    }


    return 0;
}

调度函数更新

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 调度函数更新顺序
    {
        // 对于重定义的函数,我们该如何实现向量化、并行化等操作呢

        // Consider the definition:
        Func f;
        f(x, y) = x * y;
        // Set row zero to each row 8
        f(x, 0) = f(x, 8);
        // Set column zero equal to column 8 plus 2
        f(0, y) = f(8, y) + 2;

        // 原始的 f(x,y) = x * y向量化后再并行化
        f.vectorize(x, 4).parallel(y);
        // 使用Func::update(int upid)来指定第upid次的函数更新的后续优化手段,这里upid=0只能对x那个轴进行向量化
        f.update(0).vectorize(x, 4);
        // upid = 1,表示在第二update公式上进行操作
        Var yo, yi;
        f.update(1).split(y, yo, yi, 4).parallel(yo);//分割后向量化,再接并行

        Buffer<int> halide_result = f.realize({16, 16});

        // See figures/lesson_09_update_schedule.mp4 for a visualization.

        // Here's the equivalent (serial) C:
        int c_result[16][16];

        // Pure step. Vectorized in x and parallelized in y.
        for (int y = 0; y < 16; y++) {  // Should be a parallel for loop
            for (int x_vec = 0; x_vec < 4; x_vec++) {
                int x[] = {x_vec * 4, x_vec * 4 + 1, x_vec * 4 + 2, x_vec * 4 + 3};
                c_result[y][x[0]] = x[0] * y;
                c_result[y][x[1]] = x[1] * y;
                c_result[y][x[2]] = x[2] * y;
                c_result[y][x[3]] = x[3] * y;
            }
        }

        // First update. Vectorized in x.
        for (int x_vec = 0; x_vec < 4; x_vec++) {
            int x[] = {x_vec * 4, x_vec * 4 + 1, x_vec * 4 + 2, x_vec * 4 + 3};
            c_result[0][x[0]] = c_result[8][x[0]];
            c_result[0][x[1]] = c_result[8][x[1]];
            c_result[0][x[2]] = c_result[8][x[2]];
            c_result[0][x[3]] = c_result[8][x[3]];
        }

        // Second update. Parallelized in chunks of size 4 in y.
        for (int yo = 0; yo < 4; yo++) {  // Should be a parallel for loop
            for (int yi = 0; yi < 4; yi++) {
                int y = yo * 4 + yi;
                c_result[y][0] = c_result[y][8] + 2;
            }
        }

        // Check the C and Halide results match:
        for (int y = 0; y < 16; y++) {
            for (int x = 0; x < 16; 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;
                }
            }
        }
    }


    return 0;
}

多阶段重定义操作

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 当涉及到多个计算图的时候,重定义操作又是如何进行的呢
    {
        // 对于默认的多计算内敛模式而言,重定义似乎没有意义
        Func producer, consumer;
        producer(x) = x * 2;//纯定义
        producer(x) += 10;//更新
        consumer(x) = 2 * producer(x);//内敛
        Buffer<int> halide_result = consumer.realize({10});

        // See figures/lesson_09_inline_reduction.gif for a visualization.

        // The equivalent C is:
        int c_result[10];
        for (int x = 0; x < 10; x++) {
            int producer_storage[1];
            // Pure step for producer
            producer_storage[0] = x * 2;
            // Update step for producer
            producer_storage[0] = producer_storage[0] + 10;
            // Pure step for consumer
            c_result[x] = 2 * producer_storage[0];
        }

        // Check the results match
        for (int x = 0; x < 10; x++) {
            if (halide_result(x) != c_result[x]) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    // 对comsumer这个后计算图进行重定义操作
    // 情况1: consumer只引用producer的纯定义
    {
        Func producer, consumer;

        // The producer is pure.
        producer(x) = x * 17;
        consumer(x) = 2 * producer(x);
        consumer(x) += 50;

        // The valid schedules for the producer in this case are
        // the default schedule - inlined, and also:
        //
        // 1) producer.compute_at(x), which places the computation of
        // the producer inside the loop over x in the pure step of the
        // consumer.
        //
        // 2) producer.compute_root(), which computes all of the
        // producer ahead of time.
        //
        // 3) producer.store_root().compute_at(x), which allocates
        // space for the consumer outside the loop over x, but fills
        // it in as needed inside the loop.
        //
        // Let's use option 1.

        producer.compute_at(consumer, x);

        Buffer<int> halide_result = consumer.realize({10});

        // See figures/lesson_09_compute_at_pure.gif for a visualization.

        // The equivalent C is:
        int c_result[10];
        // compute_at使得prudcer先算,又因为内建的原因,consumer的内建操作使得完成consumer的纯定义的计算
        for (int x = 0; x < 10; x++) {
            // Pure step for producer
            int producer_storage[1];
            producer_storage[0] = x * 17;
            c_result[x] = 2 * producer_storage[0];
        }
        // 然后重定义计算
        for (int x = 0; x < 10; x++) {
            c_result[x] += 50;
        }
        // Check the results match
        for (int x = 0; x < 10; x++) {
            if (halide_result(x) != c_result[x]) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    {
        // 情况 2: consumer只在重定义的时候引用producer
        Func producer, consumer;
        producer(x) = x * 17;
        consumer(x) = 100 - x * 10;
        consumer(x) += producer(x);
        producer.compute_at(consumer, x);

        // Note however, that we didn't say:
        //
        // producer.compute_at(consumer.update(0), x).
        //
        // Scheduling is done with respect to Vars of a Func, and
        // the Vars of a Func are shared across the pure and
        // update steps.

        Buffer<int> halide_result = consumer.realize({10});

        // See figures/lesson_09_compute_at_update.gif for a visualization.

        // The equivalent C is:
        int c_result[10];
        // 纯定义
        for (int x = 0; x < 10; x++) {
            c_result[x] = 100 - x * 10;
        }
        // 重定义部分
        for (int x = 0; x < 10; x++) {
            // producer的纯定义
            int producer_storage[1];
            producer_storage[0] = x * 17;
            c_result[x] += producer_storage[0];
        }

        // Check the results match
        for (int x = 0; x < 10; x++) {
            if (halide_result(x) != c_result[x]) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    {
        // 情况 3: 在纯定义和重定义的时候都使用了producer的情况
        Func producer, consumer;
        producer(x) = x * 17;
        consumer(x) = 170 - producer(x);
        consumer(x) += producer(x) / 2;
        producer.compute_at(consumer, x);

        Buffer<int> halide_result = consumer.realize({10});

        // See figures/lesson_09_compute_at_pure_and_update.gif for a visualization.
        // 不管如何,中间计算过程都没有专门为了保存producer而设置的中间Buffer
        int c_result[10];
        // Pure step for the consumer
        for (int x = 0; x < 10; x++) {
            // Pure step for producer
            int producer_storage[1];
            producer_storage[0] = x * 17;
            c_result[x] = 170 - producer_storage[0];
        }
        // Update step for the consumer
        for (int x = 0; x < 10; x++) {
            // Another copy of the pure step for producer
            int producer_storage[1];
            producer_storage[0] = x * 17;
            c_result[x] += producer_storage[0] / 2;
        }

        // Check the results match
        for (int x = 0; x < 10; x++) {
            if (halide_result(x) != c_result[x]) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }

    {
        // 情况 4:在多次重定义的过程中才使用producer
        Func producer, consumer;
        producer(x, y) = (x * y) / 10 + 8;
        consumer(x, y) = x + y;
        consumer(x, 0) += producer(x, x);
        consumer(0, y) += producer(y, 9 - y);

        // In this case neither producer.compute_at(consumer, x)
        // nor producer.compute_at(consumer, y) will work, because
        // either one fails to cover one of the uses of the
        // producer. So we'd have to inline producer, or use
        // producer.compute_root().

        // Let's say we really really want producer to be
        // compute_at the inner loops of both consumer update
        // steps. Halide doesn't allow multiple different
        // schedules for a single Func, but we can work around it
        // by making two wrappers around producer, and scheduling
        // those instead:

        // Attempt 2:
        Func producer_1, producer_2, consumer_2;
        producer_1(x, y) = producer(x, y);
        producer_2(x, y) = producer(x, y);

        consumer_2(x, y) = x + y;
        consumer_2(x, 0) += producer_1(x, x);
        consumer_2(0, y) += producer_2(y, 9 - y);

        // The wrapper functions give us two separate handles on
        // the producer, so we can schedule them differently.
        producer_1.compute_at(consumer_2, x);
        producer_2.compute_at(consumer_2, y);

        Buffer<int> halide_result = consumer_2.realize({10, 10});

        // See figures/lesson_09_compute_at_multiple_updates.mp4 for a visualization.

        // The equivalent C is:
        int c_result[10][10];
        // Pure step for the consumer
        for (int y = 0; y < 10; y++) {
            for (int x = 0; x < 10; x++) {
                c_result[y][x] = x + y;
            }
        }
        // First update step for consumer
        for (int x = 0; x < 10; x++) {
            int producer_1_storage[1];
            producer_1_storage[0] = (x * x) / 10 + 8;
            c_result[0][x] += producer_1_storage[0];
        }
        // Second update step for consumer
        for (int y = 0; y < 10; y++) {
            int producer_2_storage[1];
            producer_2_storage[0] = (y * (9 - y)) / 10 + 8;
            c_result[y][0] += producer_2_storage[0];
        }

        // Check the results match
        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;
                }
            }
        }
    }

    {
        // Case 5: 当使用了RDom限定下标的时候,使用compute_at

        // We are not just restricted to scheduling producers at
        // the loops over the pure variables of the consumer. If a
        // producer is only used within a loop over a reduction
        // domain (RDom) variable, we can also schedule the
        // producer there.

        Func producer, consumer;

        RDom r(0, 5);
        producer(x) = x % 8;
        consumer(x) = x + 10;
        consumer(x) += r + producer(x + r);

        producer.compute_at(consumer, r);

        Buffer<int> halide_result = consumer.realize({10});

        // See figures/lesson_09_compute_at_rvar.gif for a visualization.

        // The equivalent C is:
        int c_result[10];
        // Pure step for the consumer.
        for (int x = 0; x < 10; x++) {
            c_result[x] = x + 10;
        }
        // Update step for the consumer.
        for (int x = 0; x < 10; x++) {
            // The loop over the reduction domain is always the inner loop.
            for (int r = 0; r < 5; r++) {
                // We've schedule the storage and computation of
                // the producer here. We just need a single value.
                int producer_storage[1];
                // Pure step of the producer.
                producer_storage[0] = (x + r) % 8;

                // Now use it in the update step of the consumer.
                c_result[x] += r + producer_storage[0];
            }
        }

        // Check the results match
        for (int x = 0; x < 10; x++) {
            if (halide_result(x) != c_result[x]) {
                printf("halide_result(%d) = %d instead of %d\n",
                       x, halide_result(x), c_result[x]);
                return -1;
            }
        }
    }


    return 0;
}

使用归化函数进行模糊操作

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 下面是一个模糊的归化函数
    {
        // 下述是一个边界限定的函数
        Func clamped = BoundaryConditions::repeat_edge(input);
        // 定义一个从该点开始的相对位置为-2的点
        RDom r(-2, 5, -2, 5);
        // 计算每个像素值的附近的和
        Func local_sum;
        local_sum(x, y) = 0;// 纯定义
        local_sum(x, y) += clamped(x + r.x, y + r.y); // 归约求和
        // 求平均值,然后int32->uint8
        Func blurry;
        blurry(x, y) = cast<uint8_t>(local_sum(x, y) / 25);
        Buffer<uint8_t> halide_result = blurry.realize({input.width(), input.height()});

        // 默认的调度就是把clamped进行了内联
        // 等效C代码如下
        Buffer<uint8_t> c_result(input.width(), input.height());
        for (int y = 0; y < input.height(); y++) {
            for (int x = 0; x < input.width(); x++) {
                int local_sum[1];
                // Pure step of local_sum
                local_sum[0] = 0;
                // Update step of local_sum
                for (int r_y = -2; r_y <= 2; r_y++) {
                    for (int r_x = -2; r_x <= 2; r_x++) {
                        // The clamping has been inlined into the update step.
                        int clamped_x = std::min(std::max(x + r_x, 0), input.width() - 1);
                        int clamped_y = std::min(std::max(y + r_y, 0), input.height() - 1);
                        local_sum[0] += input(clamped_x, clamped_y);
                    }
                }
                // Pure step of blurry
                c_result(x, y) = (uint8_t)(local_sum[0] / 25);
            }
        }

        // Check the results match
        for (int y = 0; y < input.height(); y++) {
            for (int x = 0; x < input.width(); x++) {
                if (halide_result(x, y) != c_result(x, y)) {
                    printf("halide_result(%d, %d) = %d instead of %d\n",
                           x, y, halide_result(x, y), c_result(x, y));
                    return -1;
                }
            }
        }
    }


    return 0;
}

内置的归化函数sum

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 归约函数sum
    {
        // Halide.h中包含了sum归约函数
        Func f1;
        RDom r(0, 100);
        f1(x) = sum(r + x) * 7;
        // 上述操作等下等效
        Func f2;
        Func anon;
        anon(x) = 0;
        anon(x) += r + x;//固定x,那么对于每一个r[0,100],我们必须计算anon[x] = [x + r for r in [0,100]]
        f2(x) = anon(x) * 7;

        // 所以尽管f1使用了sum函数,但是这个依然是个纯定义函数
        Buffer<int> halide_result_1 = f1.realize({10});
        Buffer<int> halide_result_2 = f2.realize({10});

        // The equivalent C is:
        int c_result[10];
        for (int x = 0; x < 10; x++) {
            int anon[1];
            anon[0] = 0;
            for (int r = 0; r < 100; r++) {
                anon[0] += r + x;
            }
            c_result[x] = anon[0] * 7;
        }

        // Check they all match.
        for (int x = 0; x < 10; x++) {
            if (halide_result_1(x) != c_result[x]) {
                printf("halide_result_1(%d) = %d instead of %d\n",
                       x, halide_result_1(x), c_result[x]);
                return -1;
            }
            if (halide_result_2(x) != c_result[x]) {
                printf("halide_result_2(%d) = %d instead of %d\n",
                       x, halide_result_2(x), c_result[x]);
                return -1;
            }
        }
    }


    return 0;
}

使用归约函数的一个复杂例子

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

using namespace Halide;

int main(int argc, char **argv) {
    Var x("x"), y("y");
    Buffer<uint8_t> input = load_image("../../images/gray.png");
    // 一个使用归约函数的复杂例子
    {
        // 其他的归约函数如下所示:
        // 累乘"product", 最小值"minimum",
        // 最大值"maximum", 反解最小"argmin", 反解最大"argmax". 

        // 首先进行边界限定
        Func clamped;
        Expr x_clamped = clamp(x, 0, input.width() - 1);
        Expr y_clamped = clamp(y, 0, input.height() - 1);
        clamped(x, y) = input(x_clamped, y_clamped);

        RDom box(-2, 5, -2, 5);
        // 然后进行如下归约
        Func spread;
        spread(x, y) = (maximum(clamped(x + box.x, y + box.y)) -
                        minimum(clamped(x + box.x, y + box.y)));

        // 循环拆分、并行化
        Var yo, yi;
        spread.split(y, yo, yi, 32).parallel(yo);

        // 使用向量化,maximum和minimum归约函数也会受到向量化影响
        spread.vectorize(x, 16);

        // 为了放置边界被多次计算,这里使用了store_at和compute_at
        clamped.store_at(spread, yo).compute_at(spread, yi);

        Buffer<uint8_t> halide_result = spread.realize({input.width(), input.height()});

// 这里使用优化的C代码和Halide代码进行比较,来看两者的速度差别
#define __SSE2__
#define _OPENMP

#ifdef __SSE2__
        // Don't include the time required to allocate the output buffer.
        Buffer<uint8_t> c_result(input.width(), input.height());

#ifdef _OPENMP
        double t1 = current_time();
#endif

        // Run this one hundred times so we can average the timing results.
        for (int iters = 0; iters < 100; iters++) {

#pragma omp parallel for
            for (int yo = 0; yo < (input.height() + 31) / 32; yo++) {
                int y_base = std::min(yo * 32, input.height() - 32);

                // Compute clamped in a circular buffer of size 8
                // (smallest power of two greater than 5). Each thread
                // needs its own allocation, so it must occur here.

                int clamped_width = input.width() + 4;
                uint8_t *clamped_storage = (uint8_t *)malloc(clamped_width * 8);//中间计算结果,store_at

                for (int yi = 0; yi < 32; yi++) {//循环拆分的内层循环
                    int y = y_base + yi;

                    uint8_t *output_row = &c_result(0, y);

                    // compute_at的效用,使得在向量化前计算了clamped
                    int min_y_clamped = (yi == 0) ? (y - 2) : (y + 2);
                    int max_y_clamped = (y + 2);
                    for (int cy = min_y_clamped; cy <= max_y_clamped; cy++) { //compute_at,用于进行计算操作
                        // Figure out which row of the circular buffer
                        // we're filling in using bitmasking:
                        uint8_t *clamped_row =
                            clamped_storage + (cy & 7) * clamped_width;

                        // Figure out which row of the input we're reading
                        // from by clamping the y coordinate:
                        int clamped_y = std::min(std::max(cy, 0), input.height() - 1);
                        uint8_t *input_row = &input(0, clamped_y);

                        // Fill it in with the padding.
                        for (int x = -2; x < input.width() + 2; x++) {
                            int clamped_x = std::min(std::max(x, 0), input.width() - 1);
                            *clamped_row++ = input_row[clamped_x];
                        }
                    }

                    // 现在进行最终的计算
                    for (int x_vec = 0; x_vec < (input.width() + 15) / 16; x_vec++) {
                        int x_base = std::min(x_vec * 16, input.width() - 16);

                        // Allocate storage for the minimum and maximum
                        // helpers. One vector is enough.
                        __m128i minimum_storage, maximum_storage;

                        // The pure step for the maximum is a vector of zeros
                        maximum_storage = _mm_setzero_si128();

                        // maximum操作
                        for (int max_y = y - 2; max_y <= y + 2; max_y++) {
                            uint8_t *clamped_row =
                                clamped_storage + (max_y & 7) * clamped_width;
                            for (int max_x = x_base - 2; max_x <= x_base + 2; max_x++) {
                                __m128i v = _mm_loadu_si128(
                                    (__m128i const *)(clamped_row + max_x + 2));
                                maximum_storage = _mm_max_epu8(maximum_storage, v);
                            }
                        }

                        // 定义一个0值的minimum
                        minimum_storage = _mm_cmpeq_epi32(_mm_setzero_si128(),
                                                          _mm_setzero_si128());
                        // minmum操作
                        for (int min_y = y - 2; min_y <= y + 2; min_y++) {
                            uint8_t *clamped_row =
                                clamped_storage + (min_y & 7) * clamped_width;
                            for (int min_x = x_base - 2; min_x <= x_base + 2; min_x++) {
                                __m128i v = _mm_loadu_si128(
                                    (__m128i const *)(clamped_row + min_x + 2));
                                minimum_storage = _mm_min_epu8(minimum_storage, v);
                            }
                        }

                        // sse减法
                        __m128i spread = _mm_sub_epi8(maximum_storage, minimum_storage);

                        // sse存储
                        _mm_storeu_si128((__m128i *)(output_row + x_base), spread);
                    }
                }

                free(clamped_storage);
            }
        }

// Skip the timing comparison if we don't have openmp
// enabled. Otherwise it's unfair to C.
#ifdef _OPENMP
        double t2 = current_time();

        // Now run the Halide version again without the
        // jit-compilation overhead. Also run it one hundred times.
        for (int iters = 0; iters < 100; iters++) {
            spread.realize(halide_result);
        }

        double t3 = current_time();
        printf("Halide spread took %f ms. C equivalent took %f ms\n",
               (t3 - t2) / 100, (t2 - t1) / 100);

#endif  // _OPENMP

        // Check the results match:
        for (int y = 0; y < input.height(); y++) {
            for (int x = 0; x < input.width(); x++) {
                if (halide_result(x, y) != c_result(x, y)) {
                    printf("halide_result(%d, %d) = %d instead of %d\n",
                           x, y, halide_result(x, y), c_result(x, y));
                    return -1;
                }
            }
        }

#endif  // __SSE2__
    }


    return 0;
}
# 结果如下所示
Halide spread took 0.668330 ms. C equivalent took 68.161500 ms
 类似资料: