optimizing_cpp(10)

十、向量操作

现代处理器基本都支持向量指令了，这种技术叫SIMD(single instruction multiple data)。向量大小从64位（MMX）、128位(XMM)、256位(YMM)、512位(ZMM)

向量操作适用于对多个数据执行相同的操作且运行并行，比如图像处理、声音处理、矩阵运算。

下面附上各个向量指令集支持的数据类型

Type of elements	Size of each elements, bits	Number of elements	Total size of vector, bits	Instruction set
`char`	8	8	64	MMX
`short int`	16	4	64	MMX
`int`	32	2	64	MMX
`int64_t`	64	1	64	MMX

`char`	8	16	128	SSE2
`short int`	16	8	128	SSE2
`int`	32	4	128	SSE2
`int64_t`	64	2	128	SSE2
`float`	32	4	128	SSE
`double`	64	2	128	SSE2

`char`	8	32	256	AVX2
`short int`	16	16	256	AVX2
`int`	32	8	256	AVX2
`int64_t`	64	4	256	AVX2
`float`	32	8	256	AVX
`double`	64	4	256	AVX

`char`	8	64	512	AVX512BW
`short int`	16	32	512	AVX512BW
`int`	32	16	512	AVX512
`int64_t`	64	8	512	AVX512
`float`	32	16	512	AVX512
`double`	64	8	512	AVX512

这里忽略了anger关于avx256的介绍，因为查阅资料后发现有不少错误，比如他强调的_mm256_zeroupper() ，这个函数会由编译器自动加上，程序员不应当手动调用。

另外附上各阶段指令集支持类型

MMX	64 位整型
SSE	128 位浮点运算，整数运算仍然要使用 MMX 寄存器，只支持单精度浮点运算
SSE2	对整型数据的支持，支持双精度浮点数运算，CPU 快取的控制指令
SSE3	扩展的指令包含寄存器的局部位之间的运算，例如高位和低位之间的加减运算；浮点数到整数的转换，以及对超线程技术的支持。
SSE4
AVX	256 位浮点运算
AVX2	对 256 位整型数据的支持，三运算指令（3-Operand Instructions）
AVX512	512 位运算

另外发现一个挺不错的入门介绍，贴上来一文读懂SIMD指令集目前最全SSE/AVX介绍_Axurq的博客-CSDN博客_simd指令集

1.自动向量化

许多编译器有自动向量化的功能，比如gcc中-O3就开启了自动向量化，除此之外，还可以通过-fopt-info-vec-optimized或-fopt-info-vec命令来查看是否开启向量化。

自动向量化最重要的方面就是对齐，尽量保证在内存中是地址对齐的。

内存对齐：

内存对齐是一个常见的现象，网上查阅资料其实说的都不是很详尽，个人整理理解如下

为什么要内存对齐？

1.现代计算机中内存空间都是按照byte划分的，从理论上讲似乎对任何类型的变量的访问可以从任何地址开始，但实际情况是在访问特定变量的时候经常在特定的内存地址访问，比如2、4、8的倍数。

2.由于1的原因，如果内存没有对齐在读内存时可能本来只需要读一次，但却需要读两次。就像下面这张图（太丑了凑合看），如果不内存对齐，读double时先读前八个字节，再读后八个字节，需要读两次，但是右边的图内存对齐后只需要读一次就行了。

并且经过测试，编译器除了会自动内存对齐，还会将内存地址也对齐，其元素的内部地址总是能被结构体中最大的元素的大小整除，方便cpu进行内存访问。

仍然存疑关于avx256要32位对齐：64位cpu一次访存读取的块是多大的？是64位吗，那为什么向量化从内存载入时要求还32字节对齐呢（8字节对齐不就行了吗），cpu可以一次读取32字节数据吗？

2.手动向量化

一般来说尽量让编译器自动向量化，但是有些情况编译器无法自动向量化，而向量化带来的收益很可观就可以手动向量化，例如下面这个栗子，当存在分支时，编译器很难向量化，此时可以利用内部函数手动向量化。

// Example 12.4a. Loop with branch

// Loop with branch
void SelectAddMul(short int aa[], short int bb[], short int cc[])
{
    for (int i = 0; i < 256; i++)
    {
        aa[i] = (bb[i] > 0) ? (cc[i] + 2) : (bb[i] * cc[i]);
    }
}

可以改成

// Example 12.4b. Vectorized with SSE2

