修改函数以使用 SSE 内在函数

Posted 2023-02-16

【中文标题】修改函数以使用 SSE 内在函数

【英文标题】：Modifying a function to use SSE intrinsics

【发布时间】：2015-01-30 14:04:46

【问题描述】：

我正在尝试计算激进的近似值：sqrt(i + sqrt(i + sqrt(i + ...))) 使用 SSE 以从矢量化中获得加速（我还读到 SIMD 平方根函数的运行速度比先天 FPU 平方快大约 4.7 倍-根函数）。但是，我在矢量化版本中获得相同的功能时遇到问题；我得到的值不正确，我不确定

我原来的功能是这样的：

template <typename T>
T CalculateRadical( T tValue, T tEps = std::numeric_limits<T>::epsilon() )

    static std::unordered_map<T,T> setResults;

    auto it = setResults.find( tValue );
    if( it != setResults.end() )
    
        return it->second;
    

    T tPrev = std::sqrt(tValue + std::sqrt(tValue)), tCurr = std::sqrt(tValue + tPrev);

    // Keep iterating until we get convergence:
    while( std::abs( tPrev - tCurr ) > tEps )
    
        tPrev = tCurr;
        tCurr = std::sqrt(tValue + tPrev);
    

    setResults.insert( std::make_pair( tValue, tCurr ) );
    return tCurr;

我写的 SIMD 等价物（当这个模板函数用T = float 实例化并给出tEps = 0.0005f）是：

// SSE intrinsics hard-coded function:
__m128 CalculateRadicals( __m128 values )

    static std::unordered_map<float, __m128> setResults;

    // Store our epsilon as a vector for quick comparison:
    __declspec(align(16)) float flEps[4] =  0.0005f, 0.0005f, 0.0005f, 0.0005f ;
    __m128 eps = _mm_load_ps( flEps );

    union U 
        __m128 vec;
        float flArray[4];
    ;

    U u;
    u.vec = values;

    float flFirstVal = u.flArray[0];
    auto it = setResults.find( flFirstVal );
    if( it != setResults.end( ) )
    
        return it->second;
    

    __m128 prev = _mm_sqrt_ps( _mm_add_ps( values, _mm_sqrt_ps( values ) ) );
    __m128 curr = _mm_sqrt_ps( _mm_add_ps( values, prev ) );

    while( _mm_movemask_ps( _mm_cmplt_ps( _mm_sub_ps( curr, prev ), eps ) ) != 0xF )
    
        prev = curr;
        curr = _mm_sqrt_ps( _mm_add_ps( values, prev ) );
    

    setResults.insert( std::make_pair( flFirstVal, curr ) );
    return curr;

我正在使用以下代码循环调用该函数：

long long N;
std::cin >> N;

float flExpectation = 0.0f;
long long iMultipleOf4 = (N / 4LL) * 4LL;
for( long long i = iMultipleOf4; i > 0LL; i -= 4LL )

    __declspec(align(16)) float flArray[4] =  static_cast<float>(i - 3), static_cast<float>(i - 2), static_cast<float>(i - 1), static_cast<float>(i) ;
    __m128 arg = _mm_load_ps( flArray );
    __m128 vec = CalculateRadicals( arg );

    float flSum = Sum( vec );
    flExpectation += flSum;


for( long long i = iMultipleOf4; i < N; ++i )

    flExpectation += CalculateRadical( static_cast<float>(i), 0.0005f );


flExpectation /= N;

对于输入5，我得到以下输出：

With SSE version: 2.20873
With FPU verison: 1.69647

差异来自哪里，我在 SIMD 等效项中做错了什么？

---

编辑：我意识到 Sum 函数在这里是相关的：

float Sum( __m128 vec1 )

    float flTemp[4];
    _mm_storeu_ps( flTemp, vec1 );
    return flTemp[0] + flTemp[1] + flTemp[2] + flTemp[3];

【问题讨论】：

我想知道，您比较中的 0xF 是否被符号扩展为与 0x0F 不同的值？

@BrianCain 将值更改为0x0F 不幸的是并没有改变结果！不过感谢您的建议！

为了稳健性，这个比较std::abs( tPrev - tCurr ) > tEps 可能应该关心isnan。不过可能不是你的问题。

@BrianCain 是的，但我认为合同可以与参数传递者有关；只要0 <= tValue <= FL_MAX isnan 应该不是问题吧？

当然，除非isnan(tValue)。足够公平，调用者应该知道isnan(tValue)这个函数是否永远不会返回。

【参考方案1】：

SSE 内在函数有时可能相当乏味...

但不是在这里。你只是搞砸了你的循环：

for( long long i = iMultipleOf4; i > 0LL; i -= 4LL )

我怀疑它是否符合您的预期。如果iMultipleOf4 是 4，那么您的函数将使用 4,3,2,1 而不是 0 进行计算。然后您的第二个循环使用 4 重做计算。

这两个函数对我来说给出了相同的结果，而循环在更正后给出了相同的flExpectation。尽管仍然存在细微差别，可能是因为 FPU 在计算方式上略有不同。

【讨论】：

非常感谢；这似乎解决了我的问题！然而;我很好奇为什么值会有如此大的差异，我没想到前5个有效数字会有差异；你知道我能做些什么来帮助纠正这个问题吗？ （我需要精确到 10^-3 的解）？

@Shaktal 我检查了（在我的机器上）到底发生了什么，似乎std::sqrt 正在使用非向量双精度 SSE 进行计算。您的 SSE 代码使用矢量单精度，因此您可能会得到不太精确的结果。匹配计算的唯一方法可能是使用双精度。

@Shaktal 与您的问题无关，自 2011 年以来的 CPU 具有 AVX 指令，可启用 256 位宽的向量，而不是您使用的 128 位。如果你想使用它，你应该在运行时检查 AVX 支持。此外，2015 年，新的 CPU 将支持 AVX-512。

那么您是说std::sqrt 会将我传递给它的float 值转换为double，然后使用内部双精度sqrt FPU 函数，然后再转换回来？

@Shaktal 是的。当我编译你的代码时，程序集将浮点数加载到 SSE 寄存器中，转换为双精度，SSE 双精度平方，然后转换为单精度。这可能因编译器、版本和架构而异，但这就是我得到的。 （我在 Windows 8.1 amd64 上使用过 VS2013）