SIMD

合集 - 并行(7)

1.nku并行01-132.并行硬件和并行软件01-13

3.SIMD01-13

4.Pthread01-135.OpenMP01-136.MPI01-137.CUDA01-13

收起

SIMD概念#

当前的SIMD架构#

多媒体扩展：SSE、AVX
图形和游戏处理器：CUDA

SIMD并行的问题#

我的理解：

• 通常在向量寄存器上执行，这些寄存器比普通的CPU寄存器宽，可以存储多个数据元素。例如，一个128位的SIMD寄存器可以存储四个32位的浮点数。
• 数据被组织成向量，每个向量包含多个数据元素。CPU执行一条指令时，这条指令会应用于向量中的所有元素。（向量化）
• 为了最大化SIMD的性能，数据通常需要以连续块的形式存储在内存中，这样可以一次性加载到SIMD寄存器中。

额外开销#

打包/解包数据的开销#

打包源运算对象：拷贝到连续内存区域
解包目的运算对象：拷贝回内存

对齐#

• 最坏情况需要计算地址，动态对齐
• 编译器（程序员）可分析确认对齐
  • 一般而言数据是从起始地址处对齐的
  • 如果在一个循环中顺序访问数据，起始位置固定，则对齐特性是不变的
• 未对齐的，可调整算法，先串行处理到对齐边界，然后进行SIMD计算

有时对齐开销会完全抵消SIMD的并行收益

控制流#

如条件分支、循环、同步等

如果有控制流，要求执行所有路径

for (i=0;i<16;i++)
	if (a[i]!=0)
		b[i]++;

对所有元素计算$a[i]!=0$
计算所有元素的$b[i]++$存入临时寄存器$t1$
将原$b[i]$拷贝到寄存器$t2$
根据a$[i]!=0$的结果合并$t1$和$t2$

一般观点，当存在控制流问题时，SIMD不是一个好的编程模型

SSE/AVX编程#

概念#

SSE（Streaming SIMD Extensions），是由英特尔在其 x86 处理器中引入的一组 SIMD（单指令多数据） 指令集扩展，首次出现在 1999 年的 Pentium III 处理器中。
8 个 128 位宽的寄存器

AVX(Advanced Vector eXtensions)，2011年
16 个 256 位宽的寄存器

intrinsics内部函数#

支持多种数据类型

• 整数(16字节、8 short、4int、2 longlong、1dqword)
• 单精度浮点数(4 float)
• 双精度浮点数(2 double)

可以计算，比如int是32bit的，$128/32=4$，所以可以有4个int
(以上默认都是SSE的)

常用指令#

在C++中，intrinsic（内部函数）是一种特殊的函数，它提供了对底层硬件指令的直接访问

头文件：#include <immintrin.h>
这个头文件包含了SSE、AVX等指令集的内部函数。

数据类型命名规律#

三部分

1. __m，两个_
2. 位宽
3. 数据类型—— i：整数，d：双精度浮点数，空：单精度浮点数
   例子：
   float : __m128 , int：__m128i, double：__m128d

函数命名规律#

三部分

1. _mm：SSE，_mm256：AVX，_mm512：AVX-512
2. 操作——_add：加法， _mul：乘法，_load：加载到寄存器中，_loadu 将 16 位未对齐的操作数加载到寄存器中。
3. 数据类型——ps：(四个）单精度浮点数，pd：（两个）双精度浮点数，pixx：xx 为长度，有符号整数，64 位寄存器，epixx： 有符号整数，128 位寄存器，epuxx 无符号整数，ss 操作第一个（低位）单精度浮点数（标量）。
   (加u才是不对齐的)

如_mm_mul_epi32 对参数中所有的 32 位有符号整数进行乘法运算。

矩阵乘法#

重点

串行版本#

void mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]) {
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            c[i][j] = 0.0;
            for (int k = 0; k < n; ++k) {
                c[i][j] += a[i][k] * b[k][j];//
            }
        }
    }
}

cache优化#

转置矩阵b，让a,b都是行优先，提高缓存命中率

