引

• 从本篇开始，我们将进入一个新的板块——parallelism （并行计算）。学完之前的课程，我们可以说已经完全理解一台计算机了，现在我们要让这台计算机更加强大。
• parallel有很多类型，比如一个cpu核，可以进行多线程并行；对于一个流水线，我们可以多指令并行。然而，随着流水线的不断加深，cpu的功耗也会增加，已经达到了不可接受的程度，所以人们开始思考是否可以并行计算多个数据来提高性能。比如两个向量相加，我们可以并行执行所有元素的加法，从而提高效率。

Application: Machine Learning

• 在机器学习或者更具体的深度学习中，我们通常需要进行大量矩阵运算，比如矩阵乘法、卷积等，这些操作都可以通过并行计算来加速。
• 比如矩阵乘法，作为一种在image filter，machien learning中很常见的操作，早在C语言尚没有流行时（当时流行一种叫FORTRAN的语言），就有了专门的dgemm指令(double-precision floating-point matrix multiplication), 足以显示该操作的重要性。

Matrices

• 对于两个矩阵A和B(均为NxN)的乘法，假设结果是C矩阵，那么C矩阵中的每个元素C[i][j]都是A矩阵的第i行和B矩阵的第j列的点积，这里进行的N个乘法显然是可以同时计算的，同理，C中的每一个元素的计算也是可以同时进行的，因此矩阵乘法是一个非常适合并行计算的操作。理想情况下(无穷加法乘法器和缓存)，矩阵乘法的时间复杂度可以从$O(N^3)$降到$O(\log N)$ ?! 原因是即便有足够的处理器核心来并行计算所有东西，最后的加法还是有下界的，当形成一个树形加法结构时，时间复杂度是$O(\log N)$的。

Parallelism

Matrix Multiplication

• 对于正常的矩阵乘法：

1
2
3
4
5
6
7
8

# python version of dgemm
# It costs 2*N^3 FLOPs
def dgemm(N, a, b, c):
    for i in range(N):
        for j in range(N):
            c[i][j] = 0
            for k in range(N):
                c[i][j] += a[i][k] * b[k][j]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

// c version of dgemm
void dgemm_scalar(int N, double *a, double *b, double *c) {
    for (int i=0; i< N; i++) {
        for (int j=0; j< N; j++) {
            double cij = 0;
            for (int k=0; k< N; k++) {
                cij += a[i+k*N] * b[k+j*N];
            }
            c[i+j*N] = cij;
        }
    }
}

• 我们可以设定timer，统计两个程序的运行时间t，然后$2\cdot N^3 / t$就是FLOPS（floating-point operations per second），从而反应两个程序的效率。结果如下表所示，C语言版本竟然比Python版本快了大约240倍！

• C比python快了244倍可以理解，但是这两列数据看起来怪怪的？对C来说，当N较小时，效率就是CPU的最大频率，所以不变；当N很大时，由于缓存放不下所有的数据，所以有一个延迟，效率就会下降。对Python来说，当N不是特别大时，时间成本主要是解释器开销，时间和N^3成正比，所以效率不变？

Flynn’s Taxonomy

• Choice of hardware and software parallelism are independent
  • Concurrent software can also run on serial hardware
  • Sequential software can also run on parallel hardware
• Flynn’s taxonomy 根据指令流和数据流的不同，将计算机系统分为4类：

SIMD: specialized function units (hardware), for handling lock-step calculatiions involving arrays. 常用于Scientific computing, machine learning, signal processing …

• 下图展示了4种类型的工作方式：

其中，SIMD通过集成多个硬件单元来实现数据并行，MIMD则通过多个处理器核心来实现指令并行，可以看成是多个SIMD的组合。MISD比较少见，mainly of historical interest only.

SIMD Architectures

• 下图展示了Inter SIMD架构的进化史：

SIMD本质上是将多个数据同时加载到寄存器中，然后并行计算。

由于Inter坚持向后兼容，所以你会发现每一代新SIMD指令集都包含了之前的指令集，这样就保证了旧软件可以在新硬件上运行（即在本就很复杂的指令集上增加新指令集）。同时，随着SIMD指令集的不断发展，所需要的寄存器也不断增加和扩展。

