我应该何时以及如何在我的 simd 例程中执行浮点转换？
Posted
技术标签:
【中文标题】我应该何时以及如何在我的 simd 例程中执行浮点转换？【英文标题】：When and how should I perform floating point conversion in my simd routine? 【发布时间】：2020-09-16 04:29:34 【问题描述】：我正在计算 2 张图像的双向（水平和垂直）前缀和（扫描），产生像素总和、平方和以及两幅图像的叉积。所有的计算都是在 32 位整数中完成的，最后我需要将 32 位整数转换为双精度数以计算窗口函数中两个图像的均值、方差和协方差。
首先——这是执行此操作的最佳方式吗？我可以用双精度构建整个前缀和数组，并且没有转换步骤。
第二——如果这是正确的方法，我是否会从使用打包的双 simd 操作中获得很多好处？我只能放心地假设我一次会拿出 2 个单位。
第三——我应该将数据单元打包在一起还是保留其当前的平面格式？ [平面格式是一种像素被“组件”分解的格式。如果您获得 32 位 RGBA 输入，即 8 位 R、8 位 G、8 位 B 和 8 位 A，则打包格式将为 RGBARGBA，而平面格式将为 RRRRRRRRRRRRR....GGGGGGGGGGGG....BBBBB.. ...AAAAA...等等。]
以下是我迄今为止完成的与该主题相关的三个功能。前 2 个是标量版本，因此更容易阅读和理解正在发生的事情。第三个是功能 1 的当前 SIMD 实现。第四个功能（缺少且尚未完成）是这个问题的主题，可能是第二个的 SIMD 实现。
std::unique_ptr<uint32_t[],boost::alignment::aligned_delete> computeSumMatrixForwardScalar2PassAll(uint8_t const* pImgData1, uint8_t const* pImgData2,
unsigned width, unsigned height)
using namespace simdpp;
std::unique_ptr<uint32_t[], boost::alignment::aligned_delete> sumArray((uint32_t*)boost::alignment::aligned_alloc(64, 5*width*height*sizeof(uint32_t)));
auto pSumArray = sumArray.get();
BOOST_ALIGN_ASSUME_ALIGNED(pImgData1, 64);
BOOST_ALIGN_ASSUME_ALIGNED(pImgData2, 64);
BOOST_ALIGN_ASSUME_ALIGNED(pSumArray, 64);
//#pramga omp parallel for private(h) shared(pImgData, pSumArray, w )
#pragma omp for simd
for (unsigned h = 0; h < height; ++h)
uint32_t lastValX = 0;
uint32_t lastValY = 0;
uint32_t lastValXX = 0;
uint32_t lastValYY = 0;
uint32_t lastValXY = 0;
for (unsigned w = 0; w < width; ++w)
uint32_t imgValX = pImgData1[h * width + w];
uint32_t newValX = lastValX + imgValX;
uint32_t newValXX = lastValXX + imgValX*imgValX;
uint32_t imgValY = pImgData2[h*width + w];
uint32_t newValY = lastValY + imgValY;
uint32_t newValYY = lastValYY + imgValY*imgValY;
uint32_t newValXY = lastValXY + imgValX*imgValY;
pSumArray[h*width + w]= newValX;
pSumArray[width*height+h*width + w] = newValY;
pSumArray[2*width*height+ h*width + w] = newValXX;
pSumArray[3*width*height+h*width + w] = newValYY;
pSumArray[4*width*height+h*width + w] = newValXY;
lastValX = newValX;
lastValXX = newValXX;
lastValY = newValY;
lastValYY = newValYY;
lastValXY = newValXY;
for (unsigned i = 0; i < 5; ++i)
for (unsigned h = 0; h+1 < height; ++h)
for (unsigned w = 0; w < width; ++w)
uint32_t above = pSumArray[i*width*height + h * width + w];
uint32_t current = pSumArray[i*width*height+ (h+1) *width +w];
pSumArray[i*width*height + (h+1) * width +w]= above+current;
return sumArray;
第二个：SSIM 转换——注意不同的语言——因为我还没有完成它的 C++ 实现。注意，它在里面调用了weberSumMatrix，和上面的函数是一样的。
export function weberSsim(
pixels1: ImageMatrix,
pixels2: ImageMatrix,
options: Options
): MSSIMMatrix
// console.time("weber ssim");
const bitDepth, k1, k2, windowSize = options
const L = (1 << bitDepth) - 1
const c1 = k1 * L * (k1 * L)
const c2 = k2 * L * (k2 * L)
const windowSquared = windowSize * windowSize
const pixels1Data = pixels1.data;
const pixels2Data = pixels2.data;
const width = pixels1.width;
const height = pixels1.height;
// Produces exactly the same output as the C++ prefix sum above.
const sumMatrix = weberSumMatrix(pixels1Data, pixels2Data, width, height);
const windowHeight = height-windowSize;
const windowWidth = width-windowSize;
const imageSize = width*height;
const ssims = new Array(windowHeight*windowWidth);
// lets handle w = 0 h = 0 first and initialize mssim
let cumulativeSsim;
const reciprocalWindowSquared = 1 / windowSquared;
const windowOffset = windowSize - 1;
let bottomOffset = windowOffset*width;
const meanx = (sumMatrix[bottomOffset+ windowOffset]) * reciprocalWindowSquared;
const meany = (
sumMatrix[imageSize + bottomOffset+ windowOffset]) * reciprocalWindowSquared;
const varx = (
sumMatrix[2*imageSize + bottomOffset+ windowOffset]) * reciprocalWindowSquared - meanx*meanx ;
const vary = (
sumMatrix[3*imageSize + bottomOffset+ windowOffset]) * reciprocalWindowSquared - meany*meany;
const cov = (
sumMatrix[4*imageSize + bottomOffset+ windowOffset]) * reciprocalWindowSquared - meanx*meany;
const na = 2 * meanx * meany + c1
const nb = 2 * cov + c2
const da = meanx * meanx + meany * meany + c1
const db = varx + vary + c2
const ssim = (na * nb) / (da * db)
ssims[0] = ssim
// mssim = ssim
cumulativeSsim = ssim;
// next handle all of the h = 0, w > 0 cases first
for (let w = 1; w < windowWidth; ++w)
// in h =0 cases, there is no top left or top right
let leftOffset = w - 1;
const rightx = sumMatrix[bottomOffset+leftOffset];
const leftx = sumMatrix[bottomOffset+(windowOffset+w)];
const meanx = (leftx-rightx)* reciprocalWindowSquared;
const righty= sumMatrix[imageSize + bottomOffset+ leftOffset];
const lefty = sumMatrix[imageSize + bottomOffset+ (windowOffset+w)];
const meany = (lefty-righty) * reciprocalWindowSquared;
const rightxx = sumMatrix[2*imageSize + bottomOffset+leftOffset];
const leftxx = sumMatrix[2*imageSize + bottomOffset+ (windowOffset+w)];
const varx = (leftxx-rightxx) * reciprocalWindowSquared - meanx*meanx ;
const rightyy = sumMatrix[3*imageSize + bottomOffset+leftOffset];
const leftyy = sumMatrix[3*imageSize + bottomOffset+ (windowOffset+w)]
const vary = (leftyy - rightyy) * reciprocalWindowSquared - meany*meany;
const rightxy = sumMatrix[4*imageSize + bottomOffset+leftOffset];
const leftxy = sumMatrix[4*imageSize + bottomOffset+ (windowOffset+w)];
const cov = (leftxy-rightxy) * reciprocalWindowSquared - meanx*meany;
const na = 2 * meanx * meany + c1
const nb = 2 * cov + c2
const da = meanx * meanx + meany * meany + c1
const db = varx + vary + c2
const ssim = (na * nb) / (da *db)
ssims[w] = ssim
// mssim = mssim + (ssim - mssim) / (i + 1)
cumulativeSsim += ssim;
const windowOffset = windowSize - 1;
// There will be lots of branch misses if we don't split the w==0 and h==0 cases
for (let h = 1; h < windowHeight; ++h)
// now the w=0 on each line
let bottomOffset = (h+windowSize-1)*width;
let topOffset = (h-1)*width;
// since there is no left side we can skip two operations
const topx = sumMatrix[topOffset+ windowOffset];
const bottomx = sumMatrix[bottomOffset+ windowOffset];
const meanx = (bottomx - topx) * reciprocalWindowSquared;
const topy = sumMatrix[imageSize + topOffset+ windowOffset];
const bottomy = sumMatrix[imageSize + bottomOffset+ windowOffset];
const meany = (bottomy - topy) * reciprocalWindowSquared;
const topxx = sumMatrix[2*imageSize + topOffset+ windowOffset];
const bottomxx = sumMatrix[2*imageSize + bottomOffset+ windowOffset];
const varx = (bottomxx-topxx) * reciprocalWindowSquared - meanx*meanx ;
const topyy = sumMatrix[3*imageSize + topOffset+ windowOffset];
const bottomyy = sumMatrix[3*imageSize + bottomOffset+ windowOffset];
const vary = (bottomyy-topyy) * reciprocalWindowSquared - meany*meany;
const topxy = sumMatrix[4*imageSize + topOffset+ windowOffset];
const bottomxy = sumMatrix[4*imageSize + bottomOffset+ windowOffset];
const cov = (bottomxy-topxy) * reciprocalWindowSquared - meanx*meany;
const na = 2 * meanx * meany + c1
const nb = 2 * cov + c2
const da = meanx * meanx + meany * meany + c1
const db = varx + vary + c2
const ssim = (na * nb) / (da *db)
ssims[h*windowWidth] = ssim
// mssim = mssim + (ssim - mssim) / (i + 1)
cumulativeSsim += ssim;
for (let w = 1; w < windowWidth; ++w)
// add top left sub top right sub bottom left add bottom right
const rightOffset = w + windowSize - 1;
const leftOffset = w - 1;
const meanx = (sumMatrix[topOffset + leftOffset]
- sumMatrix[topOffset+ rightOffset]
- sumMatrix[bottomOffset+leftOffset]
+ sumMatrix[bottomOffset+ rightOffset]) * reciprocalWindowSquared;
const meany = (sumMatrix[imageSize+ topOffset + leftOffset]
- sumMatrix[imageSize + topOffset+ rightOffset]
- sumMatrix[imageSize + bottomOffset+leftOffset]
+ sumMatrix[imageSize + bottomOffset+ rightOffset]) * reciprocalWindowSquared;
const varx = (sumMatrix[2*imageSize+ topOffset + leftOffset]
- sumMatrix[2*imageSize + topOffset+ rightOffset]
- sumMatrix[2*imageSize + bottomOffset+leftOffset]
+ sumMatrix[2*imageSize + bottomOffset+ rightOffset]) * reciprocalWindowSquared - meanx*meanx ;
const vary = (sumMatrix[3*imageSize+ topOffset + leftOffset]
- sumMatrix[3*imageSize + topOffset+ rightOffset]
- sumMatrix[3*imageSize + bottomOffset+leftOffset]
+ sumMatrix[3*imageSize + bottomOffset+ rightOffset]) * reciprocalWindowSquared - meany*meany;
const cov = (sumMatrix[4*imageSize+ topOffset + leftOffset]
- sumMatrix[4*imageSize + topOffset+ rightOffset]
- sumMatrix[4*imageSize + bottomOffset+leftOffset]
+ sumMatrix[4*imageSize + bottomOffset+ rightOffset]) * reciprocalWindowSquared - meanx*meany;
const na = 2 * meanx * meany + c1
const nb = 2 * cov + c2
const da = meanx * meanx + meany * meany + c1
const db = varx + vary + c2
const ssim = (na * nb) / (da * db)
ssims[h*windowWidth+w] = ssim
cumulativeSsim += ssim;
// mssim = mssim + (ssim - mssim) / (i + 1)
const mssim = cumulativeSsim / (windowHeight*windowWidth);
return data: ssims, width, height, mssim
第三个：当前 SIMD 前缀总和。
std::unique_ptr<uint32_t[],boost::alignment::aligned_delete> computeSumMatrixForwardSimd2PassAll(uint8_t const* pImgData1, uint8_t const* pImgData2,
unsigned width, unsigned height)
using namespace simdpp;
std::unique_ptr<uint32_t[], boost::alignment::aligned_delete> sumArray((uint32_t*)boost::alignment::aligned_alloc(64, 5*width*height*sizeof(uint32_t)));
auto pSumArray = sumArray.get();
BOOST_ALIGN_ASSUME_ALIGNED(pImgData1, 64);
BOOST_ALIGN_ASSUME_ALIGNED(pImgData2, 64);
BOOST_ALIGN_ASSUME_ALIGNED(pSumArray, 64);
//#pramga omp parallel for private(h) shared(pImgData, pSumArray, w )
uint32x4 zero = make_zero();
for (unsigned h = 0; h < height; ++h)
uint32x4 lastValSplatX = zero;
uint32x4 lastValSplatY = zero;
uint32x4 lastValSplatXX = zero;
uint32x4 lastValSplatYY = zero;
uint32x4 lastValSplatXY = zero;
for (unsigned w = 0; w < width; w += 16)
// starting left value
// previous line values..
prefetch_read(pImgData1+(w+1)*64);
prefetch_read(pImgData2+(w+1)*64);
uint32v4 imgDataX = to_uint32(uint8x16(load(pImgData1 + h * width + w)));
uint32v4 imgDataY = to_uint32(uint8x16(load(pImgData2 + h * width + w)));
static_assert(uint32v4::vec_length == 4);
static_assert(sizeof(uint32v4::base_vector_type::native_type) == 16);
for (unsigned i = 0 ; i < uint32v4::vec_length; ++i)
// a_0 a_1 a_2 a_3
uint32v4::base_vector_type x = imgDataX.vec(i);
uint32v4::base_vector_type y = imgDataY.vec(i);
uint32v4::base_vector_type xx = mul_lo(x,x);
uint32v4::base_vector_type yy = mul_lo(y, y);
uint32v4::base_vector_type xy = mul_lo(x, y);
// a_0 a_0+a_1 a_1+a_2 a_2+a_3
x = add(x, move4_r<1>(x));
x = add(x, move4_r<2>(x));
x = add(x, lastValSplatX);
lastValSplatX = permute4<3,3,3,3>(x);
store(pSumArray+h*width+w+i*4, x);
y = add(y, move4_r<1>(y));
y = add(y, move4_r<2>(y));
y = add(y, lastValSplatY);
lastValSplatY = permute4<3,3,3,3>(y);
store(width*height+pSumArray+h*width+w+i*4, y);
xx = add(xx, move4_r<1>(xx));
xx = add(xx, move4_r<2>(xx));
xx = add(xx, lastValSplatXX);
lastValSplatXX = permute4<3,3,3,3>(xx);
store(2*width*height+pSumArray+h*width+w+i*4, xx);
yy = add(yy, move4_r<1>(yy));
yy = add(yy, move4_r<2>(yy));
yy = add(yy, lastValSplatYY);
lastValSplatYY = permute4<3,3,3,3>(yy);
store(3*width*height+pSumArray+h*width+w+i*4, yy);
xy = add(xy, move4_r<1>(xy));
xy = add(xy, move4_r<2>(xy));
xy = add(xy, lastValSplatXY);
lastValSplatXY = permute4<3,3,3,3>(xy);
store(4*width*height+pSumArray+h*width+w+i*4, xy);
// 16 bit 8s for grins...
// a_0 a_1 a_2 a_3 a_4 a_5 a_6 a_7
// a_0 a_0+a_1 a_1+a_2 a_2+a_3 a_3+a_4 a_4+a_5 a_5+a_6 a_6+a_7 (>>1)
// d d -a_0 -a_0+a_1 -a_0+a_1+a_2 -a_0+a_1+a_2+a_3 -a_0+a_1+a_2+a_3+a_4 - -a_0+a_1+a_2+a_3+a_4+a_5
// d d shuffle shuffle shuffle+add shuffle+add shuffle+add+shuffle+add shuffle+add+shuffle+add
// a_0 a_1 a_2 a_3 a_4 a_5 a_6 a_7
// a_0 a_1 a_2 a_3 a_4 a_5
// a_1 a_2+a_0 a_3+a_1-a_0 a_4+a_2-a1-a0 a_5+a_3-a_2-a_0 a_6+a_4...
for (unsigned i = 0; i < 5; ++i)
for (unsigned h = 0; h+1 < height; ++h)
for (unsigned w = 0; w < width; w += 16)
uint32v4 above = load(i*width*height +pSumArray + h * width + w);
uint32v4 current = load(i*width*height+pSumArray +(h+1) * width +w);
store(i*width*height+pSumArray +(h+1) * width +w, add(above,current));
return sumArray;
【问题讨论】：
您同时提出三个不同的问题。问题 1 和 2 是相关的，问题 3 完全不同，更适合单独的问题。如果您确实发布它，请务必准确解释您所拥有的“平面格式”究竟是什么。事实上，很难说出你在问什么。 【参考方案1】：首先——这是执行此操作的最佳方式吗？我可以用双精度构建整个前缀和数组，并且没有转换步骤。
通常，整数计算比每个向量具有相同元素数的浮点计算要快得多。例如，在 Skylake 上，paddd 的延迟为 1，倒数吞吐量为 0.33 个周期，addps - 4 和 0.5。对于整数计算，如果有乘法，有时会因为必须在不同元素大小之间转换或组合乘积的上半部分和下半部分而产生开销。但通常它最终仍然比 FP 快，如果您能够合并 pmaddwd 或 pmaddubsw 指令（它们是 FMA 的 INT 版本），您也可以减少它，或者您可以丢弃产品的下半部分或上半部分。
说到 FMA，FP 因为 AVX2 具有 FMA 指令的优势，它允许每次乘法免费进行一次加法。在您的情况下这是否有益取决于您的算法和输入数据，但如果您有整数输入，我仍然希望尽可能地对其进行 INT 计算。
INT 计算相对于 FP 的另一个优势是您可以对较小的元素进行计算，这意味着您可以在每条指令中处理更多数据。当然，这只有在您的输入和算法允许的情况下才有可能，但特别是在图像处理中，情况往往如此。
关于这部分的最后一点是，您应该考虑在算法的每个阶段需要处理的数据量。 INT 到 FP 的转换不是免费的，因此您需要转换的数据越少越好。如果您的中间数据量少于输入，那么推迟转换将是有益的。
第二——如果这是正确的方法，我是否会从使用打包的双 simd 操作中获得很多好处？我只能放心地假设我一次会拿出 2 个单位。
嗯，每条指令 2 个双打比一个要好，所以在我的书中这是值得的。与标量相比，您还可以对向量执行更精细的操作，例如 FMA、屏蔽、混合、最小/最大值（尽管编译器甚至可能在标量代码上为您生成一些指令）。如果您在运行时检测到 AVX，您可以机会性地使吞吐量翻倍。
此外，您应该考虑 32 位浮点精度是否足以满足您的情况。您可能不需要双精度结果，并且通过使用一些 FP 技术，您可以减少使用 32 位 FP 数学累积的错误，并且仍然具有比 64 位计算更好的性能。
第三——我应该将数据单元打包在一起还是保留其当前的平面格式？
一般来说，您应该更喜欢平面输入数据。一般而言，SIMD，尤其是 SSE/AVX 更适合垂直操作（即在不同向量的对应元素之间执行操作时）而不是水平操作（当在同一向量的元素之间执行操作时）。对于打包输入，您可能必须对输入数据执行解包和混洗，这将增加开销。现代 CPU 能够跟踪多个从/到内存的读取或写入流，因此硬件预取器应该能够处理不同数据平面上的线性内存访问。
【讨论】：
关于您在中间数据量较少时延迟转换的评论，很好。在我的情况下，如果你仔细观察，你会发现我的情况正好相反。中间数据是原始数据数量的 5 倍。由于我需要将它们全部转换为 uint32，因此 5x 一直持续到存在到 double 的转换。就我而言，我已经确定我可以使用浮动到 16x16 窗口大小的点。在那之后，它似乎是唯一的选择。但是，float 不是前缀总和表的选项。 @DanWeber 如果您的中间数据大于输入，那么您将不会从减少转换操作的数量中受益，但您仍然可以从输入处理的某些部分是 INT 的事实中受益，这比 FP 快。这取决于你在 INT 能做多少。 关于转换为 uint32_t 的需要，请注意，在累积阶段，您将添加 8 位整数。总和始终适合 16 位，这意味着您可以在每条指令中累积两倍的元素。并且“lastValXX + imgValX*imgValX”表达式映射到 pmaddwd 指令。如果您可以管理其中一个向量以使其具有签名的 8 位元素，甚至可以使用 pmaddubsw。 使用 pmaddubsw 的一种方法是在其中一个输入为 128 或更大的情况下添加单独的 shift+addition。请注意，x * (128 + y) = x * 128 + x * y = (x 【参考方案2】：对于第 3 部分——至少，英特尔® 64 和 IA-32 架构优化参考手册指定平面格式通常更理想。所以我可以一次性过关。看起来我可以为高达 16x16 (255*256**2) 的窗口执行浮点值而不是双精度值，因为这在溢出之前是可以表示的。从这个角度来看，我得到了每个向量 4 个浮点数而不是 2 个双精度数的好处。这很有帮助。看来我可以从四个浮点数中得到一个滚动总和，并在最后计算它。因此，如果我根据最大窗口大小限制函数，我可以确定何时使用浮点数以及何时使用双精度数。实现并没有太大的不同。
我想在我尝试看看之前我不会知道它的真实性能。
【讨论】：
【参考方案3】：关于第 1 部分，根据整数的数量及其范围，您可能会使用 64 位整数而不是双精度数。通常稍微快一点，尤其是延迟，也更精确。这是 SSE2 中的一些片段。
// Return acc64 += a32. The a32 is viewed as uint32_t lanes, the accumulator is uint64_t
inline __m128i integerAdd( __m128i a32, __m128i acc64 )
const __m128i low = _mm_and_si128( a32, _mm_set1_epi64x( UINT_MAX ) );
acc64 = _mm_add_epi64( acc64, low );
const __m128i high = _mm_srli_epi64( a32, 32 );
acc64 = _mm_add_epi64( acc64, high );
return acc64;
// Compute a32 * b32, add to the accumulator. The two inputs are viewed as uint32_t lanes, the accumulator is uint64_t
inline __m128i integerFma( __m128i a32, __m128i b32, __m128i acc64 )
const __m128i low = _mm_mul_epu32( a32, b32 );
a32 = _mm_srli_si128( a32, 4 );
b32 = _mm_srli_si128( b32, 4 );
const __m128i high = _mm_mul_epu32( a32, b32 );
acc64 = _mm_add_epi64( acc64, low );
acc64 = _mm_add_epi64( acc64, high );
return acc64;
// Add both 64-bit lanes of the accumulator, convert to double
inline double accumulatedValue( __m128i acc64 )
acc64 = _mm_add_epi64( acc64, _mm_unpackhi_epi64( acc64, acc64 ) );
const uint64_t v = (uint64_t)_mm_cvtsi128_si64( acc64 );
return (double)v;
更新：平均值为sum(x) / count。你不需要分数来计算整数的总和，如果你有很多值，你只需要 64 位整数来避免溢出。最后一条指令只需要非整数。
对于您需要sum((x-mean)^2) 的方差。这可以重写为sum(x^2) + mean * ( mean - 2 * sum(x) )。现在很明显，您只需要对输入数据进行一次传递，计算平方和和值的总和。如果输入太多且太大，计算平方和可能会溢出 64 位整数。但对于许多实际应用来说，这很好，例如如果您的输入在 [0 .. 65535] 范围内，64 位整数仅在 40 亿次样本后溢出，您可能没有那么多。
【讨论】：
我认为这没有帮助，因为我对 uint32_t 没有任何问题（即使溢出或翻转）。如果我必须转换为浮点类型的数字来执行最终计算，那么使用 64 位整数对我没有任何帮助，除非您提出一种更有效地进行小数计算的方法。 如何保证mean是整数或整数而不损失精度？我不担心溢出，因为舍入算法的工作方式可以处理额外的 65535 个值，该值大于任何合理的窗口大小。 @DanWeber 平均值不是整数。但是，元素之和及其平方和都保证为整数。处理完完整的输入数据集后，才需要双重数学。在两个累加器上调用accumulatedValue()，然后进行双重数学运算。性能不再重要，因为这些只是 2 个数字，即使您输入了十亿个元素。
你明白我正在处理图像并且这是一个双向前缀和吗？以上是关于我应该何时以及如何在我的 simd 例程中执行浮点转换？的主要内容，如果未能解决你的问题，请参考以下文章
如何从使用 pg_cron 运行的清除例程中链接 VACUUM？