向量化嵌套循环 - SIMD答案

【问题标题】：Vectorizing Nested Loop - SIMD向量化嵌套循环 - SIMD
【发布时间】：2012-05-09 00:44:01
【问题描述】：

有谁知道使用 SIMD 对这样的东西进行矢量化：

for(size_t i = 0; i < refSeq.length() / 4; i++){

    for(size_t j = 0; j<otherSeq.length(); j++){
    if(refSeq[i] == otherSeq[j]){
        if(i == 0 || j == 0)
            L[i][j] = 1;
       else
        L[i][j] = L[i-1][j-1] + 1;
    }
       else
        L[i][j] = 0;
    }
}

【问题讨论】：

当您说“SIMD”时，您需要指定您正在谈论的 CPU 系列，例如x86（在这种情况下您需要指定可以假设的 SSE/AVX 级别）、PowerPC (AltiVec)、POWER (VMX/VSX)、ARM (Neon)、Cell 等。
另外，refSeq[]、otherSeq[] 和 L[][] 的数据类型是什么？
您确实非常坚持使这个最长的子字符串算法并行:) 再次 - 数据依赖。 SIMD 适用于独立的数据块。在这里，循环中有if（坏）和if（更糟）。您需要重新设计算法以使用屏蔽而不是分支，我不确定它是否会运行得更快。
它适用于 x86，所有这些都是字符串。 SSE intel 的说明。
@Hristo 我是，我尝试了一切，但我能找到的只是一些有趣的例子，没有更复杂的更大的例子。 :) 是的，掩蔽......虽然一个:) 谢谢:)

标签： c++ parallel-processing x86 vectorization simd

【解决方案1】：

由于某种原因，@arunmoezhi 的解决方案没有自动矢量化。 https://godbolt.org/z/KM4fh8TY7

.L8:
        mov     edx, DWORD PTR [rax+4+r8]
        cmp     DWORD PTR [rcx], edx
        je      .L5
        mov     DWORD PTR [rsi+4+rax], 0
        add     rax, 4
        cmp     rdi, rax
        jne     .L8

这是一个带有 eve 库的实现：https://godbolt.org/z/j8Kva8caP

主循环看起来像这样 (avx2) 只是展开：

.LBB0_5:                                # =>This Inner Loop Header: Depth=1
        vmovdqu ymm2, ymmword ptr [rax + rcx]
        vpcmpeqd        ymm3, ymm0, ymmword ptr [rsi + rcx]
        vpsubd  ymm2, ymm2, ymm1
        vpand   ymm2, ymm3, ymm2
        vmovdqu ymmword ptr [rdi + rcx], ymm2
        lea     rbx, [rsi + rcx]
        add     rbx, 32
        add     rcx, 32
        cmp     rbx, r11
        jne     .LBB0_5

【讨论】：

【解决方案2】：

让我尝试提出一个解决方案。最初计算 L[i][0] 和 L[0][j] 的值。现在从 i=1 和 j=1 开始迭代。现在可以删除循环每次迭代中对 i==0 或 j==0 的检查。这样做的另一个优点是，对于一行中每次迭代中的每个 L[i][j]，L[i-1][j-1] 的值都是可用的。现在，假设向量寄存器是否可以容纳数组的 4 个元素。现在我们可以加载 refSeq、otherSeq、L(previous row) 和 L(current row) 的 4 个元素。理论上我们现在得到了矢量化。我假设自动矢量化器不会识别这一点。所以我们必须手动完成。如果我错了，请纠正我。

for(size_t i=0;i<refSeq.length()/4;i++)
{
    if(refSeq[i]==otherSeq[0])
        L[i][0]=1;
    else
        L[i][0]=0;
}
for(size_t j=0; j<otherSeq.length();j++)
{
    if(refSeq[0]==otherSeq[j])
        L[0][j]=1;
    else
        L[0][j]=0;
}

for(size_t i=1;i<refSeq.length()/4;i++)
{
    for(size_t j=1; j<otherSeq.length();j++)
    {
        if(refSeq[i]==otherSeq[j])
            L[i][j] = L[i-1][j-1] + 1;
        else
            L[i][j]=0;
    }
}

一个缺点是，现在我们正在加载前一行，无论 refSeq[i] 是否等于 otherSeq[j] 与原始代码一样，其中仅当序列相等时才访问对角元素。

【讨论】：

【解决方案3】：

这是一个动态规划问题，海峡前向实现具有太多的数据依赖性，不适合 SIMD 计算。

但是，如果您将算法从逐行迭代更改为对角迭代，则可以并行计算整个对角线。见下图。

下面的“伪”代码使用具有 1 个额外行/列的矩阵来简化“内部”计算。这个额外的行/列在每次对角迭代之前初始化。

int i, j, k;
for (k = 1; ; k++) {
    int minI = k > refLen ? k - refLen : 1;
    int maxI = k > otherLen ? otherLen : k - 1;

    for (i = maxI; i >= minI; ) {
        j = k - i;

        // vectorized calculation 256 bit (AVX2)
        if (i >= 32 && otherLen - j >= 32) {
            // calculate 32 values of the diagonal with SIMD
            i -= 32;
            continue;
        }

        // vectorized calculation 128 bit (SSE)
        if (i >= 16 && otherLen - j >= 16) {
            // calculate 16 values of the diagonal with SIMD
            i -= 16;
            continue;
        }

        // scalar calculation
        if (refSeq[i - 1] == otherSeq[j - 1]) {
            L[i][j] = L[i - 1][j - 1] + 1;
        } else {
            L[i][j] = 0;
        }
        i--;
    }

    if (k == otherLen + refLen) {
        break;
    }

    // initialize next j-endpoint in diagonal
    if (k <= refLen) {
        L[0][k] = 0;
    }
    // initialize next i-endpoint in diagonal
    if (k <= otherLen) {
        L[k][0] = 0;
    }
}

不确定您是否需要用于计算的实际 SIMD 指令，或者只是知道如何并行化/矢量化计算。

【讨论】：