AVX2 赢家通吃视差搜索

Posted 2023-02-23

【中文标题】AVX2 赢家通吃视差搜索

【英文标题】：AVX2 Winner-Take-All Disparity Search

【发布时间】：2015-08-02 13:13:56

【问题描述】：

我正在使用 AVX2 优化视差估计算法的“赢家通吃”部分。我的标量例程是准确的，但在 QVGA 分辨率和 48 个视差下，运行时间在我的笔记本电脑上慢到 14 毫秒，令人失望。我创建了 LR 和 RL 视差图像，但为了简单起见，我将只包含用于 RL 搜索的代码。

我的标量例程：

int MAXCOST = 32000;
for (int i = maskRadius; i < rstep-maskRadius; i++) 

    // WTA "RL" Search:
    for (int j = maskRadius; j+maskRadius < cstep; j++) 
        int minCost = MAXCOST;
        int minDisp = 0;
        for (int d = 0; d < numDisp && j+d < cstep; d++) 
            if (asPtr[(i*numDisp*cstep)+(d*cstep)+j] < minCost) 
                minCost = asPtr[(i*numDisp*cstep)+(d*cstep)+j];
                minDisp = d;
            
        
        dRPtr[(i*cstep)+j] = minDisp;
    

我尝试使用 AVX2：

int MAXCOST = 32000;
int* dispVals = (int*) _mm_malloc( sizeof(int32_t)*16, 32 );

for (int i = maskRadius; i < rstep-maskRadius; i++) 

    // WTA "RL" Search AVX2:
    for( int j = 0; j < cstep-16; j+=16) 

        __m256i minCosts = _mm256_set1_epi16( MAXCOST );
        __m128i loMask   = _mm_setzero_si128();
        __m128i hiMask   = _mm_setzero_si128();

        for (int d = 0; d < numDisp && j+d < cstep; d++) 
            // Grab 16 costs to compare
            __m256i costs = _mm256_loadu_si256((__m256i*) (asPtr[(i*numDisp*cstep)+(d*cstep)+j]));

            // Get the new minimums
            __m256i newMinCosts = _mm256_min_epu16( minCosts, costs );

            // Compare new mins to old to build mask to store minDisps
            __m256i mask   = _mm256_cmpgt_epi16( minCosts, newMinCosts );
            __m128i loMask = _mm256_extracti128_si256( mask, 0 );
            __m128i hiMask = _mm256_extracti128_si256( mask, 1 );
            // Sign extend to 32bits
            __m256i loMask32 = _mm256_cvtepi16_epi32( loMask );
            __m256i hiMask32 = _mm256_cvtepi16_epi32( hiMask );

            __m256i currentDisp = _mm256_set1_epi32( d );
            // store min disps with mask
            _mm256_maskstore_epi32( dispVals, loMask32, currentDisp );    // RT error, why?
            _mm256_maskstore_epi32( dispVals+8, hiMask32, currentDisp );  // RT error, why?

            // Set minCosts to newMinCosts
            minCosts = newMinCosts;
        

        // Write the WTA minimums one-by-one to the RL disparity image
        int index = (i*cstep)+j;
        for( int k = 0; k < 16; k++ ) 
            dRPtr[index+k] = dispVals[k];
        
    

_mm_free( dispVals );

视差空间图像 (DSI) 的大小为 HxWxD (320x240x48)，我将其水平布局以便更好地访问内存，这样每一行的大小都是 WxD。

视差空间图像具有每像素匹配成本。这个汇总 用一个简单的盒子过滤器来制作另一张完全相同大小的图像， 但成本总和超过 3x3 或 5x5 的窗口。这种平滑使 结果更“稳健”。当我使用 asPtr 访问时，我正在索引 到这个汇总成本图像中。

另外，为了节省不必要的计算，我已经开始 并以掩码半径偏移的行结束。这个遮罩半径就是半径 我的人口普查面具。我可以做一些花哨的边界反射，但它是 更简单，更快，只是为了不打扰这个边界的差异。 这当然也适用于开头和结尾的 cols，但是会弄乱 当我强制我的整个算法只运行时，这里的索引不好 在列是 16 的倍数的图像上（例如 QVGA：320x240），这样我 可以简单地索引并使用 SIMD 命中所有内容（无残留标量处理）。

另外，如果您认为我的代码一团糟，我鼓励您查看 高度优化的 OpenCV 立体算法。我发现它们是不可能的，并且几乎没有使用它们。

我的代码编译但在运行时失败。我正在使用 VS 2012 Express Update 4。当我使用调试器运行时，我无法获得任何见解。我对使用内在函数比较陌生，所以我不确定调试时应该看到哪些信息、寄存器数量、__m256i 变量是否应该可见等。

