用于通用 SIMD（SSE、AVX、NEON）测试零匹配的高效 C 向量。 （求FP最大绝对值和指数）

我想看看是否可以编写一些可以高效编译的通用 SIMD 代码。主要用于 SSE、AVX 和 NEON。该问题的简化版本是：找到浮点数数组的最大绝对值并返回该值和索引。导致问题的是最后一部分，即最大值的索引。似乎没有一种很好的方法来编写具有分支的代码。

请参阅最后的更新，以获取使用一些建议答案的完成代码。

这是一个示例实现（更完整的版本godbolt https://godbolt.org/#z:OYLghAFBqd5QCxAYwPYBMCmBRdBLAF1QCcAaPECAMzwBtMA7AQwFtMQByARg9KtQYEAysib0QXACx8BBAKoBnTAAUAHpwAMvAFYTStJg1DIApACYAQuYukl9ZATwDKjdAGFUtAK4sGISRqkrgAyeAyYAHI%2BAEaYxCAArAmkAA6oCoRODB7evv6BaRmOAqHhUSyx8Um2mPbFDEIETMQEOT5%2BATV1WY3NBKWRMXGJyQpNLW15nWN9A%2BWVIwCUtqhexMjsHOYAzGHI3lgA1CbbbmP4ggB0CCfYJhoAgvdPZttYNOGHAGrB2BGHAA5DgB6YGHeZxQ6oKiHWqYNiCBSHMKHABumAcJFIgKhxEOkmeoLRGKIeOhhyotFQTAICmeBAAnilMO8KVSaWiIkiAPrcmkEYh4aJeAiYXlQdGY4jcjIAL0wEB%2BfwAVHLMNDqOyCIsdYsTlZHkTJaSoTDzGZoqhPJhDHSzGZ6UyWZgYWECJzoodefzBcLReKIMaSDK8PLFb8IqrQ%2BqqBA3bq9dsDU9DWCPAwZojDgQEByAO4kADWFJIh0tOZxhnQ%2BNh9ARtOeOzwMKV/xOABEO4DG693mFMN8IExSNFSMhSOggnxSMBSAhFodh6Px5PMNPZzdHjtakoe28Xf3B0uxxOp1QZ3OF8eV3vXM2933PkIAOoPZRcMxAr4QaTJABspAAOzYoEXCkGYpDbHqW69geT6vsof6SIOEHbCBpBgYchwAcB%2BKkAk0EvPuHwDi%2Bb7bGYg6YYEWFoYcEFYckeFYbhf6Ec8MyOMgZZWrQi4MAygYRJ6qgLiYgHJlhWHNoc4Z/Mc2x3Ns7YAoRUnqYcqgKe2mnHIBbhetywp0I4DAyggXhUJSmBBsQECqKQDmHGR76fomklSVpXZefphnGbQpnmZZ1m2fZjnYi5SHuc86lecpuniQZvL%2BYFCgWVZ9ChQ5TkuRR0WPOpxCYAQawMJpJgJBYGgVZ2SaNoBnZbqmhwAJIMFgWnkjQxBjNmxBeAOqJiANyJlaolyEmCAAqCB4EiKTEBgXgbEi6LEAyZZMNWaBYNiCioIceYDtoXi9egB2EFCDAbIu7o5gOrimtmCADlSqApBNzXJV4JlhNyDiysclheiwLBmAkf7ciwqDoiwTAKIW3IpEic2HCwy0IGWxWimSZUPF8AAan0PJxeDcW6h2zcgCBCSJYkSTFUn8HZXgZngwDhNWeDaYc1VJsiCkGa2%2BrItY1gLjJ9kVRYeA1QuRUlcQZWy3VBVSUSrOPiyyIwgwB1wvWSLNAOAoDfVjUpg8lLUu6tloAwqKajbyqGT6QoimK3JQGI7OcxAeU6oc6A0sOhysxkHPaww2Lh772vKmEHX0x5WGolyaNMHFOniRY36BPnpAF0XeoNfqjOp%2BnBBw6oZdq1hseR1z1c83zKclniEAN5zAtdhoCkWD3pyHAwIvc2L8Wtsn5caWnSLul2ypCQoyqLBAwdNEDA94PlDwaVJ8/xQAtAfBkHwAYtmhwgIcx%2B17ve8eocyAsCkPMH9gGc16r9975LhiCc/FIgcc7TwftmZuXYCCjwgfFFWbcwGS1bApJSOlVJ6XgWAqSs9swKDzEwFIH4gRdm%2Br9MyaVgqZRJCQCABBSC0OcghQhO9MEPwPjpE%2BOC8EEM/ApC%2BtIuGEKvtmO%2BLCsLiQtqIzkc9cH4KQjzEhAU/rkIyjZKhdlaH0MipIZhoi2GXxOKfGRKQ5EnD4UYuR18oHf0kdg/h%2BCKLyKMj9RRZD0ohTUTQuhEUEJ5REaImB7DBacPsZRUxwSUgOMsX4ve4jp6xLrocBO7VMBZ2RBA4GeYqY0yrpnbSXZq46KwgrUqn9pbVQSLVZM8SHgcGWLQTgCReB%2BA4FoUgqBOBuDFsDfaawbo7B4HQzQtTliFn8FwS42xAJmA0GYAAnBoJCGgASAUkMBepHBJBNKGW0zgvAFAgECAQIZyw4CwCQGgF%2BdA4jkEoBclIVz4j7EMMAbkpsGCFj4CZOI%2ByIDRE4AM6IYRmgMn%2BbwC59YADyDBaAgpabwLAcMjDiDhROPARUHB4HRPslFKSMQe1BeQQQtRtm0CFMQYFHgsAEoFHgFgxy%2BAGGAAoL4eBMB5ghcyZpAz%2BCCBEGIdgUgZCCEUCodQKLdBgQMEYFAXSbCkuiPsyAyx3r1GxYfCF2wb5w3WDcZSsoHZxE1YfKy8MCCHzhhWc1etD5MFRKoPZcIMXOAgK4SYfgwIhH7BCeIYFCiZAEG6vQfr6hzCGD6roJIejjFaJ4doeg7CRoEL0FooaKjDDAjMCYsa8gZujamhYXBlg9PWAKuhRVNg8DqQ0rZKL2kcFUACP8h85FPKMLJN5hYFwQE6ZYaw2JcCEFLP07EHhLn0DxP0xYvAjlwp1KQUZn5LgJG2AkKZqzAIJEkH%2BDQq79CcE2aQZprS617IOYM2dpBTmIBQKgMd1yKAQDuQ8kAzQWCogBIfVtwBj79XeZ8gK3zKB/JRYC5g60CXgsYAQKFMLtkIueci1p%2BB0WOCxds3FyB8XcF4G6YlKL5XkvWpSit07BR0ovZSJgTKWVso5YwAlPLhCiHEIKxjIq1DbN0BBKVxhZX6CFIqiAyqUiqs4Oqo12rqYdlfe%2By4TamA3xNWMc1NIEAOu6M6112b3VBHavm9NqR0j%2BuyNpoNRmQ1erDfGx19Rk0xtyDphNTqGh5ss2m8Nmb7NxtzbMNzBai2rBLXoAUmASNVo4I0w92y60yY/Qpr97bf2dtkj2qwlh%2B34BNMOw4o77njqBnladxyRkgG2BoS4GhZn2i3VwbdkgEhcFmWs/dNbj27NsGemdWg53rLMK13gJ7z3deWGtDIzhJBAA%3D):

