4. SIMD向量点积

任务目标

实现一个使用NEON（ARM）或SSE（x86）指令集加速的浮点向量点积函数，并与普通标量实现对比验证正确性。本案例以x86 SSE为例，代码可在ARM上通过替换内建函数适配NEON。

核心实现

#include <cstdint>
#include <cstring>
#include <vector>
#include <random>
#include <iostream>
#include <cmath>

// SSE 向量点积（4路浮点）
float dot_sse(const float* a, const float* b, size_t n) {
    __m128 sum = _mm_setzero_ps();
    size_t i = 0;
    // 主循环：每次处理4个浮点
    for (; i + 4 <= n; i += 4) {
        __m128 va = _mm_loadu_ps(a + i);
        __m128 vb = _mm_loadu_ps(b + i);
        __m128 prod = _mm_mul_ps(va, vb);
        sum = _mm_add_ps(sum, prod);
    }
    // 水平相加
    float tmp[4];
    _mm_storeu_ps(tmp, sum);
    float result = tmp[0] + tmp[1] + tmp[2] + tmp[3];
    // 处理剩余元素
    for (; i < n; ++i) {
        result += a[i] * b[i];
    }
    return result;
}

// 标量点积（参考实现）
float dot_scalar(const float* a, const float* b, size_t n) {
    float sum = 0.0f;
    for (size_t i = 0; i < n; ++i) {
        sum += a[i] * b[i];
    }
    return sum;
}

int main() {
    const size_t N = 1000;
    std::vector<float> a(N), b(N);
    std::mt19937 rng(42);
    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
    for (size_t i = 0; i < N; ++i) {
        a[i] = dist(rng);
        b[i] = dist(rng);
    }

    float sse_result = dot_sse(a.data(), b.data(), N);
    float scalar_result = dot_scalar(a.data(), b.data(), N);

    // 验证正确性
    float diff = std::fabs(sse_result - scalar_result);
    float tolerance = 1e-4f * std::max(1.0f, std::fabs(scalar_result));
    bool pass = diff < tolerance;

    std::cout << "SSE result:    " << sse_result << "\n";
    std::cout << "Scalar result: " << scalar_result << "\n";
    std::cout << "Difference:    " << diff << "\n";
    std::cout << (pass ? "PASS" : "FAIL") << std::endl;

    return pass ? 0 : 1;
}

编译与运行

# x86 平台（需支持SSE，现代CPU默认支持）
g++ -O2 -msse4.2 -o dot_sse dot_sse.cpp
./dot_sse

验证方法

• 程序输出SSE与标量两种实现的计算结果及差值。
• 当差值小于1e-4 × max(1, |标量结果|)时输出PASS，退出码为0；否则输出FAIL，退出码为1。
• 由于浮点运算顺序不同，两种实现结果存在微小舍入误差属正常现象，容差已考虑该因素。

关键说明

• 使用_mm_loadu_ps支持未对齐内存访问，避免对齐要求带来的额外约束。
• 主循环每次处理4个float，尾部不足4个元素由标量循环处理，保证任意长度输入的正确性。
• 若需移植到ARM平台，将__m128替换为float32x4_t，_mm_loadu_ps替换为vld1q_f32，_mm_mul_ps替换为vmulq_f32，_mm_add_ps替换为vaddq_f32，水平相加使用vaddvq_f32即可。