• 对于之前矩阵乘法例子，我们使用AVX指令集来加速：

可以发现 AVX 版本的效率比普通 C 代码高了大约 4 倍，在 N 较大时同样受到缓存限制，性能下降。理论上，这颗 Intel i7-5557U 的单核浮点峰值约为 25 GFLOPS：3.1 GHz × 2 instructions/cycle × 4 mults/inst = 24.8 GFLOPS。其中 3.1 GHz 是最大睿频频率，2 instructions/cycle 是超标量发射数（这里特指两个 AVX 向量执行单元），4 mults/inst 是指一条 AVX 指令可以同时做 4 个双精度乘法。实际测得的性能（3.6–5.5 GFLOPS）远低于理论峰值，可能的原因包括：CPU 没有跑满最高频率、没有同时利用两个 AVX 单元、内存带宽受限，或者编译器生成的指令流没有完美填满流水线。但无论如何，AVX 相比普通标量 C 代码确实带来了接近 4 倍的加速，SIMD 的优势是明确的。

• 在SEE中，有专门的向量寄存器xmm0-xmm7，每一个都是128位，可以同时存储4个32位的单精度浮点数或者2个64位的双精度浮点数。在AVX中，有ymm0-ymm15，每一个都是256位。在AVX512中，有zmm0-zmm31，每一个是512位！如果你拥有一台linux电脑/服务器，或者是虚拟机，你可以使用lscpu命令在flags中查看是否支持avx512.

SIMD Array Processing

• 下面，让我们看看SIMD具体如何加速Array处理的吧！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

// pseudocode
for each f in array:
    f = sqrt(f)

// SISD
for each f in array {
    load f to the floating-point register
    calculate the square root
    write the result from the register to memory
}

// SIMD
// SEE register is 128-bit, can hold 4 single-precision floats,
// so we can process 4 elements at a time:
for every 4 members in array {
    load 4 members to the SEE register 
    calculate the square roots in one operation
    write the result from the register to memory 
}

// 回到我们之前的矩阵乘法例子，由于AVX寄存器是256位，可以同时存放4个double数据，
// 所以我们可以同时计算4个元素的乘积，这就是为什么AVX版本快了大约4倍。

• 再举一个SIMD的汇编例子：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

// C code , all data is single precision float 
vec_res.x = v1.x + v2.x;
vec_res.y = v1.y + v2.y;
vec_res.z = v1.z + v2.z;
vec_res.w = v1.w + v2.w;

// SSE assembly code
// movaps : memory aligned, packed single precision
movaps address-of-v1, %xmm0

// addps是inter风格指令，能直接对内存中的数据和寄存器中的数据进行操作
// RISC-V中是没有的，必须先load到寄存器中才能进行操作。
addps address-of-v2, %xmm0

movaps %xmm0, address-of-vec_res

• 通常来说，编译器会自动将一些简单的循环自动向量化（auto-vectorization），也就是自动将标量代码转换为SIMD代码，但对于更复杂的代码，编译器可能无法正确识别并行机会，这时我们就可以使用intrinsics或者直接手写汇编来利用SIMD指令。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

// 什么是Intrinsics? 它是一种特殊的C函数，和SIMD指令一一对应
// 使得程序员在可以在C代码中直接使用SIMD指令，而不需要写汇编代码。

// Vector data type 
    _m128d 

// Load and Store operations:
    _mm_load_pd   // MOVAPD/aligned,packed double
    _mm_store_pd  // MOVAPD/aligned,packed double 

// Arithmetic:
    _mm_add_pd   // ADDPD/add, packed double
    _mm_mul_pd    // MULPD/multiple, packed double

/* 
所以，如果你手上有一份c代码(x86-64), 你可以大量使用intrinsics来加速这个代
码, 往往可以获得更好的性能。
*/

Matrix Multiply Example