#define VLEN 8
typedef float vNs __attribute__((vector_size(VLEN*sizeof(float))));
typedef int vNb __attribute__((vector_size(VLEN*sizeof(int))));
#define SWAP128 4,5,6,7, 0,1,2,3
#define SWAP64 2,3, 0,1,  6,7, 4,5
#define SWAP32 1, 0,  3, 2,  5, 4,  7, 6

static bool any(vNb x) {
    x = x | __builtin_shufflevector(x,x, SWAP128);
    x = x | __builtin_shufflevector(x,x, SWAP64);
    x = x | __builtin_shufflevector(x,x, SWAP32);
    return x[0];
}

float maxabs(float* __attribute__((aligned(32))) data, unsigned n, unsigned *index) {
    vNs max = {0,0,0,0,0,0,0,0};
    vNs tmax;
    unsigned imax = 0;
    for (unsigned i = 0 ; i < n; i += VLEN) {
        vNs t = *(vNs*)(data + i);
        t = -t < t ? t : -t;  // Absolute value
        vNb cmp = t > max;
        if (any(cmp)) {
            tmax = t; imax = i;
            // broadcast horizontal max of t into every element of max
            vNs tswap128 = __builtin_shufflevector(t,t, SWAP128);
            t = t < tswap128 ? tswap128 : t;
            vNs tswap64 = __builtin_shufflevector(t,t, SWAP64);
            t = t < tswap64 ? tswap64 : t;
            vNs tswap32 = __builtin_shufflevector(t,t, SWAP32);
            max = t < tswap32 ? tswap32 : t;
        }
    }
    // To simplify example, ignore finding index of true value in tmax==max
    *index = imax; // + which(tmax == max);
    return max[0];
}