#include <emmintrin.h> // Define SSE2 intrinsic functions
// Function to load unaligned integer vector from array
static inline __m128i LoadVector(void const * p)
{
    return _mm_loadu_si128((__m128i const*)p);
}
// Function to store unaligned integer vector into array
static inline void StoreVector(void * d, __m128i const & x)
{
    _mm_storeu_si128((__m128i *)d, x);
}
// Branch/loop function vectorized:
void SelectAddMul(short int aa[], short int bb[], short int cc[])、
{
    // Make a vector of (0,0,0,0,0,0,0,0)
    __m128i zero = _mm_set1_epi16(0);
    // Make a vector of (2,2,2,2,2,2,2,2)
    __m128i two = _mm_set1_epi16(2);
    // Roll out loop by eight to fit the eight-element vectors:
    for (int i = 0; i < 256; i += 8)
    {
        // Load eight consecutive elements from bb into vector b:
        __m128i b = LoadVector(bb + i);
        // Load eight consecutive elements from cc into vector c:
        __m128i c = LoadVector(cc + i);
        // Add 2 to each element in vector c
        __m128i c2 = _mm_add_epi16(c, two);
        // Multiply b and c
        __m128i bc = _mm_mullo_epi16 (b, c);
        // Compare each element in b to 0 and generate a bit-mask:
        __m128i mask = _mm_cmpgt_epi16(b, zero);
        // AND each element in vector c2 with the bit-mask:
        c2 = _mm_and_si128(c2, mask);
        // AND each element in vector bc with the inverted bit-mask:
        bc = _mm_andnot_si128(mask, bc);
        // OR the results of the two AND operations:
        __m128i a = _mm_or_si128(c2, bc);
        // Store the result vector in eight consecutive elements in aa:
        StoreVector(aa + i, a);
    }
}

手动向量化如何对齐数据

1
2
3

// Example 12.5. Aligned arrays
const int size = 256; // Array size
alignas(16) int16_t aa[size]; // Make aligned array

查找表向量化

查找表是一种很高效的优化手段，以空间换时间，后面会详细讲解如何使用。查找表通常是向量化的障碍，从avx2开始支持对查找表有利的函数如gather等。

3.将串行代码向量化

有许多代码执行相同的操作，但是却因为循环依赖无法向量化和并行，我们可以通过修改代码结构来使其能够向量化，例如：

// Example 12.8a. Sum of a list
float a[100];
float sum = 0;
for (int i = 0; i < 100; i++) 
    sum += a[i];

上述的代码是串行的，因为每次迭代 sum 的值都依赖于前一次迭代后 sum 的值。诀窍是将循环按 n 展开并重新组织代码，每个值依赖于 n 个位置之前的值，其中 n 是向量中元素的数量。如果 n = 4，我们得到：

// Example 12.8b. Sum of a list, rolled out by 4
float a[100];
float s0 = 0, s1 = 0, s2 = 0, s3 = 0, sum;
for (int i = 0; i < 100; i += 4)
{
    s0 += a[i];
    s1 += a[i+1];
    s2 += a[i+2];
    s3 += a[i+3];
}
sum = (s0+s1)+(s2+s3);

现在，s0、s1、s2 和 s3 可以组合成一个128位的向量，这样我们就可以在一个操作中做4个加法。如果我们使用 fast math 选项并指定SSE 或更高指令集的选项，一个好的编译器会自动将例 12.8a转换为12.8b，并将代码向量化。

再一些更复杂的情况下不能自动向量化。例如，让我们看看泰勒级数的例子。指数函数可由级数计算：

$$
e^x=\sum_{n=0}^\infty\frac{x^n}{n!}
$$
用C++ 实现看起来可能是这样的：

// Example 12.9a. Taylor series
float Exp(float x)
{
    // Approximate exp(x) for small x
    float xn = x; // x^n
    float sum = 1.f; // sum, initialize to x^0/0!
    float nfac = 1.f; // n factorial
    for (int n = 1; n <= 16; n++)
    {
        sum += xn / nfac;
        xn *= x;
        nfac *= n+1;
    }
return sum;
}

在这里每个 x^n 的值由前一个值计算而来，即 x^n = x*(x^{n-1})，每个 n! 的值也由前一个值计算而来，即 n!= n*(n-1)!。如果我们想要将循环按 4展开，那我们必须要用 4个位置之前的值来计算当前的值。因此，我们将用 x^4*x^{n-4}来计算x^n。没有简单的方法来展开阶乘的计算，但是这个并不是必需的，因为阶乘并不依赖 x ，我们可以将值预先计算好，存一个表中。更好的方法是存储阶乘的倒数，这样我们就不需要除法了（如你所知，除法是很慢的）。现在上述的的代码可以按如下的方式向量化（使用 Intel vector classes）：

