sse - 使用数组进行特征向量化

Question

我正在处理点云数据(每个云 150k 点)。我想，对于每个 (x,y) 点，计算到引用点 O 的距离和方位角:

for each point p in points
    dx = p.x - ox
    dy = p.y - oy
    d = hypot(dx, dy)
    az = atan2(dy, dx)

我有一个手动 SSE 实现。我希望使用 eigen 使代码更清晰:

ArrayXf x(points.size()), y(points.size());
for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}
const ArrayXf d = (dx.square() + dy.square()).sqrt();
// implement a polynomial approximation to atan (same as the SSE)

但是，从我的计时实验来看，这似乎根本没有矢量化，因为时间与基线实现相同。我知道 SSE2 已启用，因为我正在同一文件中编译一些 SSE2 代码。

但是，根据文档，Eigen 在支持时确实利用了 SSE2(以及 3.3 中的 AVX)。它仅适用于向量和矩阵运算吗？

编辑:我研究了生成的汇编代码，它确实包含一些 SSE 指令。但还是很慢

编辑:这里有更多计时信息。我循环了 100 多个帧，每帧大约 150k 点。

• 没有atan2的简单实现:150ms
• sse 实现(处理 4 x 4 点并丢弃最后几个未填满数据包的点):30ms
• 使用特征图实现特征值:90ms(diff:36ms，hypot:16ms，index:17ms)

这是我的特征码:

const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > px(&(points[0].x), points.size());
const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > py(&(points[0].y), points.size());

// difference with the origin (ox and oy are floats)
const Eigen::ArrayXf dx = px - ox, dy = py - oy;

// distance and index
const Eigen::ArrayXf d = sqrt(dx.square() + dy.square());

static const float r_res_mult = 1.0f / r_res; //2x faster than div
const Eigen::ArrayXi didx = (d * r_res_mult).cast<int>();

Answer 1

您的主要问题是您的数据格式不适合 SIMD。您正在使用结构数组(xyxyxyxyxyxy ...)，然后对您所做的代码进行矢量化

for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}

转换为数组结构(xxxxxxxx....yyyyyyy...)。这种转换的成本很高。

更好的解决方案是将您的点存储为数组结构。更好的解决方案是使用数组的混合结构，也称为数组结构的数组。对于 SSE，假设您使用单个浮点，那么您将执行 xxxxyyyyxxxxyyyy....

接下来我建议您使用 SIMD 数学库。 Intel 提供了昂贵且闭源的 SVML。 AMD 优惠libm这是免费但闭源的。但这些库在竞争对手的硬件上都不能很好地运行。最好的 SIMD 库是 Agner Fog 的 Vector Class Library (VCL) 。它是开源、免费的，并且经过优化，可以在 Intel 和 AMD 处理器上运行。它也像 Eigen 一样，只是头文件，因此，像 Eigen 一样，您不必编译和链接库。您刚刚包含了头文件。以下是您如何为 SSE 或 AVX 实现 float (VLC 将在没有 AVX 的系统上模拟 AVX)。

//    g++ -O3 -Ivectorclass -msse4.2 foo.cpp
// or g++ -O3 -Ivectorclass -mavx foo.cpp
#include <vectorclass.h>
#include <vectormath_trig.h>

struct Point2DBlock {
    float x[8];
    float y[8];
};

int main(void) {
    const int nblocks = 10; //each block contains eight points
    Point2DBlock aosoa[nblocks]; //xxxxxxxxyyyyyyyy xxxxxxxxyyyyyyyy ...
    float ox = 0.0f, oy = 0.0f;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<nblocks; i++) {
        Vec8f dx = Vec8f().load(aosoa[i].x) - vox;
        Vec8f dy = Vec8f().load(aosoa[i].y) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}

如果您确实需要hypot。您可以使用 pseudo-code from wikipedia 从 VCL 构建一个.

static inline Vec8f hypot(Vec8f const &x, Vec8f const &y) {
    Vec8f t;
    Vec8f ax = abs(x), ay = abs(y);
    t  = min(ax,ay);
    ax = max(ax,ay);
    t  = t/ax;
    return ax*sqrt(1+t*t);
}

编辑:

这是一个使用结构数组的方法。这需要一些改组，但与其他计算相比，这可能可以忽略不计。 VLC 使用模板元编程来确定有效的洗牌方法。

#include <vectorclass.h>
#include <vectormath_trig.h>

int main(void) {
    const int npoints=80;
    float points[2*npoints]; //xyxyxyxyxyxy...
    float ox = 0.0, oy = 0.0;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<npoints; i+=16) {
        Vec8f l1 = Vec8f().load(&points[i+0]);
        Vec8f l2 = Vec8f().load(&points[i+8]);
        Vec8f dx = blend8f<0, 2, 4, 6, 8, 10, 12, 14>(l1,l2) - vox;
        Vec8f dy = blend8f<1, 3, 5, 7, 9, 11, 13, 15>(l1,l2) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}