SIMD-高性能计算入门

SIMD

SIMD(Single Instruction Multiple Data)即单指令流多数据流，是一种采用一个控制器来控制多个处理器，同时对一组数据（又称“数据向量”）中的每一个分别执行相同的操作从而实现空间上的并行性的技术。常见的已经实现的指令集如下：

指令集	寄存器位宽	寄存器名	一次能处理的 float (32-bit)	一次能处理的 double (64-bit)
SSE	128 bit	xmm0-15	4 × float	2 × double
AVX	256 bit	ymm0-15	8 × float	4 × double
AVX2	256 bit	ymm0-15	8 × float	4 × double
AVX-512	512 bit	zmm0-31	16 × float	8 × double

Intrinsics 101：头文件 & 命名规则

ISA	头文件	前缀	例
SSE	`<xmmintrin.h>`	`_mm_*`	`_mm_add_ps`
AVX	`<immintrin.h>`	`_mm256_*`	`_mm256_add_ps`
AVX-512	`<immintrin.h>`	`_mm512_*`	`_mm512_add_ps`

只需 #include <immintrin.h> 即可覆盖 SSE → AVX-512。

最常用的 Intrinsics 指令

类别	Scalar	AVX2 (256-bit)	语义
置零	`float z = 0.0f`	`_mm256_setzero_ps()`	8×0
加载	`a[i]`	`_mm256_load_ps(ptr)`	从内存 → ymm
存储	`a[i] = v`	`_mm256_store_ps(ptr, v)`	ymm → 内存
加法	`a + b`	`_mm256_add_ps(va, vb)`	8 并行加
乘法	`a * b`	`_mm256_mul_ps(va, vb)`	8 并行乘
水平求和	手写 for	`_mm256_hadd_ps` → `_mm_add_ps` → `_mm_store_ss`	把 8 个数累加成一个

实践

这里使用 Intel Intrinsics，可以直接使用 C/C++ 代码调用 SIMD 指令集。下面使用一个向量 L2 距离计算的例子来展示 SIMD 的使用方式、可能遇到的浮点误差以及性能提升。

首先需要引入头文件

#include <immintrin.h> // AVX, AVX2, AVX-512
#include <iostream>
#include <vector>
#include <chrono>
#include <random>
#include <cmath>

定义计时函数并且定义生成随机向量函数以便测试：

// ---------- 计时 ----------
template <typename Func>
double timeit(Func f, const float* a, const float* b, size_t n, int rounds = 10) {
    double best = 1e20;
    for (int r = 0; r < rounds; ++r) {
        auto t0 = high_resolution_clock::now();
        volatile auto res = f(a, b, n);
        auto t1 = high_resolution_clock::now();
        best = min(best, duration<double, milli>(t1 - t0).count());
    }
    return best;
}

// 生成随机浮点向量
vector<float> rand_vec(size_t n) {
    vector<float> v(n);
    mt19937 rng{random_device{}()};
    uniform_real_distribution<float> dist(-1.f, 1.f);
    for (auto& x : v) x = dist(rng);
    return v;
}

下面分别使用 float 和 double 来进行数值计算对比：

// ---------- 标量版本 (float 累加) ----------
float l2_scalar(const float* a, const float* b, size_t n) {
    float sum = 0.0f;
    for (size_t i = 0; i < n; ++i) {
        float d = a[i] - b[i];
        sum += d * d;
    }
    return sum;
}
// ---------- AVX2 版本 (float 累加) ----------
float l2_avx2(const float* a, const float* b, size_t n) {
    __m256 vsum = _mm256_setzero_ps();
    size_t i = 0;
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(&a[i]);
        __m256 vb = _mm256_loadu_ps(&b[i]);
        __m256 diff = _mm256_sub_ps(va, vb);
        __m256 sqr  = _mm256_mul_ps(diff, diff);
        vsum = _mm256_add_ps(vsum, sqr);
    }

    // 把 256-bit 拆成两个 128-bit
    __m128 vlow  = _mm256_castps256_ps128(vsum);   // 低 128
    __m128 vhigh = _mm256_extractf128_ps(vsum, 1); // 高 128
    __m128 vsum128 = _mm_add_ps(vlow, vhigh);

    // 水平加和 128-bit (4个float)
    __m128 shuf = _mm_movehdup_ps(vsum128);        // (x3,x3,x1,x1)
    __m128 sums = _mm_add_ps(vsum128, shuf);
    shuf        = _mm_movehl_ps(shuf, sums);
    sums        = _mm_add_ss(sums, shuf);

    float final_sum;
    _mm_store_ss(&final_sum, sums);
    return final_sum;
}

// ---------- 标量版本 (double 累加) ----------
double l2_scalar_double(const float* a, const float* b, size_t n) {
    double sum = 0.0;
    for (size_t i = 0; i < n; ++i) {
        double d = (double)a[i] - (double)b[i];
        sum += d * d;
    }
    return sum;
}