void mul2(int n, float a[][maxN], float b[][maxN], float c[][maxN]) {
    for (int i = 0; i < n; i++)
        for (int j = 0; j < i; j++)
            swap(b[i][j], b[j][i]);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            c[i][j] = 0.0;
            for (int k = 0; k < n; ++k) {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
}

SSE优化#

_mm_hadd_ps水平加法，对于a(A3,A2,A1,A0)和b(B3,B2,B1,B0)，`_mm_hadd_ps(a,b)=（A3+A2,A1+A0,B3+B2,B1+B0)

void mul3(int n, float a[][maxN], float b[][maxN], float c[][maxN]) {
    for (int i = 0; i < n; i++)
        for (int j = 0; j < i; j++)
            swap(b[i][j], b[j][i]);
    __m128 sum, t1, t2;//
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {//n
            c[i][j] = 0.0;
            sum = _mm_setzero_ps();//
            for (int k = n - 4; k >= 0; k -= 4) {
                //注意是loadu而不是load
                t1 = _mm_loadu_ps(a[i] + k);//a[i] + k=&a[i][k]
                t2 = _mm_loadu_ps(b[j] + k);
                t1 = _mm_mul_ps(t1, t2);
                sum = _mm_add_ps(sum, t1);
            }
            sum = _mm_hadd_ps(sum, sum);
            sum = _mm_hadd_ps(sum, sum);
            _mm_store_ss(c[i] + j, sum);//只存一个标量，没有u
            //处理剩下的
            for (int k = (n % 4) - 1; k >= 0; k--) {
                c[i][j] += a[i][k] * b[j][k];
            }
        }
    }
}

分片策略#

SSE的问题是，A的每个元素进入cache一次，但B要进N次（先A矩阵某一行和B矩阵每一行计算，B每次取的都不同，所以要N次进入cache，命中率低），对整个矩阵进$N+N\ast N$次
优化方案是把矩阵分成几个子矩阵，对于宽为$T$的子矩阵，A、B每个元素都要进入$N/T$次，那么对整个矩阵只用进$T\ast T\ast 2$次

void mul4(int n, float a[][maxN], float b[][maxN], float c[][maxN]) {
    __m128 t1, t2, sum;
    float t;
    for (int i = 0; i < n; ++i) for (int j = 0; j < i; ++j) swap(b[i][j], b[j][i]);
    for (int r = 0; r < n / T; ++r)
        for (int q = 0; q < n / T; ++q) {
            //分块，以下是对每一个块进行计算,r和c分别表示a,b块的索引
            for (int i = 0; i < T; ++i)
                for (int j = 0; j < T; ++j) 
                    c[r * T + i][q * T + j] = 0.0;
            for (int p = 0; p < n / T; ++p) {
                //p 表示矩阵 k 的块索引,
                for (int i = 0; i < T; ++i)
                    for (int j = 0; j < T; ++j) {
                        sum = _mm_setzero_ps();
                        for (int k = 0; k < T; k += 4){
                            t1 = _mm_loadu_ps(a[r * T + i] + p * T + k);
                            t2 = _mm_loadu_ps(b[q * T + j] + p * T + k);
                            t1 = _mm_mul_ps(t1, t2);
                            sum = _mm_add_ps(sum, t1);
                        }
                        sum = _mm_hadd_ps(sum, sum);
                        sum = _mm_hadd_ps(sum, sum);
                        _mm_store_ss(&t, sum);
                        c[r * T + i][q * T + j] += t;
                    }
            }
        }
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < i; ++j) swap(b[i][j], b[j][i]);
}

AVX#

AVX就是加宽了的SSE，一次处理8个float

void avx_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]) {
    __m256 t1, t2, sum;
    _m128 s1,s2;
    for (int i = 0; i < n; ++i) for (int j = 0; j < i; ++j) swap(b[i][j], b[j][i]);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            sum = _mm256_setzero_ps();
            for (int k = 0; k < n; k += 8) {
                t1 = _mm256_load_ps(&a[i][k]);
                t2 = _mm256_load_ps(&b[k][j]);
                sum = _mm256_add_ps(sum, _mm256_mul_ps(t1, t2));
            }
            //从一个 256-bit 的 __m256 类型向量中提取出低128-bit 的 __m128 向量,前四个
            s1 = _mm256_extractf128_ps(sum, 0);
            //从一个 256-bit 的 __m256 类型向量中提取出高128-bit 的 __m128 向量,后四个
            s2 = _mm256_extractf128_ps(sum, 1);
            //将原来的8个合成4个
            s1 = _mm_hadd_ps(s1, s2);
            //4个自己水平加法两次合成1个
            s1 = _mm_hadd_ps(s1, s1);
            s1 = _mm_hadd_ps(s1, s1);
            _mm_store_ss(&c[i][j], s1);//取第一个数存入c[i][j]
    }
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < i; ++j) swap(b[i][j], b[j][i]);
}

其他指令：

• shuffle洗牌

__m128 a = _mm_set_ps(1.0f, 2.0f, 3.0f, 4.0f); 
__m128 b = _mm_set_ps(5.0f, 6.0f, 7.0f, 8.0f); // 使用 _mm_shuffle_ps 将 a 和 b 中的元素重新排列 
__m128 result = _mm_shuffle_ps(a, b, _MM_SHUFFLE(2, 1, 0, 3));

• blend选择性合并

__m128 a = _mm_set_ps(1.0f, 2.0f, 3.0f, 4.0f); 
b = _mm_set_ps(5.0f, 6.0f, 7.0f, 8.0f); 
// 使用掩码来选择性地将 a 和 b 的元素合并 
__m128 result = _mm_blend_ps(a, b, 0b1100); // 掩码：0b1100 -> 选择 a 的前两个元素，b 的后两个元素