听取下面的评论建议，我通过使用更智能的索引将标量时间从 ~14 提高到 ~8。我的 CPU 是 i7-4980HQ，我在同一个文件的其他地方成功使用了 AVX2 内部函数。

【问题讨论】：

你确定这个循环是正确的吗：for( int j = 0; j < cstep-16; j+=16) - 看起来你可能错过了最后 16 个值？

为了加快速度，特别是在最内层的循环中，重复的表达式应该被提取并且只执行一次。例如，d * cstep 表达式只需要在每个内部循环中执行一次，与 i * numDisp * cstep 类似

现在你有一个工作版本并不重要，但我会试图通过查看反汇编来找出导致运行时错误的确切原因，以及错误发生时的寄存器状态触发。因此，您可以判断用于_mm256_maskstore_epi32 的指针是否正确，以及指向的内存是否可通过调试器访问。可能是错误的指针（C 指针数学给出了意外的结果？）或 malloc 失败的问题。 “运行时错误”可以是任何东西，从非法指令到错误指针，所以它没有那么有用。

【参考方案1】：

我仍然没有发现问题，但我确实看到了一些您可能想要更改的内容。不过，您没有检查_mm_malloc 的返回值。如果它失败了，那就可以解释了。 （也许它不喜欢分配 32 字节对齐的内存？）

如果您在内存检查器或其他东西下运行代码，那么它可能不喜欢从未初始化的内存中读取dispVals。 （_mm256_maskstore_epi32 可能算作读-修改-写，即使掩码是全一。）

在调试器下运行您的代码并找出问题所在。 “运行时错误”意义不大。

_mm_set1* 函数很慢。 VPBROADCASTD 需要它的源在内存或向量 reg，而不是 GP reg，因此编译器可以 movd 从 GP reg 到向量 reg 然后广播，或者存储到内存然后广播。无论如何，这样做会更快

const __m256i add1 = _mm256_set1_epi32( 1 );
__m256i dvec = _mm256_setzero_si256();
for (d;d...;d++) 
    dvec = _mm256_add_epi32(dvec, add1);

其他内容： 如果您不将内部循环的每次迭代都存储到内存中，这可能会运行得更快。使用混合指令 (_mm256_blendv_epi8) 或类似指令来更新最小成本的位移向量。 Blend = 带有注册目标的蒙版移动。

此外，您的位移值应该适合 16b 整数，因此在找到它们之前不要将它们符号扩展为 32b。英特尔 CPU 可以即时将 16b 内存位置签名扩展至 gp 寄存器，而不会降低速度（movsz 与mov 一样快），所以很有可能。只需将您的dRPtr 数组声明为uint16_t。然后，您根本不需要矢量代码中的符号扩展内容（更不用说在您的内部循环中了！）。希望_mm256_extracti128_si256( mask, 0 ) 编译为空，因为您想要的 128 已经是 low128，所以只需使用 reg 作为vmovsx 的 src，但仍然。

您还可以通过不先加载来保存指令（和融合域 uop）。 （除非编译器足够聪明，不会省略 vmovdqu 并将 vpminuw 与内存操作数一起使用，即使您使用了 load 内部函数）。

所以我在想这样的事情：

// totally untested, didn't even check that this compiles.
for(i)  for(j) 
// inner loop, compiler can hoist these constants.
const __m256i add1 = _mm256_set1_epi16( 1 );
__m256i dvec = _mm256_setzero_si256();
__m256i minCosts = _mm256_set1_epi16( MAXCOST );
__m256i minDisps = _mm256_setzero_si256();

for (int d=0 ; d < numDisp && j+d < cstep ;
     d++, dvec = _mm256_add_epi16(dvec, add1))

    __m256i newMinCosts = _mm256_min_epu16( minCosts, asPtr[(i*numDisp*cstep)+(d*cstep)+j]) );
    __m256i mask   = _mm256_cmpgt_epi16( minCosts, newMinCosts );
    minDisps = _mm256_blendv_epi8(minDisps, dvec, mask); // 2 uops, latency=2
    minCosts = newMinCosts;


// put sign extension here if making dRPtr uint16_t isn't an option.
int index = (i*cstep)+j;
_mm256_storeu_si256 (dRPtr + index, __m256i minDisps);

如果有两个并行的依赖链：minCosts0 / minDisps0 和 minCosts1 / minDisps1，您可能会获得更好的性能，然后在最后将它们组合起来。 minDisps 是一个循环携带的依赖项，但循环只有 5 条指令（包括 vpadd，它看起来像循环开销，但不能通过展开来减少）。它们解码为 6 微指令（blendv 为 2），加上循环开销。它应该在 haswell 上以 1.5 个周期/迭代（不计算循环开销）运行，但 dep 链会将其限制为每 2 个周期进行一次迭代。 （假设展开以摆脱循环开销）。并行执行两个 dep 链可以解决此问题，并且与展开循环具有相同的效果：减少循环开销。