godbolt 上的代码允许将 VLEN 更改为 8 或 4。

这大多效果很好。对于 AVX/SSE，绝对值变为t & 0x7fffffff用一个(v)andps，即清除符号位。对于 NEON 来说是这样完成的vneg + fmaxnm。查找和广播水平最大值的块变成了排列和最大值指令的有效序列。 gcc 能够使用 NEONfabs为绝对值。

4 元素 SSE/NEON 目标上的 8 元素向量在 clang 上运行良好。它在两组寄存器上使用一对指令，对于 SWAP128 水平操作将max or or两个寄存器没有任何不必要的排列。另一方面，gcc 确实无法处理这个问题，并且生成的大部分是非 SIMD 代码。如果我们将向量长度减少到 4，gcc 对于 SSE 和 NEON 就可以正常工作。

但有一个问题if (any(cmp))。对于 clang + SSE/AVX，效果很好，vcmpltps + vptest， 与orps在 SSE 上从 8 到 4。

但是 NEON 上的 gcc 和 clang 会执行所有排列和 OR 运算，然后将结果移至 gp 寄存器进行测试。

除了架构特定的内在函数之外，是否还有一些代码可以获取ptest与海湾合作委员会和vmaxvq使用 clang/gcc 和 NEON？

我尝试了一些其他方法，比如if (x[0] || x[1] || ... x[7])但他们更糟。

Update

我创建了一个更新的例子 https://godbolt.org/z/eEh6G4YW9显示了两种不同的实现，即原始方法和 chtz 建议的“向量中的索引”方法以及 Aki Suihkonen 的答案中所示的方法。我们可以看到生成的 SSE 和 NEON 输出。

尽管有些人可能会持怀疑态度，但编译器确实从通用 SIMD（不是自动矢量化！）C++ 代码中生成了非常好的代码。在 SSE/AVX 上，我发现循环中的代码改进空间很小。 NEON 版本仍然受到“any()”的次优实现的困扰。

除非数据通常按升序排列，或几乎如此，否则我的原始版本在 SSE/AVX 上仍然是最快的。我还没有在 NEON 上测试过。这是因为大多数循环迭代都找不到新的最大值，最好针对这种情况进行优化。 “向量中的索引”方法产生更紧密的循环，编译器也做得更好，但常见情况只是在 SSE/AVX 上慢一点。常见情况在 NEON 上可能相同或更快。

关于编写通用 SIMD 代码的一些注意事项。

浮点向量的绝对值可以通过以下方式找到。它在 SSE/AVX（并带有清除符号位的掩码）和 NEON（fabs 指令）上生成最佳代码。

static vNs vabs(vNs x) {
    return -x < x ? x : -x;
}

这将在 SSE/AVX/NEON 上有效地实现垂直最大。它不进行比较；它产生架构的“max”指令。在 NEON 上，将其更改为使用>代替<导致编译器生成非常糟糕的标量代码。我猜是有非规范或异常的东西。

template <typename v>  // Deduce vector type (float, unsigned, etc.)
static v vmax(v a, v b) {
    return a < b ? b : a; // compiles best with "<" as compare op
}

