有效地移位或大位向量

我有一個很大的記憶體陣列作為一些指標uint64_t * arr（加上大小），它代表普通位。我需要非常有效地（最高性能/最快）將這些位向右移動一些從 0 到 63 的量。

通過移動整個陣列，我的意思是不移動每個元素（如a[i] <<= Shift），而是將其作為單個大位向量進行移動。換句話說，對于每個中間位置i（第一個和最后一個元素除外），我可以在回圈中執行以下操作：

dst[i] = w | (src[i] << Shift);
w = src[i] >> (64 - Shift);

wherew是一些臨時變數，保存前一個陣列元素的右移值。

上面的這個解決方案簡單明了。但我需要更高效的東西，因為我有千兆位元組的資料。

理想情況下是使用一些 SIMD 指令，所以我正在尋找專家的 SIMD 建議。我需要為所有四種流行的指令集實作移位代碼 - SSE-SSE4.2 / AVX / AVX-2 / AVX-512。

但據我所知，例如對于 SSE2，只存在_mm_slli_si128()內在/指令，它僅按 8 的倍數移動（換句話說，位元組移動）。而且我需要按任意位大小進行移位，而不僅僅是位元組移位。

如果沒有 SIMD，我也可以通過使用shld reg, reg, reg指令一次移位 128 位，這允許進行 128 位移位。它在 MSVC 中作為內在的__shiftleft128()實作，并生成可以在這里看到的匯編代碼。

順便說一句，我需要所有 MSVC/GCC/CLang 的解決方案。

同樣在單回圈迭代中，我可以在順序操作中移動 4 或 8 個字，這將使用 CPU 流水線來加速多條指令的并行亂序執行。

如果需要，我的位向量可以與記憶體中的任意數量的位元組對齊，如果這有助于例如通過對齊讀/寫來提高 SIMD 速度。源和目標位向量存盤器也不同（非重疊）。

換句話說，我正在尋找有關如何在不同的 Intel CPU 上最有效（最高效）解決我的任務的所有建議。

注意，澄清一下，我實際上必須做幾個班次，而不僅僅是單班。我有大位向量X和數百個移位大小s0, s1, ..., sN，其中每個移位大小不同并且也可能很大（例如移位 100K 位），然后我想計算生成的大位向量Y = (X << s0) | (X << s1) | ... | (X << sN)。我只是將 StackOverflow 的問題簡化為移動單個向量。但可能這個關于原始任務的細節非常重要。

根據@Jake'Alquimista'LEE 的要求，我決定實作一個現成的玩具最小可重復示例，說明我想要做什么，計算輸入位向量的移位src或生成或最終dst位向量。這個例子根本沒有優化，只是我的任務如何解決的一個簡單的簡單變體。為簡單起見，此示例的輸入向量很小，而不是像我這樣的千兆位元組。這是一個玩具示例，我沒有檢查它是否正確解決了任務，它可能包含小錯誤：

在線試試吧！

#include <cstdint>
#include <vector>
#include <random>

#define bit_sizeof(x) (sizeof(x) * 8)

using u64 = uint64_t;
using T = u64;

