AVX2 赢者通吃差异搜索
AVX2 Winner-Take-All Disparity Search
我正在使用 AVX2 优化视差估计算法的 "winner-take-all" 部分。我的标量例程是准确的，但在 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 );
视差 Space 图像 (DSI) 的大小为 HxWxD (320x240x48)，为了更好的内存访问，我将其水平放置，以便每一行的大小为 WxD。
视差 Space 图像具有每像素匹配成本。这个聚合
用一个简单的盒子过滤器来制作另一个大小完全相同的图像，
但是成本总和，比如 3x3 或 5x5 window。这种平滑使得
结果更多'robust'。当我使用 asPtr 访问时，我正在编制索引
进入此汇总成本图像。
此外，为了节省不必要的计算，我一直在开始
并以掩码半径偏移的行结束。这个mask radius就是半径
我的人口普查面具。我可以做一些花哨的边界反射，但它是
更简单，更快，只是为了不为这个边界的差异而烦恼。
这当然也适用于开头和结尾的 cols，但是弄乱了
当我将整个算法强制为 运行 时，此处的索引并不好
在列数为 16 的倍数的图像上（例如 QVGA：320x240），这样我
可以简单地索引并使用 SIMD 命中所有内容（无残差标量处理）。
此外，如果您认为我的代码一团糟，我鼓励您查看
高度优化的 OpenCV 立体算法。我发现它们是不可能的，而且几乎没有用到它们。
我的代码编译但在 运行 时失败。我正在使用 VS 2012 Express Update 4。当我 运行 使用调试器时，我无法获得任何见解。我对使用内在函数比较陌生，所以我不确定在调试时应该看到什么信息、寄存器数量、__m256i 变量是否应该可见等。
听从下面的评论建议，我通过使用更智能的索引将标量时间从 ~14 提高到 ~8。我的 CPU 是 i7-4980HQ，我在同一文件的其他地方成功使用了 AVX2 内在函数。
我仍然没有找到问题所在，但我确实看到了您可能想要更改的一些内容。不过，您没有检查 _mm_malloc 的 return 值。如果它失败了，那就可以解释了。 （也许它不喜欢分配 32 字节对齐的内存？）
如果您 运行 在内存检查器或其他东西下检查您的代码，那么它可能不喜欢从未初始化的内存中读取 dispVals。 （即使掩码是全一，_mm256_maskstore_epi32 也可能算作读取-修改-写入。）
运行 您的代码在调试器下运行并找出问题所在。 "runtime error"意义不大。
_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。 Intel CPUs 可以在运行时将 16b 内存位置符号扩展到 gp 寄存器中，而不会造成速度损失（movsz 与 mov 一样快），所以很有可能。只需将 dRPtr 数组声明为 uint16_t。然后，您根本不需要矢量代码中的符号扩展内容（更不用说在您的内部循环中了！）。希望 _mm256_extracti128_si256( mask, 0 ) 编译为空，因为你想要的 128 已经是 low128，所以只需使用 reg 作为 vmovsx 的 src，但仍然。
您还可以通过不首先加载来保存指令（和融合域 uop）。 （除非编译器足够聪明，不会省略 vmovdqu 并将 vpminuw 与内存操作数一起使用，即使您使用了 load intrinsic）。
所以我在想这样的事情：
// 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 uops（blendv 为 2），加上循环开销。它应该 运行 在 1.5 个周期/迭代（不计算循环开销）在 haswell 上，但 dep 链将它限制为每 2 个周期一次迭代。 （假设展开以消除循环开销）。并行执行两个 dep 链可以解决此问题，并且与展开循环具有相同的效果：减少循环开销。
嗯，实际上是在 Haswell，
pminuw 可以在 p1/p5 上 运行。 （以及 p2/p3 上的负载部分）pcmpgtw 可以在 p1/p5vpblendvb 是 p5 的 2 微指令。padduw 可以在 p1/p5movdqa reg,reg 可以在 p0/p1/p5 上 运行（并且可能根本不需要执行单元）。展开应该消除 minCosts = newMinCosts 的任何开销，因为编译器可以从下一次迭代的第一个循环体的右寄存器中的最后一个展开的循环体中得到 newMinCosts。- fused
sub / jge（循环计数器）可以在 p6 上 运行。 （在 dvec 上使用 PTEST + jcc 会更慢）。当不与 jcc. 融合时，add/sub 可以在 p0/p1/p5/p6 上 运行
好的，所以实际上循环每次迭代需要 2.5 个周期，受限于只能在 p1/p5 上 运行 的指令。按 2 或 4 展开将减少循环 / movdqa 开销。由于 Haswell 每个时钟可以发出 4 微指令，因此它可以更有效地排队微指令以进行乱序执行，因为循环不会有超高的迭代次数。 （48 是你的例子。）有很多 uops 排队会让 CPU 在离开循环后有事可做，并隐藏缓存未命中等的任何延迟。
_mm256_min_epu16 (PMINUW) 是另一个循环携带的依赖链。将它与内存操作数一起使用使其具有 3 或 4 个周期的延迟。然而，指令的加载部分可以在地址已知后立即开始，因此将加载折叠到修改操作中以利用微融合不会使 dep 链比使用单独的加载更长或更短。
有时您需要为未对齐的数据使用单独的加载（AVX 删除了内存操作数的对齐要求）。我们受执行单元的限制比 4 微指令/时钟问题限制更多，因此使用专用加载指令可能没问题。
在您进行特定于平台的优化之前，可以执行大量可移植的优化。提取循环不变量，将索引乘法转换为增量加法等...
这可能不准确，但大致理解了：
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" 实际上是对顺序数据进行操作的变量。 ...下一步将相应地重新排序循环，如果您仍需要进一步优化，请应用特定于平台的内在函数。
我正在使用 AVX2 优化视差估计算法的 "winner-take-all" 部分。我的标量例程是准确的，但在 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 );
视差 Space 图像 (DSI) 的大小为 HxWxD (320x240x48)，为了更好的内存访问，我将其水平放置，以便每一行的大小为 WxD。
视差 Space 图像具有每像素匹配成本。这个聚合 用一个简单的盒子过滤器来制作另一个大小完全相同的图像， 但是成本总和，比如 3x3 或 5x5 window。这种平滑使得 结果更多'robust'。当我使用 asPtr 访问时，我正在编制索引 进入此汇总成本图像。
此外，为了节省不必要的计算，我一直在开始 并以掩码半径偏移的行结束。这个mask radius就是半径 我的人口普查面具。我可以做一些花哨的边界反射，但它是 更简单，更快，只是为了不为这个边界的差异而烦恼。 这当然也适用于开头和结尾的 cols，但是弄乱了 当我将整个算法强制为 运行 时，此处的索引并不好 在列数为 16 的倍数的图像上（例如 QVGA：320x240），这样我 可以简单地索引并使用 SIMD 命中所有内容（无残差标量处理）。
此外，如果您认为我的代码一团糟，我鼓励您查看 高度优化的 OpenCV 立体算法。我发现它们是不可能的，而且几乎没有用到它们。
我的代码编译但在 运行 时失败。我正在使用 VS 2012 Express Update 4。当我 运行 使用调试器时，我无法获得任何见解。我对使用内在函数比较陌生，所以我不确定在调试时应该看到什么信息、寄存器数量、__m256i 变量是否应该可见等。
听从下面的评论建议，我通过使用更智能的索引将标量时间从 ~14 提高到 ~8。我的 CPU 是 i7-4980HQ，我在同一文件的其他地方成功使用了 AVX2 内在函数。
我仍然没有找到问题所在，但我确实看到了您可能想要更改的一些内容。不过，您没有检查 _mm_malloc 的 return 值。如果它失败了，那就可以解释了。 （也许它不喜欢分配 32 字节对齐的内存？）
如果您 运行 在内存检查器或其他东西下检查您的代码，那么它可能不喜欢从未初始化的内存中读取 dispVals。 （即使掩码是全一，_mm256_maskstore_epi32 也可能算作读取-修改-写入。）
运行 您的代码在调试器下运行并找出问题所在。 "runtime error"意义不大。
_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。 Intel CPUs 可以在运行时将 16b 内存位置符号扩展到 gp 寄存器中，而不会造成速度损失（movsz 与 mov 一样快），所以很有可能。只需将 dRPtr 数组声明为 uint16_t。然后，您根本不需要矢量代码中的符号扩展内容（更不用说在您的内部循环中了！）。希望 _mm256_extracti128_si256( mask, 0 ) 编译为空，因为你想要的 128 已经是 low128，所以只需使用 reg 作为 vmovsx 的 src，但仍然。
您还可以通过不首先加载来保存指令（和融合域 uop）。 （除非编译器足够聪明，不会省略 vmovdqu 并将 vpminuw 与内存操作数一起使用，即使您使用了 load intrinsic）。
所以我在想这样的事情：
// 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 uops（blendv 为 2），加上循环开销。它应该 运行 在 1.5 个周期/迭代（不计算循环开销）在 haswell 上，但 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会更慢）。当不与jcc. 融合时，
add/sub 可以在 p0/p1/p5/p6 上 运行
好的，所以实际上循环每次迭代需要 2.5 个周期，受限于只能在 p1/p5 上 运行 的指令。按 2 或 4 展开将减少循环 / movdqa 开销。由于 Haswell 每个时钟可以发出 4 微指令，因此它可以更有效地排队微指令以进行乱序执行，因为循环不会有超高的迭代次数。 （48 是你的例子。）有很多 uops 排队会让 CPU 在离开循环后有事可做，并隐藏缓存未命中等的任何延迟。
_mm256_min_epu16 (PMINUW) 是另一个循环携带的依赖链。将它与内存操作数一起使用使其具有 3 或 4 个周期的延迟。然而，指令的加载部分可以在地址已知后立即开始，因此将加载折叠到修改操作中以利用微融合不会使 dep 链比使用单独的加载更长或更短。
有时您需要为未对齐的数据使用单独的加载（AVX 删除了内存操作数的对齐要求）。我们受执行单元的限制比 4 微指令/时钟问题限制更多，因此使用专用加载指令可能没问题。
在您进行特定于平台的优化之前，可以执行大量可移植的优化。提取循环不变量，将索引乘法转换为增量加法等...
这可能不准确，但大致理解了：
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" 实际上是对顺序数据进行操作的变量。 ...下一步将相应地重新排序循环，如果您仍需要进一步优化，请应用特定于平台的内在函数。