此代码将在寄存器中广播水平最大值。它在 SSE/AVX 上编译得很好。在 NEON 上，如果编译器可以使用水平最大指令然后广播结果，可能会更好。令我印象深刻的是，如果在 SSE/NEON 上使用 8 个元素向量（只有 4 个元素寄存器），编译器就会足够聪明，只使用一个寄存器来广播结果，因为顶部 4 个元素和底部 4 个元素是相同的。

template <typename v>
static v hmax(v x) {
    if (VLEN >= 8)
        x = vmax(x, __builtin_shufflevector(x,x, SWAP128));
    x = vmax(x, __builtin_shufflevector(x,x, SWAP64));
    return vmax(x, __builtin_shufflevector(x,x, SWAP32));
}

这是我发现的最好的“any()”。它在 SSE/AVX 上是最佳的，使用单个 ptest 指令。在 NEON 上，它执行排列和 OR，而不是水平最大指令，但我还没有找到一种方法可以在 NEON 上获得更好的效果。

static bool any(vNb x) {
    if (VLEN >= 8)
        x |= __builtin_shufflevector(x,x, SWAP128);
    x |= __builtin_shufflevector(x,x, SWAP64);
    x |= __builtin_shufflevector(x,x, SWAP32);
    return x[0];
}

同样有趣的是，在 AVX 上的代码i = i + 1将被编译为vpsubd ymmI, ymmI, ymmNegativeOne，即减去-1。为什么？因为 -1s 的向量是由vpcmpeqd ymm0, ymm0, ymm0这比广播 1 向量更快。

这是最好的which()我已经想出来了。这为您提供了布尔向量中第一个真值的索引（0 = false，-1 = true）。使用 movemask 在 AVX 上可以做得更好一些。我不知道最好的NEON。

// vector of signed ints
typedef int vNi __attribute__((vector_size(VLEN*sizeof(int))));
// vector of bytes, same number of elements, 1/4 the size
typedef unsigned char vNb __attribute__((vector_size(VLEN*sizeof(unsigned char))));
// scalar type the same size as the byte vector
using sNb = std::conditional_t<VLEN == 4, uint32_t, uint64_t>;
static int which(vNi x) {
    vNb cidx = __builtin_convertvector(x, vNb);
    return __builtin_ctzll((sNb)cidx) / 8u;
}

正如 chtz 所评论的，最通用和典型的方法是使用另一个掩码来收集索引：

Vec8s indices = { 0,1,2,3,4,5,6,7};
Vec8s max_idx = indices;
Vec8f max_abs = abs(load8(ptr)); 

for (auto i = 8; i + 8 <= vec_length; i+=8) { 
    Vec8s data = abs(load8(ptr[i]));
    auto mask = is_greater(data, max_abs);
    max_idx = bitselect(mask, indices, max_idx);
    max_abs = max(max_abs, data);    
    indices = indices + 8;
}

另一种选择是交错值和索引：

auto data = load8s(ptr) & 0x7fffffff; // can load data as int32_t
auto idx = vec8s{0,1,2,3,4,5,6,7};
auto lo = zip_lo(idx, data);
auto hi = zip_hi(idx, data);

for (int i = 8; i + 8 <= size; i+=8) {
    idx = idx + 8;
    auto d1 = load8s(ptr + i) & 0x7fffffff;
    auto lo1 = zip_lo(idx, d1);
    auto hi1 = zip_hi(idx, d1);
    lo = max_u64(lo, lo1);
    hi = max_u64(hi, hi1);
}

如果输入范围足够小，可以将输入左移，同时将索引中的一些位附加到同一字的 LSB 位，则此方法尤其有利可图。

即使在这种情况下，我们也可以重新利用浮点数中的 1 位，从而节省一半的位/索引选择操作。

auto data0 = load8u(ptr) << 1; // take abs by shifting left 
auto data1 = (load8u(ptr + 8) << 1) + 1;  // encode odd index to data
auto mx = max_u32(data0, data1);  // the LSB contains one bit of index

貌似可以用double作为存储，因为甚至 SSE2 也支持_mm_max_pd（需要注意 Inf/Nan 处理，当重新解释为 64 位双精度数的高位部分时，它们不再编码为 Inf/Nan）。