// ---------- AVX2 版本 (double 累加) ----------
double l2_avx2_double(const float* a, const float* b, size_t n) {
    __m256d vsum = _mm256_setzero_pd(); // double 累加
    size_t i = 0;

    // 每次处理 8 个 float → 8 个 diff²
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(&a[i]);
        __m256 vb = _mm256_loadu_ps(&b[i]);
        __m256 diff = _mm256_sub_ps(va, vb);
        __m256 sqr  = _mm256_mul_ps(diff, diff);

        // 拆成两个 128 位 float，转成 double
        __m128 lo = _mm256_castps256_ps128(sqr);
        __m128 hi = _mm256_extractf128_ps(sqr, 1);

        __m256d dlo = _mm256_cvtps_pd(lo); // 4 float → 4 double
        __m256d dhi = _mm256_cvtps_pd(hi); // 4 float → 4 double

        vsum = _mm256_add_pd(vsum, dlo);
        vsum = _mm256_add_pd(vsum, dhi);
    }

    // 水平加和 vsum (4 doubles)
    __m128d vlow  = _mm256_castpd256_pd128(vsum);
    __m128d vhigh = _mm256_extractf128_pd(vsum, 1);
    __m128d sum128 = _mm_add_pd(vlow, vhigh);
    __m128d shuf   = _mm_unpackhi_pd(sum128, sum128);
    __m128d sums   = _mm_add_sd(sum128, shuf);

    double final_sum;
    _mm_store_sd(&final_sum, sums);

    // 处理尾部 (不足 8 个 float)
    for (; i < n; i++) {
        double d = (double)a[i] - (double)b[i];
        final_sum += d * d;
    }
    return final_sum;
}

测试 main 函数

int main() {
    const size_t N = (1 << 20); // 1M
    auto a = rand_vec(N);
    auto b = rand_vec(N);

    // 测试 double 类型
    double ref = l2_scalar_double(a.data(), b.data(), N);
    double avx2_val = l2_avx2_double(a.data(), b.data(), N);

    // 校验 double 类型正确性
    double abs_diff = fabs(avx2_val - ref);
    double rel_diff = abs_diff / ref;

    cout << "ref = " << ref << ", avx2 = " << avx2_val
         << ", abs diff = " << abs_diff
         << ", rel diff = " << rel_diff << '\n';

    if (rel_diff > 1e-4) { 
        cerr << "SIMD 结果错误！\n";
    }

    // 测试 float 类型
    float ref_f = l2_scalar(a.data(), b.data(), N);
    float avx2_val_f = l2_avx2(a.data(), b.data(), N);

    // 校验 float 类型正确性
    float abs_diff_f = fabs(avx2_val_f - ref_f);
    float rel_diff_f = abs_diff_f / ref_f;

    cout << "ref(float) = " << ref_f << ", avx2(float) = " << avx2_val_f
         << ", abs diff = " << abs_diff_f
         << ", rel diff = " << rel_diff_f << '\n';

    if (rel_diff_f > 1e-4) { 
        cerr << "SIMD float 结果错误！\n";
    }

    // 计时
    double t_scalar_f = timeit(l2_scalar, a.data(), b.data(), N);
    double t_scalar_d = timeit(l2_scalar_double, a.data(), b.data(), N);
    double t_avx2_f   = timeit(l2_avx2, a.data(), b.data(), N);
    double t_avx2_d   = timeit(l2_avx2_double, a.data(), b.data(), N);

    cout << "===== L2 distance (1M floats) =====\n";
    cout << "scalar(float)  : " << t_scalar_f << " ms\n";
    cout << "scalar(double) : " << t_scalar_d << " ms\n";
    cout << "AVX2(float)    : " << t_avx2_f   << " ms  (加速 "
         << t_scalar_f / t_avx2_f << "×)\n";
    cout << "AVX2(double)   : " << t_avx2_d   << " ms  (加速 "
         << t_scalar_d / t_avx2_d << "×)\n";

    return 0;
}

编译/运行

1 2	`g++ -O3 -mavx2 -mfma -std=c++17 l2_simd_double.cpp -o l2_simd_double ./l2_simd_double`

结果分析

ref = 699061, avx2 = 699061, abs diff = 8.89073e-05, rel diff = 1.27181e-10
ref(float) = 698656, avx2(float) = 699043, abs diff = 386.812, rel diff = 0.000553652
SIMD float 结果错误！
===== L2 distance (1M floats) =====
scalar(float)  : 0.477931 ms
scalar(double) : 1.15133 ms
AVX2(float)    : 0.144433 ms  (加速 3.30902×)
AVX2(double)   : 0.130193 ms  (加速 8.84323×)

可以看到 float 版本的 SIMD 结果误差较大，达到了 0.05% 量级，而 double 版本的误差非常小，为 1e-10 量级。原因在于 float 的有效位数只有 23 位，而在累加过程中会丢失精度。同时可以看出对于标量运行的方式，使用 double 进行计算相比 float 会慢大约 2.4 倍，而使用 SIMD 指令集后，double 版本的性能提升更为显著，达到了 8.8 倍（与 AVX2 float 相近）。

HPC

#HPC

SIMD-高性能计算入门

https://github.com/Cookiecoolkid/Cookiecoolkid.github.io/2025/09/02/高性能计算入门/

作者

Cookiecoolkid

发布于

2025年9月2日

许可协议

一些学习经验和反思上一篇