• 至此，我们知道了SIMD扩展是如何工作的了，它们是如何被指令集架构支持的，如何在硬件中实现，以及如何在汇编中使用这些指令，或者在C代码中使用内联函数。所以，让我们看一个例子，将所有东西融合，用SIMD实现一个Matrix Multiply吧！

• 如图, 假设内存是列优先，我们首先计算黄框部分。

 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

// Initialization 
C_1 [0,0]
C_2 [0,0]

// First iteration
/* 
_mm_load_pd : load 2 double values from memory, note that 
memory in Column order 
*/
A [A_11, A_21] 

/*
_mm_load1_pd : load 1 double value and stores it in the high and 
low double words of the XMM register
*/
B_1 [B_11, B_11]
B_2 [B_12, B_12]

C_1 = _mm_add_pd(C_1, _mm_mul_pd(a, b_1))
C_2 = _mm_add_pd(C_2, _mm_mul_pd(a, b_2))

// result after first iteration
C_1 [A_11*B_11, A_21*B_11]
C_2 [A_11*B_12, A_21*B_12]

// Second iteration
A [A_12, A_22]

B_1 [B_21, B_21]
B_2 [B_22, B_22]

C_1 = _mm_add_pd(C_1, _mm_mul_pd(a, b_1))
C_2 = _mm_add_pd(C_2, _mm_mul_pd(a, b_2))

// result after second iteration
C_1 [A_11*B_11 + A_12*B_21, A_21*B_11 + A_22*B_21]
C_2 [A_11*B_12 + A_12*B_22, A_21*B_12 + A_22*B_22]

/*
以上就是我们如何使用SIMD指令来计算矩阵乘法的一个例子，看似复杂，但确实充分
利用了SIMD的并行计算能力，显著提升了性能。
*/

这里提一嘴RISC-V，课程视频中（2020年）它的SIMD指令集处于草案阶段，但是由于还没有实际硬件支持，所以没有机会做RISC-V的SIMD编程😭

现在（2026），RISC-V的V系列向量指令集已经正式发布，并且有一些商用处理器已经支持！

Lab09

• 在lab09中，我们将学习如何在实际c代码中运用SIMD指令来加速，在exercise 1中，我们需要到inter官方文档中找到要求的intrinsics，熟悉它们的用法。在exercise 2中，我们将利用这些intrinsics来加速一个数组的求和，即每次都使用一个SIMD指令来同时处理4个元素，从而提高效率。在exercise 3中，我们将运用Loop Unrolling技巧一次循环做多个SIMD指令，处理更多元素，减少循环开销，进一步提升效率。

• 以下内容纯个人分析，未经验证，仅供参考🥵：

你也许会觉得exercise 2中的SIMD版本和普通版本Loop Unrolling是一个道理，即每次循环都是处理了4个元素，然而并非，SIMD版本下我们可以每次独立对4个元素进行加法，但是普通Loop Unrolling中我们每次获取4个元素后必须将它们加到一起，有一个依赖关系，所以仍是长路径。

试想在sum_unrolled函数中，初始化4个变量，sum1,sum2,sum3,sum4, 每次循环分别将4个元素加到这4个变量中，最后再将这4个变量加在一起，这样就完全消除了依赖关系，此时的sum_unrolled_optimal确实和SIMD没有逻辑上的区别。然而，这只是理论上没有区别，实际上由于硬件执行层面是一条一条指令发射的，所以4个加法仍然有先后顺序，无法完全并行执行，所以性能上不会有很大提升（有可能是add比较快，依赖延迟小）。

而SIMD真正做到了发射一条指令，同时加载4个元素，同时做4个加法，真正实现了并行计算，所以性能提升是非常显著的。如果再用上Loop Unrolling的话，一次循环发射多条SIMD指令，性能就更好了。此时，再初始化4个128位的寄存器，分为temp1,temp2,temp3,temp4，消除依赖关系，竟然有0.2s的提升！可能是128位SIMD加法的延迟大于普通加法？所以在向量加法器充足情况下，消除依赖关系的效果就更明显了。

其实这个就是data hazards, 只不过普通add延迟小，所以无法提现消除data hazards的效果，而SIMD加法延迟大，所以消除data hazards的效果就很明显了。