// Example 12.9b. Taylor series, vectorized
#include <dvec.h> // Define vector classes (Intel)
#include <pmmintrin.h> // SSE3 required
// This function adds the elements of a vector, uses SSE3.
// (This is faster than the function add_horizontal)
static inline float add_elements(__m128 const & x)
{
    __m128 s;
    s = _mm_hadd_ps(x, x);
    s = _mm_hadd_ps(s, s);
    return _mm_cvtss_f32(s);
}
float Exp(float x)
{
    // Approximate exp(x) for small x
    __declspec(align(16)) // align table by 16
    const float coef[16] = { // table of 1/n!
    1., 1./2., 1./6., 1./24., 1./120., 1./720., 1./5040.,
    1./40320., 1./362880., 1./3628800., 1./39916800.,
    1./4.790016E8, 1./6.22702E9, 1./8.71782E10,
    1./1.30767E12, 1./2.09227E13};
    float x2 = x * x; // x^2
    float x4 = x2 * x2; // x^4
    // Define vectors of four floats
    F32vec4 xxn(x4, x2*x, x2, x); // x^1, x^2, x^3, x^4
    F32vec4 xx4(x4); // x^4
    F32vec4 s(0.f, 0.f, 0.f, 1.f); // initialize sum
    for (int i = 0; i < 16; i += 4)
    {
        // Loop by 4
        s += xxn * _mm_load_ps(coef+i); // s += x^n/n!
        xxn *= xx4; // next four x^n
    }
    return add_elements(s); // add the four sums
}

这个循环在一个向量中计算四个连续的项。如果循环很长，那么进一步展开循环可能是值得的，因为这里的速度可能受到 xxn 相乘的延迟而不是吞吐量的限制（参见11 乱序执行）。这里的系数表是在编译时计算的。在运行时计算系数表可能更方便，只要你能确保系数表只被计算一次，而不是每次调用函数时都会被计算一次。

4.向量运算的数学函数

除了英特尔官方的向量基本运算函数外，还有许多函数库支持对数、指数、三角等数学函数的运算。

有两种向量数学库，长向量库和短向量库，长向量库在每一步等待所有值计算完成后再进入下一步，这意味着在内存中会有中间数组存储中间值，而短向量库对一组数据（比如avx256就是4个double）执行全部步骤后得到最终结果再算下一组数据，这样做就不会产生中间值。一般来说短向量库要高效一些因为没有中间值，但是当数学函数的计算存在依赖链时，短向量库由于阻止了乱序执行效率会显著降低。

下面是一些长向量数学库的列表：

Intel vector math library (VML, MKL)。支持所有x86 平台。这个库降低了非英特尔CPU 的性能，除非你重写了英特尔的 CPU分派程序。见13.7 Intel 编译器中的 CPU分派。
Intel Performance Primitives (IPP)。支持所有x86 平台。在非英特尔 CPU 上也能工作的很好。包括许多用于统计学，信号处理和图像处理的函数。
Yeppp。开源库。支持 x86 和 ARM 平台以及多种编程语言。www.yeppp.info。

下面是一些短向量数学库的列表：

Sleef library。支持多种不同的平台。开源。www.sleef.org。
Intel short vector math library (SVML)。这是由 Intel 编译器 提供的，并通过自动向量化调用。Gnu编译器 可以通过选项 -mveclibabi=svml 使用这个库。如果不使用 Intel 编译器话，这个库通常可以很好地处理 非Intel CPU。见13.7 Intel 编译器中的 CPU分派。
AMD LIBM library。只在64位Linux 和Windows 平台上可用。这个库在没有FMA4 指令集的情况下降低了CPU 的性能（这个指令集最初是由英特尔设计的，但目前只有 AMD 的CPU 支持）。在 Gnu编译器 中，可以通过选项 -mveclibabi=acml 使用这个库。
VCL vector class library。支持所有x86 平台。支持 Microsoft、Intel、Gnu 和Clang 编译器。代码是内联的，不需要链接外部库。www.agner.org/optimize/#vectorclass。

所有这些库都具有很好的性能和精度。速度比任何非向量库快很多倍。

SVML 和 LIBM 库中的函数名没有很好的说明文档。如果你想直接调用库函数，可以参考下表中的例子：

Library	exp function of 4 floats	exp function of 2 double
Intel SVML v.10.2 & earlier	`vmlsExp4`	`vmldExp2`
Intel SVML v.10.3 & later	`__svml_expf4`	`__svml_exp2`
Intel SVML + ia32intrin.h	`_mm_exp_ps`	`_mm_exp_pd`
AMD Core Math Library	`__vrs4_expf`	`__vrd2_exp`
AMD LIBM Library	`amd_vrs4_expf`	`amd_vrd2_exp`
VCL vector class library	`exp`	`exp`