SSE 真的有那么神奇吗?

发表于

最近在知乎上看到一个问题:  为什么 C++ 只比 VBA 快 4倍 里面题主提到了一个简单的估算PI的程序, 发现运行时候C++ 需要150ms, VBA需要570ms, 发现只有4倍的差异.

其中一个miloyip 的回答使用了SSE加速计算, 结果最后估计比没有使用SSE版本快了6,7倍左右.

这个速度差异令我感到非常吃惊, 回想到当年初高中时期看电脑报\微型计算机等杂志时候都会提到新发布处理器带来新的指令集等等变化, 一直没能力测试, 于是现在打算写程序测试一下速度差异到底有多大.

完整代码还是最后贴吧, 先简单介绍一下测试方式:

计算方式是首先生成 N 个随机浮点数(float), 再对他们进行开平方运算, 看看所需的时间.

C++自带的随机数生成函数还是比较慢的, 但是由于这次测试重点不在这一方面, 因此计算时间时候把这一部分去掉了.

时间计算使用了C++11提供的新的 system_clock 类, 可以提供极高精度的时间, 我将其换算到ms.

首先是默认计算的函数:

1
2
3
4
void normalSqrt(float *a, int N) {
    for (int i = 0; i != N; ++i)
        a[i] = sqrt(a[i]);
}

然后是使用SSE指令集的计算, SSE提供了128位寄存器, 可以一次保存4个float计算.

1
2
3
4
5
void SSESqrt(float *a, int N) {
    __m128* ptr = (__m128*) a;
    for (int i = 0; i != N / 4; ++i, ++ptr, a += 4)
        _mm_store_ps(a, _mm_sqrt_ps(*ptr));
}

得出的结果跟源答案的差距差不多, 对于6400W个浮点数:

普通计算方式需要大约491ms, 而SSE只需要71ms.差距非常明显.

写到这里基本上应该已经结束了, 但是我突然脑洞一开, 回想到当年OI时期, 有一个基本的优化思路, 就是对于大的数组循环, 每次i+4, 然后在循环内部展开4个操作, 据说这样可以更快.于是又写了两个函数测试分别在循环展开4个和2个的情况下的速度:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void OISqrt(float *a, int N) {
    for (int i = 0; i < N; i += 4) {
        a[i] = sqrt(a[i]);
        a[i+1] = sqrt(a[i+1]);
        a[i+2] = sqrt(a[i+2]);
        a[i+3] = sqrt(a[i+3]);
    }
}

void OISqrt2(float *a, int N) {
    for (int i = 0; i < N; i += 2) {
        a[i] = sqrt(a[i]);
        a[i+1] = sqrt(a[i+1]);
    }
}

测试结果, 对于展开4层的循环, 速度是438ms, 2层的循环, 需要457ms. 这个结果说不上多好, 距离当年所期待的将近4倍\2倍来说还是差很多, 但是也有50ms+的提升, 这对于计算密集的应用来说已经算是很不错了. 特别是如果限时1秒的话, 这可是1/20 的提升.

写到这里我脑洞又一大开, 回想起来以前看到Intel从Core I 第二代开始加入的AVX指令集, 又看了看, 似乎提供了256bit的操作宽度, 这样一来岂不是意味着一次可以同时处理8个float了. 于是又查了查资料, 写了一个利用AVX指令集的函数:

1
2
3
4
5
void AVXSqrt(float *a, int N) {
    for (int i = 0; i != N / 8; ++i, a += 8) {
        _mm256_store_ps(a, _mm256_sqrt_ps(_mm256_load_ps(a)));
    }
}

这次的运行结果是63ms, 相对于SSE只有非常小的提升, 我运行了很多次, 结果都是有很小提升. 网上有人提到对于Sandy Bridge 来说, AVX的sqrt是使用SSE模拟出来的.

 

写到这里就是真的结束了, 其实本来还打算测试一下AVX2 指令集, 但是这个只支持到Haswell 且添加的内容似乎跟这个测试没什么关系, 而AVX-512现在还没有支持的处理器.

于是最终对比结果是:

1
2
3
4
5
6
Rand generate:667
Normal:491
SSE:71
OI4:438
OI2:457
AVX:63

附全部代码, 随手写的, 欢迎吐槽:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include <immintrin.h>
#include <emmintrin.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <chrono>
#include <cmath>
#include <iostream>

using namespace std;
using namespace std::chrono;


inline milliseconds getCurrentTime() {
    return duration_cast<milliseconds>(chrono::system_clock::now().time_since_epoch());
}

inline unsigned long getDiffTime(milliseconds &t) {
    return (getCurrentTime() - t).count();
}

inline void echoDiff (string name, unsigned long t) {
    cout<<name<<":"<<t<<endl;
}

void normalSqrt(float *a, int N) {
    for (int i = 0; i != N; ++i)
        a[i] = sqrt(a[i]);
}

void SSESqrt(float *a, int N) {
    __m128* ptr = (__m128*) a;
    for (int i = 0; i != N / 4; ++i, ++ptr, a += 4)
        _mm_store_ps(a, _mm_sqrt_ps(*ptr));
}

void OISqrt(float *a, int N) {
    for (int i = 0; i < N; i += 4) {
        a[i] = sqrt(a[i]);
        a[i+1] = sqrt(a[i+1]);
        a[i+2] = sqrt(a[i+2]);
        a[i+3] = sqrt(a[i+3]);
    }
}

void OISqrt2(float *a, int N) {
    for (int i = 0; i < N; i += 2) {
        a[i] = sqrt(a[i]);
        a[i+1] = sqrt(a[i+1]);
    }
}

void AVXSqrt(float *a, int N) {
    for (int i = 0; i != N / 8; ++i, a += 8) {
        _mm256_store_ps(a, _mm256_sqrt_ps(_mm256_load_ps(a)));
    }
}
int main(int argv, char **argc) {
    srand(time(NULL));
    int N = 64000000;
    float *a;
    posix_memalign((void**)&a, 32,  N * sizeof(float));

    auto t = getCurrentTime();
    for (int i = 0; i != N; ++i)
        a[i] = rand();
    echoDiff("Rand generate", getDiffTime(t));

    t = getCurrentTime();
    normalSqrt(a, N);
    echoDiff("Normal", getDiffTime(t));

    t = getCurrentTime();
    SSESqrt(a, N);
    echoDiff("SSE", getDiffTime(t));

    t = getCurrentTime();
    OISqrt(a, N);
    echoDiff("OI4", getDiffTime(t));
   
    t = getCurrentTime();
    OISqrt2(a, N);
    echoDiff("OI2", getDiffTime(t));

    t = getCurrentTime();
    AVXSqrt(a, N);
    echoDiff("AVX", getDiffTime(t));

    return 0;
}