嗯，实际上是在 Haswell 上，

pminuw 可以在 p1/p5 上运行。 （以及 p2/p3 上的负载部分） pcmpgtw 可以在 p1/p5 上运行 vpblendvb 是 p5 的 2 微秒。 padduw 可以在 p1/p5 上运行 movdqa reg,reg 可以在 p0/p1/p5 上运行（并且可能根本不需要执行单元）。展开应该可以消除 minCosts = newMinCosts 的任何开销，因为编译器可以在下一次迭代的第一个循环体的正确寄存器中的最后一个展开的循环体中以 newMinCosts 结束。 fused sub / jge（循环计数器）可以在 p6 上运行。 （在 dvec 上使用 PTEST + jcc 会更慢）。 add/sub 可以在未与 jcc 融合时在 p0/p1/p5/p6 上运行。

好的，所以实际上循环每次迭代需要 2.5 个周期，受只能在 p1/p5 上运行的指令的限制。展开 2 或 4 将减少循环 /movdqa 开销。由于 Haswell 可以每个时钟发出 4 个 uops，因此它可以更有效地将 uops 排队以进行乱序执行，因为循环不会有超多的迭代次数。 （48 是你的例子。）有很多 uops 排队会让 CPU 在离开循环后有事可做，并隐藏缓存未命中等的任何延迟。

_mm256_min_epu16 (PMINUW) 是另一个循环携带的依赖链。将它与内存操作数一起使用会导致 3 或 4 个周期的延迟。但是，一旦地址已知，指令的加载部分就可以开始，因此将加载折叠到修改操作中以利用微融合不会使 dep 链比使用单独的加载更长或更短。

有时您需要对未对齐的数据使用单独的加载（AVX 删除了内存操作数的对齐要求）。我们受执行单元的限制比 4 uop / 时钟发出限制更多，因此使用专用加载指令可能没问题。

source for insn ports / latencies.

【讨论】：

非常感谢！从~14ms 到~1ms。巨额储蓄。从来没有发现我原始代码的错误。

感谢您编辑您的问题以描述您的代码实际在做什么。我认为它正在索引某种成本数组，但是那里的详细信息可能对其他人有用。很高兴我能帮助收紧内环。

【参考方案2】：

在进行特定于平台的优化之前，可以执行许多可移植的优化。提取循环不变量，将索引乘法转换为增量加法等...

这可能不准确，但大致了解：

int MAXCOST = 32000, numDispXcstep = numDisp*cstep;
for (int i = maskRadius; i < rstep - maskRadius; i+=numDispXcstep) 
    for (int j = maskRadius; j < cstep - maskRadius; j++) 
        int minCost = MAXCOST, minDisp = 0;
        for (int d = 0; d < numDispXcstep - j; d+=cstep) 
            if (asPtr[i+j+d] < minCost) 
                minCost = asPtr[i+j+d];
                minDisp = d;
            
        
        dRPtr[i/numDisp+j] = minDisp;
    

一旦你这样做了，实际上发生了什么就变得很明显了。看起来“i”是最大的步骤，其次是“d”，“j”实际上是对顺序数据进行操作的变量。 ...下一步将相应地重新排序循环，如果您仍需要进一步优化，请应用特定于平台的内在函数。

【讨论】：

在循环之外提升循环不变表达式在这种情况下肯定会使代码更易于阅读，但编译器通常可以为您做到这一点。但是，如果表达式涉及局部变量以外的任何内容，则编译器不能假设它们不会改变。无论如何，数组索引之前确实是一团糟。除了检查它是否在连续的源数据上矢量化之外，我没有花太多时间试图解开它的答案。

@technosaurus 我改进了索引并大大提高了速度，当我有另一个空闲时间时，我会编辑问题以显示这一点。鉴于数据是一个 Hx(WxD) 数组，因此“图像”是并排而不是堆叠的，我认为去 i->j->d 而不将 minCosts 保存到临时是有意义的大小为 W 的数组（对于所有 j）。这样临时数据永远不会离开寄存器。否则，绝对地，索引 i->d->j 更合适，这也是我在计算成本并将其保存到 DSI 时所做的。我计划在标量例程上尝试 i->d->j 并报告。

@DerLuftmensch 我将顺序数据保持在一起的原因是内在操作往往会得到更好的优化......如果您查看它们的 typedef 和数据对齐，很明显您可以简单地做一个强制转换，而不是使用内部函数来初始化变量。