int main() {
    std::mt19937_64 rng{123};

    // Random generate source bit vector
    std::vector<T> src(100'000);
    for (size_t i = 0; i < src.size();   i)
        src[i] = rng();

    size_t const src_bitsize = src.size() * bit_sizeof(T);

    // Destination bit vector, for example twice bigger in size
    std::vector<T> dst(src.size() * 2);

    // Random generate shifts
    std::vector<u64> shifts(200);
    for (size_t i = 0; i < shifts.size();   i)
        shifts[i] = rng() % src_bitsize;

    // Right-shift that handles overflow
    auto Shr = [](auto x, size_t s) {
        return s >= bit_sizeof(x) ? 0 : (x >> s);
    };

    // Do actual Shift-Ors
    for (auto orig_shift: shifts) {
        size_t const
            word_off = orig_shift / bit_sizeof(T),
            bit_off = orig_shift % bit_sizeof(T);

        if (word_off >= dst.size())
            continue;
        
        size_t const
            lim = std::min(src.size(), dst.size() - word_off);

        T w = 0;
        
        for (size_t i = 0; i < lim;   i) {
            dst[word_off   i] |= w | (src[i] << bit_off);
            w = Shr(src[i], bit_sizeof(T) - bit_off);
        }

        // Special case of handling for last word
        if (word_off   lim < dst.size())
            dst[word_off   lim] |= w;
    }
}

My real project's current code is different from toy example above. This project already solves correctly a real-world task. I just need to do extra optimizations. Some optimizations I already did, like using OpenMP to parallelize shift-or operations on all cores. Also as said in comments, I created specialized templated functions for each shift size, 64 functions in total, and choosing one of 64 functions to do actual shift-or. Each C function has compile time value of shift size, hence compiler does extra optimizations taking into account compile time values.

uj5u.com熱心網友回復：

您可以，甚至可能不需要顯式使用 SIMD 指令。目標編譯器 GCC、CLANG 和 MSVC 以及其他編譯器如 ICC 都支持自動矢量化。雖然手動優化的程式集可以勝過編譯器生成的向量化指令，但它通常更難實作，您可能需要針對不同架構的多個版本。導致高效自動向量化指令的通用代碼是一種可以跨許多平臺移植的解決方案。

例如一個簡單的 shiftvec 版本

void shiftvec(uint64_t* dst, uint64_t* src, int size, int shift)
{
    for (int i = 0; i < size;   i,  src,  dst)
    {
        *dst = ((*src)<<shift) | (*(src 1)>>(64-shift));
    }
}

使用最新的 GCC（或 CLANG 也可以）編譯，并-O3 -std=c 11 -mavx2在程式集的核心回圈中生成 SIMD 指令

.L5:
  vmovdqu ymm4, YMMWORD PTR [rsi rax]
  vmovdqu ymm5, YMMWORD PTR [rsi 8 rax]
  vpsllq ymm0, ymm4, xmm2
  vpsrlq ymm1, ymm5, xmm3
  vpor ymm0, ymm0, ymm1
  vmovdqu YMMWORD PTR [rdi rax], ymm0
  add rax, 32
  cmp rax, rdx
  jne .L5

在 Godbolt.org 上查看：https ://godbolt.org/z/5TxhqMhnK

如果您想在 dst 中組合多個班次，這也可以概括：

void shiftvec2(uint64_t* dst, uint64_t* src1, uint64_t* src2, int size1, int size2, int shift1, int shift2)
{
    int size = size1<size2 ? size1 : size2;
    for (int i = 0; i < size;   i,  src1,  src2,  dst)
    {
        *dst = ((*src1)<<shift1) | (*(src1 1)>>(64-shift1));
        *dst |= ((*src2)<<shift2) | (*(src2 1)>>(64-shift2)); 
    }
    for (int i = size; i < size1;   i,  src1,  dst)
    {
        *dst = ((*src1)<<shift1) | (*(src1 1)>>(64-shift1));        
    }
    for (int i = size; i < size2;   i,  src2,  dst)
    {
        *dst = ((*src2)<<shift2) | (*(src2 1)>>(64-shift2));
    }
}

編譯為核心回圈：

.L38:
  vmovdqu ymm7, YMMWORD PTR [rsi rcx]
  vpsllq ymm1, ymm7, xmm4
  vmovdqu ymm7, YMMWORD PTR [rsi 8 rcx]
  vpsrlq ymm0, ymm7, xmm6
  vpor ymm1, ymm1, ymm0
  vmovdqu YMMWORD PTR [rax rcx], ymm1
  vmovdqu ymm7, YMMWORD PTR [rdx rcx]
  vpsllq ymm0, ymm7, xmm3
  vmovdqu ymm7, YMMWORD PTR [rdx 8 rcx]
  vpsrlq ymm2, ymm7, xmm5
  vpor ymm0, ymm0, ymm2
  vpor ymm0, ymm0, ymm1
  vmovdqu YMMWORD PTR [rax rcx], ymm0
  add rcx, 32
  cmp r10, rcx
  jne .L38

在一個回圈中組合多個源將減少用于加載/寫入目標的記憶體帶寬總量。您可以組合的數量限制當然受可用暫存器的限制。注意，xmm2和xmm3用于shiftvec包含所述移位值，因此具有不同版本的編譯時已知的移值可以釋放那些暫存器。

此外__restrict，為每個指標使用（由 GCC、CLANG、MSVC 支持）將告訴編譯器范圍不重疊。

我最初在 MSVC 提供適當的自動矢量化代碼時遇到了問題，但似乎添加更多類似 SIMD 的結構將使其適用于所有三個所需的編譯器 GCC、CLANG 和 MSVC：

void shiftvec(uint64_t* __restrict dst, const uint64_t* __restrict src, int size, int shift)
{
    int i = 0;
    // MSVC: use steps of 2 for SSE, 4 for AVX2, 8 for AVX512
    for (; i 4 < size; i =4,dst =4,src =4)
    {
        for (int j = 0; j < 4;   j)
            *(dst j) = (*(src j))<<shift;
        for (int j = 0; j < 4;   j)
            *(dst j) |= (*(src 1)>>(64-shift));
    }
    for (; i < size;   i,  src,  dst)
    {
        *dst = ((*src)<<shift) | (*(src 1)>>(64-shift));
    }    
}

uj5u.com熱心網友回復：

我會嘗試依靠 x64 能力來讀取未對齊的地址，并且當星星正確（未）對齊時幾乎沒有明顯的損失。一個人只需要處理 (shift % 8) 或 (shift % 16) 的幾種情況——所有這些都可以用 SSE2 指令集，用零固定余數，對資料向量有一個未對齊的偏移量，并通過memcpy.

也就是說，內部回圈看起來像：

uint16_t const *ptr;
auto a = _mm_loadu_si128((__m128i*)ptr);
auto b = _mm_loadu_si128((__m128i*)(ptr   1);
a = _mm_srl_epi16(a, c);
b = _mm_sll_epi16(a, 16 - c);
_mm_storeu_si128((__m128i*)ptr, mm_or_si128(a,b));
ptr  = 8;

將這個回圈展開幾次，也許可以_mm_alignr_epi8在 SSE3 上使用來放寬記憶體帶寬（以及那些需要結合未對齊記憶體訪問結果的流水線階段）：

auto a0 = w; 
auto a1 = _mm_load_si128(m128ptr   1);
auto a2 = _mm_load_si128(m128ptr   2);
auto a3 = _mm_load_si128(m128ptr   3);
auto a4 = _mm_load_si128(m128ptr   4);
auto b0 = _mm_alignr_epi8(a1, a0, 2);
auto b1 = _mm_alignr_epi8(a2, a1, 2);
auto b2 = _mm_alignr_epi8(a3, a2, 2);
auto b3 = _mm_alignr_epi8(a4, a3, 2);
// ... do the computation as above ...
w = a4;   // rotate the context

uj5u.com熱心網友回復：

換句話說，我正在尋找有關如何在不同的 Intel CPU 上最有效（最高效）解決我的任務的所有建議。

效率的關鍵是懶惰。懶惰的關鍵是撒謊——假裝你改變了，實際上沒有做任何改變。

對于初始示例（僅用于說明概念），請考慮：

struct Thingy {
    int ignored_bits;
    uint64_t data[];
}

void shift_right(struct Thingy * thing, int count) {
    thing->ignored_bits  = count;
}

void shift_left(struct Thingy * thing, int count) {
    thing->ignored_bits -= count;
}

int get_bit(struct Thingy * thing, int bit_number) {
    bit_number  = thing->ignored_bits;
    return !!(thing->data[bit_number / 64] & (1 << bit_number % 64));
}

對于實際代碼，您需要關心各種細節 - 您可能希望從陣列開頭的備用位（和非零位ignored_bits）開始，以便您可以假裝右移；對于每個小的移位，您可能想要清除“移入”位（否則它將表現得像浮點數 - 例如(5.0 << 8) >> 8) == 5.0）；如果/何時ignored_bits超出某個范圍，您可能需要一個大的memcpy()；等等。

為了更多的樂趣；濫用低級記憶體管理 - 使用VirtualAlloc()（Windows）或mmap()（Linux）保留一個巨大的空間，然后將您的陣列放在空間的中間，然后根據需要在陣列的開始/結束處分配/釋放頁面；這樣您只需要memcpy()在原始位向左/向右“移動”數十億位之后即可。

當然，結果是它會使代碼的其他部分復雜化——例如，將 2 個位域 OR 在一起，您必須進行棘手的“獲取 A；移位 A 以匹配 B；結果 = A OR B”調整。這不是性能的交易破壞者。

uj5u.com熱心網友回復：

#include <cstdint>
#include <immintrin.h>

template<unsigned Shift>
void foo(uint64_t* __restrict pDst, const uint64_t* __restrict pSrc, intptr_t size)
{
    uint64_t* pSrc0, * pSrc1, * pSrc2, * pSrc3, * pDst0, * pDst1, * pDst2, * pDst3;
    __m256i prev, current;
    intptr_t i, stride;

    stride = size >> 2;
    i = stride;

    pSrc0 = pSrc;
    pSrc1 = pSrc   stride;
    pSrc2 = pSrc   2 * stride;
    pSrc2 = pSrc   3 * stride;

    pDst0 = pDst;
    pDst1 = pDst   stride;
    pDst2 = pDst   2 * stride;
    pDst3 = pDst   3 * stride;

    prev = _mm256_set_epi64x(0, pSrc1[-1], pSrc2[-1], pSrc3[-1]);

    while (i--)
    {
        current = _mm256_set_epi64x(*pSrc0  , *pSrc1  , *pSrc2  , *pSrc3  );
        prev = _mm256_srli_epi64(prev, 64 - Shift);
        prev = _mm256_or_si256(prev, _mm256_slli_epi64(current, Shift));
        *pDst0   = _mm256_extract_epi64(prev, 3);
        *pDst1   = _mm256_extract_epi64(prev, 2);
        *pDst2   = _mm256_extract_epi64(prev, 1);
        *pDst3   = _mm256_extract_epi64(prev, 0);

        prev = current;
    }
}

您可以在 AVX2 上一次對最多四個 64 位元素執行操作（在 AVX512 上最多為八個）

如果 size 不是 4 的倍數，則最多可以處理 3 個剩余的。

PS：自動矢量化從來都不是一個合適的解決方案。

uj5u.com熱心網友回復：

不，你不能

NEON 和 AVX(512) 都支持高達 64 位元素的桶形移位操作。

但是，您可以ext使用 NEON 和alignrAVX上的指令將整個 128 位向量“移位”n 位元組（8 位）。

并且您應該避免使用 vector 類來提高性能，因為它只不過是對性能不利的鏈表。

標籤：c performance simd sse avx

上一篇：Openpyxl、Pandas或兩者

下一篇：如何在Python中查找日期時間范圍