高性能编程基础硬件知识

jokervTv 2021/12/15 HPC

本教程的目的是向非专业程序员简要概述现代硬件的功能，您必须了解这些功能才能编写高性能代码。翻译自 jakobnissen - hardware_introduction (opens new window)

这个笔记不是什么
计算机硬件的基本结构
避免过于频繁地访问磁盘
避免缓存未命中
保持数据与内存对齐
题外话：汇编代码
- 快指令，慢指令
分配和不变性
- 并非所有对象都需要分配
寄存器和 SIMD
- SIMD 需要固定长度的不间断迭代
- SIMD 需要一个循环，其中循环顺序无关紧要
数组结构
专用 CPU 指令
- 调用任何 CPU 指令
内联
展开
避免不可预测的分支
注意内存依赖
可变时钟速度
多线程
- 可并行性
GPUs

当今许多科学领域都在使用编程，在这些领域中，个别科学家通常必须为自己的项目编写自定义代码。然而，对于大多数科学家来说，计算机科学并不是他们的专业领域。他们必须学习编程。我把自己算作他们中的一员。虽然我们可能对编程的软件方面相当熟悉，但我们很少对计算机硬件如何影响代码性能有基本的了解。

这将是对过去几年所学知识的提炼。本教程将使用 Julia，因为它允许使用高级交互语言轻松演示这些相对低级的考虑。

# 这个笔记不是什么

# 这不是 Julia 编程语言的指南

要编写快速代码，您必须首先了解您的编程语言及其特性。但这不是 Julia 编程语言的指南。我建议阅读 Julia 文档的性能提示 (opens new window)部分。

# 这不是对特定数据结构或算法的解释

除了了解您的语言之外，您还必须了解自己的代码才能使其快速运行。您必须了解 $\mathbf{O}$ 符号背后的思想、为什么某些算法比其他算法更快，以及不同的数据结构如何在内部工作。在不知道数组是什么的情况下，您怎么可能优化使用数组的代码？

这也超出了本文的范围。但是，我认为程序员至少应该了解：

二进制整数如何在内存中表示
浮点数在内存中是如何表示的（学习这一点对于理解浮点运算的计算误差也是必要的，这在进行科学编程时是必须的）
包含 ASCII 和 UTF-8 编码的字符串的内存布局
数组结构的基础知识，以及密集数组之间的区别，例如整数和对对象的引用数组
Dict (即哈希表) 和 Set 工作原理背后的原理

此外，我还建议您熟悉：

Heaps (堆)
Deques (双端队列)
Tuples (元组)

# 这不是对代码进行基准测试的教程

为了在实践中编写快速代码，有必要分析您的代码以找出您的机器花费大部分时间的瓶颈。必须对不同的功能和方法进行基准测试，才能在实践中找到最快的。 Julia（和其他语言）有专门用于此目的的工具，但我不会在这里介绍它们。

# 计算机硬件的基本结构

现在，我们将使用简化的计算机模型。通过本文档，我将在我们的模型变得具体时添加更多细节。

$[CPU] \leftrightarrow [RAM] \leftrightarrow [DISK]$

在这个简单的“图表”中，箭头代表任一方向的数据流。该图显示了计算机的三个重要部分：

中央处理器 (CPU) 是一枚邮票大小的芯片。这是所有计算实际发生的地方，即计算机的大脑。
随机存取存储器（RAM，或简称“内存”）是计算机的短期记忆。内存需要电力来维持，并在计算机关闭时丢失。RAM 用作磁盘和 CPU 之间数据的临时存储。“加载”各种应用程序和操作系统的大部分时间实际上都花在将数据从磁盘移动到 RAM 并在那里解压。典型的消费类笔记本电脑有大约 $10^{11}$ 位 RAM 内存。
磁盘是大容量存储单元。关闭电源后，磁盘上的这些数据仍然存在，因此磁盘包含计算机的长期记忆。它的每 GB 也比 RAM 便宜得多，消费类 PC 具有大约 $10^{13}$ 位的磁盘空间。

# 避免过于频繁地访问磁盘

在讨论软件性能时，区分吞吐量和延迟很有用。延迟是从某事开始到完成所花费的时间。吞吐量是衡量在一定时间内完成多少工作的指标。

从表面上看，延迟和吞吐量之间的关系似乎很明显：如果一个操作需要 $N$ 秒钟来计算，那么每秒可以完成 $\dfrac{1}{N}$ 个操作。所以你会认为：

$朴素的方程：吞吐量 = \dfrac{1}{延时}$

实际上，事情并没有那么简单。例如，假设一个操作在开始前有 1 秒的“预热”，但之后在 0.1 秒内完成。因此延迟为 1.1 秒，但初始预热后的吞吐量为 10 操作数/秒。

或者，想象一个延迟为 1 秒的操作的情况，但可以同时运行 8 个操作。这些操作可以批量运行，吞吐量为 8 操作/秒。

区分延迟和吞吐量很有用的一个地方是程序从磁盘读取时。大多数现代计算机使用一种称为固态驱动器 (SSD) 的磁盘。当前 (2021) SSD 的延迟约为 100 µs，读/写吞吐量远远超过 1 GB/s。(取决于型号，低端消费级 SSD 并不能达到) 较旧或更便宜的大容量存储磁盘属于硬盘驱动器 (HDD) 类型。它们的延迟增加了 100 倍，接近 10 毫秒，吞吐量只有约 100 MB/s。

即使是最新、最快的 SSD，其延迟也比 RAM 慢数千倍，后者的延迟低于 100 纳秒。每当发生读取或写入时都会产生磁盘延迟。因此，要编写快速代码，您必须不惜一切代价避免重复读取或写入磁盘。

以下示例用于说明延迟的差异：第一个函数打开一个文件，访问文件中的一个字节，然后再次关闭它。第二个函数从 RAM 中随机访问 1,000,000 个整数。

begin
  # Open a file
  function test_file(path)
    open(path) do file
      # Go to 1000'th byte of file and read it
      seek(file, 1000)
      read(file, UInt8)
    end
  end

  # Randomly access data N times
  function random_access(data::Vector{UInt}, N::Integer)
    n = rand(UInt)
    mask = length(data) - 1
    @inbounds for i in 1:N
      n = (n >>> 7) ⊻ data[n & mask + 1]
    end
    return n
  end
end;

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

let
  data = rand(UInt, 2^24)
  with_terminal() do
    @time test_file("../alen/LICENSE")
    @time random_access(data, 1000000);
  end
end

1
2
3
4
5
6
7

  0.002043 seconds (13 allocations: 816 bytes)
  0.119351 seconds

1
2

对此进行基准测试有点棘手，因为第一次调用将包括两个函数的编译时间。在第二次调用中，您的操作系统将在 RAM 中存储文件的副本（或缓存文件），使文件搜索几乎立即进行。要正确计时，请运行一次，然后将文件更改为另一个最近未打开的文件，然后再次运行。所以实际上，我们应该更新我们的计算机图：

$[CPU] \leftrightarrow [RAM] \leftrightarrow [DISKCACHE] \leftrightarrow [DISK]$

在我的计算机上，查找文件中的单个字节（包括打开和关闭文件）大约需要 500 微秒，而从内存中访问 1,000,000 个整数需要 200 毫秒。因此 RAM 延迟大约比磁盘低 2,500 倍。因此，在高性能计算中必须避免重复访问文件。

仅在几年前，SSD 还不常见，HDD 吞吐量还低于今天。因此，旧文本通常会警告人们不要让您的程序完全依赖磁盘以获得高吞吐量。这个建议在今天已经过时了，因为大多数程序甚至无法在 1 GB/s 的廉价现代 SSD 的吞吐量上形成瓶颈。今天的建议仍然适用于需要频繁对磁盘进行单独读/写的程序，因为在那里会积累高延迟。在这些情况下，您确实应该将数据保存在 RAM 中。

性能最糟糕的情况是，如果您需要以小块读取/写入大文件，例如一次一个字节。在这些情况下，可以通过缓冲文件发现极大的速度改进。缓冲时，您将更大的块（缓冲区）读入内存，而当您想从文件中读取时，您会检查它是否在缓冲区中。如果没有，请从文件中将另一个大块读入缓冲区。这种方法最大限度地减少了磁盘延迟。的操作系统和您的编程语言都会使用缓存，但是，有时需要手动缓冲您的文件 (opens new window)。

# 避免缓存未命中

RAM 比磁盘快，而 CPU 又比 RAM 快。 CPU 像时钟一样滴答作响，速度约为 3 GHz，即每秒 30 亿次滴答。该时钟的一个“滴答”称为时钟周期。虽然这是一种简化，但您可以想象每个周期，CPU 都会执行一个称为 CPU 指令的简单命令，该命令对一小段数据执行一次操作。然后时钟速度可以作为计算机中其他时序的参考。值得了解一个时钟周期有多快：在一个周期中，一个光子只能行进 10 厘米左右。事实上，现代 CPU 速度如此之快，以至于对其物理布局的一个重大限制是必须考虑电力通过其内部线路所需的时间，即所谓的线路延迟。

在这种规模下，从 RAM 读取大约需要 500 个时钟周期。类似于如何通过将数据复制到更快的 RAM 来缓解磁盘的高延迟，来自 RAM 的数据被复制到物理上位于 CPU 上的较小内存芯片，称为 cache (缓存)。缓存更快，因为它在物理上位于 CPU 芯片上（减少线路延迟），并且因为它使用更快的存储类型，静态 RAM，而不是速度较慢（但制造成本较低）的动态 RAM，后者是 RAM 的主要构成。因为缓存必须放在 CPU 上，限制了它的大小，并且因为它的生产成本更高，典型的 CPU 缓存只包含大约 $10^{8}$ 位，比 RAM 少大约 1000 倍。 CPU 缓存实际上有多层，但在这里我们对其进行了简化，将“缓存”仅视为一个部件：

$[CPU] \leftrightarrow [CPUCACHE] \leftrightarrow [RAM] \leftrightarrow [DISKCACHE] \leftrightarrow [DISK]$

当 CPU 从 RAM 请求一段数据时，比如单个字节，它会首先检查内存是否已经在缓存中。如果是这样，它将从那里读取它。这比访问 RAM 要快得多，通常只有一个或几个时钟周期。如果没有，我们就会出现缓存未命中，当您的计算机将数据从 RAM 复制到缓存时，您的程序将停止大约 100 纳秒。

除了非常低级的语言之外，不可能手动管理 CPU 缓存。相反，您必须确保有效地使用缓存。

首先，您努力使用尽可能少的内存。内存越少，您的数据就越有可能在 CPU 需要时缓存。请记住，在一次缓存未命中浪费的时间内，CPU 可以执行大约 500 个小操作。

缓存的有效使用归结为局部性，时间和空间局部性：

通过时间局部性，我的意思是您最近访问的数据可能已经驻留在缓存中。因此，如果您必须多次访问一块内存，请确保您在时间上靠近在一起。
通过空间局部性，我的意思是您应该从彼此靠近的内存地址访问数据。您的 CPU 不会仅将请求的字节复制到缓存中。相反，您的 CPU 将始终以称为缓存行的更大块复制数据（通常为 512 个连续位，具体取决于 CPU 型号）。

为了说明这一点，让我们比较一下上面的 random_access 函数在短 (8 KiB) 向量和长 (16 MiB) 向量上运行时的性能。第一个足够小，只需几次加入后，所有数据都已复制到缓存中。第二个是如此之大，以至于新索引在大多数情况下都会导致缓存未命中。

with_terminal() do
  @btime random_access($(rand(UInt, 1024)), 2^20) seconds=1
  @btime random_access($(rand(UInt, 2^24)), 2^20) seconds=1
end

1
2
3
4

  1.772 ms (0 allocations: 0 bytes)
  118.476 ms (0 allocations: 0 bytes)

1
2

我们可以使用之前的函数 random_access。如果我们不是随机访问数组，而是以最坏的顺序访问它，会发生什么？

例如，我们可以使用以下函数：

function linear_access(data::Vector{UInt}, N::Integer)
    n = rand(UInt)
    mask = length(data) - 1
    for i in 1:N
        n = (n >>> 7) ⊻ data[(15 * i) & mask + 1]
    end
    return n
end;

1
2
3
4
5
6
7
8

linear_access 与 random_access 进行几乎相同的计算，但每 15 个元素访问一次。 Julia 中的 UInt 是 8 字节（64 位），因此步长为 15 意味着每个元素之间都有 $15 \times 64 + 960$ 位 -- 大于 64 字节的行缓存。这意味着每次访问都会导致缓存未命中 - 与具有大向量的 random_access 形成对比，其中只有大多数访问会强制缓存未命中。

with_terminal() do
  data = rand(UInt, 2^24)
  @btime random_access($data, 2^20) seconds=1
  @btime linear_access($data, 2^20) seconds=1
end

1
2
3
4
5

  115.172 ms (0 allocations: 0 bytes)
  5.044 ms (0 allocations: 0 bytes)

1
2

惊喜！线性访问模式快 20 倍以上！怎么可能？

CPU 的缓存旁边有一个称为预取器的小电路。该电子电路收集有关 CPU 正在访问哪些内存的数据，并寻找模式。当它检测到一个模式时，它会预取它预测将很快被访问的任何数据，以便在 CPU 请求数据时它已经驻留在缓存中。

我们的函数 linear_access 尽管缓存使用率比 random_access 差，但它以完全可预测的模式获取数据，这允许预取器完成其工作。

总之，我们已经看到

缓存未命中会导致大约 500 次 CPU 操作的损失，因此避免这些对性能至关重要
减少缓存未命中：
- 使用较小的数据，使其更容易放入缓存
- 以可预测的规则模式访问数据以允许预取器完成其工作
- 访问内存中靠近而不是相距很远的数据
- 当访问内存中靠近的数据时，在时间上尽可能靠近，所以当它第二次访问时，它仍然在缓存中。

缓存使用对您的数据结构有影响。诸如 Dicts 和 Sets 之类的哈希表本质上是缓存效率低下的，并且几乎总是导致缓存未命中，而数组则不会。因此，虽然 set 和 dics 的许多操作都是 $\mathbf{O(1)}$ ，但它们的每次操作成本很高。

本文档中的许多优化都会间接影响缓存的使用，因此请务必牢记这一点。

# 保持数据与内存对齐

正如刚才提到的，您的 CPU 将一次将通常 512 个连续位（64 字节）的整个高速缓存行移入和移出主 RAM 到高速缓存。你的整个主内存被分割成缓存行。例如，内存地址 0 到 63 是一个缓存行，地址 64 到 127 是下一个，128 到 191 是下一个，等等。您的 CPU 可能只从内存中请求这些缓存行之一，而不是例如从地址 30 到 93 的 64 个字节。

这意味着某些数据结构可以跨越缓存行之间的边界。如果我在地址 60 请求一个 64 位（8 字节）整数，CPU 必须首先从单个请求的内存地址生成两个内存请求（即获取缓存行 0-63 和 64-127），然后从两个缓存行检索整数，浪费时间。

浪费的时间可能很重要。在高速缓存内存访问证明瓶颈的情况下，减速可能接近 2 倍。在下面的示例中，我使用指针重复访问距缓存线边界给定偏移量处的数组。如果偏移量在 0:56 范围内，则整数都适合单个缓存行，并且函数很快。如果偏移量为 57:63，则所有整数都将跨越缓存行。

function alignment_test(data::Vector{UInt}, offset::Integer)
    # Jump randomly around the memory.
    n = rand(UInt)
    mask = (length(data) - 9) ⊻ 7
    GC.@preserve data begin # protect the array from moving in memory
        ptr = pointer(data)
        iszero(UInt(ptr) & 63) || error("Array not aligned")
        ptr += (offset & 63)
        for i in 1:4096
            n = (n >>> 7) ⊻ unsafe_load(ptr, (n & mask + 1) % Int)
        end
    end
    return n
end;

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

with_terminal() do
  data = rand(UInt, 256 + 8)
  @btime alignment_test($data, 0) seconds=1
  @btime alignment_test($data, 60) seconds=1
end

1
2
3
4
5

  7.043 μs (0 allocations: 0 bytes)
  7.982 μs (0 allocations: 0 bytes)

1
2

未对齐内存访问的后果非常依赖于 CPU。在我当前的 CPU 上，我看到性能下降了约 15%。在我最初编写此笔记本的旧计算机上，惩罚约为 100%。旧处理器可以做更糟糕的事情 (opens new window) - 令人难以置信的是，2001 年的 Game Boy Advance 中的 CPU 会默默地执行不同的读取！ 😱

幸运的是，编译器采用了一些技巧来降低您访问未对齐数据的可能性。首先，Julia（和其他编译语言）总是将新对象放在内存中缓存行的边界处。当一个对象正好放在边界上时，我们说它是对齐的。 Julia 还对齐较大数组的开头：

let
  memory_address = reinterpret(UInt, pointer(rand(1024)))
  @assert iszero(memory_address % 64) # should not error!
end

1
2
3
4

请注意，如果数组的开头是对齐的，那么 1 字节、2 字节、4 字节或 8 字节的对象就不可能跨越缓存线边界，并且所有内容都将对齐。

此外，一个 7 字节的对象仍然可能在数组中未对齐。在 7 字节对象的数组中，第 10 个对象将放置在字节偏移 $7 \times (10-1) = 63$ 处，并且该对象将跨越缓存行。但是，由于这个原因，编译器通常不允许具有非标准大小的结构。如果我们定义一个 7 字节的结构体：

struct AlignmentTest
    a::UInt32 # 4 bytes +
    b::UInt16 # 2 bytes +
    c::UInt8  # 1 byte = 7 bytes?
end

1
2
3
4
5

然后我们可以使用 Julia 的 introspection 来获取 AlignmentTest 对象中三个整数中每一个在内存中的相对位置：

with_terminal() do
  T = AlignmentTest
  println("Size of $T: ", sizeof(T), "bytes")
  for fieldno in 1:fieldcount(T)
    print("Name: ", fieldname(T, fieldno), '\t')
    print("Size: ", sizeof(fieldtype(T, fieldno)), '\t')
    print("Offset: ", fieldoffset(T, fieldno), '\n')
  end
end

1
2
3
4
5
6
7
8
9

Size of Main.workspace2.AlignmentTest: 8bytes
Name: a  Size: 4  Offset: 0
Name: b  Size: 2  Offset: 4
Name: c  Size: 1  Offset: 6

1
2
3
4

我们可以看到，尽管 AlignmentTest 只有 $4 + 2 + 1 = 7$ 个字节的实际数据，但它占用了 8 个字节的内存，并且从数组访问 AlignmentTest 对象将始终是对齐的。

作为编码人员，只有少数情况下您会遇到对齐问题。我可以想出两个：

如果您手动创建一个奇怪大小的对象，例如：通过使用指针访问密集整数数组。这样可以节省内存，但会浪费时间。作者的 Cuckoo 过滤器实现 (opens new window)是为了节省空间。
在矩阵运算过程中。如果你有一个矩阵，列有时是未对齐的，因为它在内存中密集存储。例如。在 Float32s 的 15x15 矩阵中，只有第一列对齐，所有其他列都不对齐。在进行矩阵运算时，这可能会产生严重影响：我见过基准测试 (opens new window)，由于对齐问题，80x80 矩阵/向量乘法比 79x79 乘法快 2 倍。

# 题外话：汇编代码

要运行，任何程序都必须被翻译或编译为 CPU 指令。 CPU 指令是您计算机上实际运行的指令，而不是用您的编程语言编写的代码，后者只是程序的描述。 CPU 指令通常以汇编形式呈现给人类。汇编是一种与 CPU 指令一一对应的编程语言。

查看汇编代码将有助于理解以下与 CPU 指令相关的部分内容。

在 Julia 中，我们可以使用 code_native 函数或等效的 @code_native 宏轻松检查编译后的汇编代码。我们可以为一个简单的函数做到这一点：

# View armasm code generated from this function call
function foo(x)
    s = zero(eltype(x))
    @inbounds for i in eachindex(x)
        s = x[i ⊻ s]
    end
    return s
end;

1
2
3
4
5
6
7
8

with_terminal(() -> @code_native foo([UInt(1)]))

  .text
; ┌ @ pluto.jl#==#a36582d4-8af0-11eb-2b5a-e577c5ed07e2:4 within `foo'
; │┌ @ abstractarray.jl:256 within `eachindex'
; ││┌ @ abstractarray.jl:109 within `axes1'
; │││┌ @ abstractarray.jl:89 within `axes'
; ││││┌ @ array.jl:133 within `size'
  movq  24(%rdi), %rcx
; │└└└└
; │┌ @ range.jl:670 within `iterate'
; ││┌ @ range.jl:519 within `isempty'
; │││┌ @ operators.jl:305 within `>'
; ││││┌ @ int.jl:83 within `<'
  testq  %rcx, %rcx
; │└└└└
  je  L55
; │ @ pluto.jl#==#a36582d4-8af0-11eb-2b5a-e577c5ed07e2:5 within `foo'
; │┌ @ abstractarray.jl:1173 within `getindex' @ array.jl:0
  movq  (%rdi), %rdx
; │└
  negq  %rcx
  movl  $1, %esi
  xorl  %eax, %eax
  nopw  %cs:(%rax,%rax)
; │┌ @ int.jl:923 within `xor' @ int.jl:333
L32:
  xorq  %rsi, %rax
; │└
; │┌ @ range.jl:674 within `iterate'
; ││┌ @ promotion.jl:410 within `=='
  leaq  1(%rcx,%rsi), %rdi
; ││└
  incq  %rsi
; │└
; │┌ @ abstractarray.jl:1173 within `getindex' @ array.jl:801
  movq  -8(%rdx,%rax,8), %rax
; │└
; │┌ @ range.jl:674 within `iterate'
; ││┌ @ promotion.jl:410 within `=='
  cmpq  $1, %rdi
; │└└
  jne  L32
; │ @ pluto.jl#==#a36582d4-8af0-11eb-2b5a-e577c5ed07e2:7 within `foo'
  retq
; │ @ pluto.jl#==#a36582d4-8af0-11eb-2b5a-e577c5ed07e2 within `foo'
L55:
  xorl  %eax, %eax
; │ @ pluto.jl#==#a36582d4-8af0-11eb-2b5a-e577c5ed07e2:7 within `foo'
  retq
  nopw  (%rax,%rax)
; └

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

让我们分解一下：

以 ; 开头的行是注释，并解释以下指令来自代码的哪一部分。它们显示了嵌套的函数调用系列，以及它们在源代码中的位置。可以看到 eachindex，调用 axes1，调用 axes，调用 size。在包含 size 调用的注释行下，我们看到了第一条 CPU 指令。指令名称在最左边，movq。该名称由两部分组成，mov，指令类型（将内容移入或移出寄存器）和后缀 q，“quad”的缩写，表示 64 位整数。有以下后缀：b（byte，8 位），w（word，16 位），l（long，32 位）和 q（quad，64 位）。

指令中的下两列 24(%rdi) 和 %rax 是 movq 的参数。这些是存储要操作的数据的寄存器的名称（稍后我们将返回到寄存器）。

您还可以看到（在汇编代码的较大显示中）代码被分割成以“L”开头的名称开头的部分，例如，当我运行它时，有一个部分 L32。这些部分在使用 if 语句或分支之间跳转。此处，L32 部分标记了实际循环。可以在 L32 部分看到以下两条指令：

; ││┌ @ promotion.jl:401 within `=='
     cmpq    $1, %rdi
; │└└
     jne     L32

1
2
3
4

第一条指令 cmpq（compare quad）将寄存器 rdi 中的数据与数字 1 进行比较，寄存器 rdi 中保存了剩余迭代次数（加一）的数据，并根据结果在 CPU 中设置某些标志（线路）。如果 CPU 中没有设置“相等”标志，则下一条指令 jne（如果不相等则跳转）进行跳转，如果还剩下一个或多个迭代，就会发生这种情况。您可以看到它跳转到 L32，这意味着该部分重复。

# 快指令，慢指令

并非所有 CPU 指令都同样快。下面是一个选定的 CPU 指令表，其中非常粗略地估计了它们执行所需的时钟周期数。您可以在这个文档 (opens new window)中找到更详细的表格。在这里，我将总结现代 Intel CPU 上的指令速度。所有现代 CPU 都非常相似。

您将看到时间以延迟和倒数吞吐量（即 $\dfrac{1}{吞吐量}$ ）的形式给出。原因是 CPU 包含多个用于某些可以并行操作的操作的电路。因此，虽然浮点乘法的延迟为 5 个时钟周期，但如果可以在 10 个不同的电路中并行计算 10 个浮点运算，则其吞吐量为 2 个运算/秒，因此倒数吞吐量为 0.5。

下表以时钟周期为单位测量时间：

Instruction	Latency	Rec. throughp.
move data	1	0.25
and/or/xor	1	0.25
test/compare	1	0.25
do nothing	1	0.25
int add/subtract	1	0.25
bitshift	1	0.5
float multiplication	5	0.5
vector int and/or/xor	1	0.5
vector int add/sub	1	0.5
vector float add/sub	4	0.5
vector float multiplic.	5	0.5
lea	3	1
int multiplic	3	1
float add/sub	3	1
float multiplic.	5	1
float division	15	5
vector float division	13	8
integer division	50	40

lea 指令接受三个输入，A、B 和 C，其中 A 必须是 2、4 或 8，并计算 $AB + C$ 。我们稍后会回到“向量”指令的作用。

为了比较，我们还可以添加一些对其他延迟来源的非常粗略的估计：

Delay	Cycles
move memory from cache	1
misaligned memory read	10
cache miss	500
read from disk	5000000

如果您有一个执行数百万次的内部循环，检查为循环生成的汇编代码并检查您是否可以根据快速 CPU 指令表达计算可能是值得的。例如，如果您有一个已知为 0 或以上的整数，并且想将其除以 8（丢弃任何余数），则可以改为进行位移，因为位移比整数除法快得多：

begin
  divide_slow(x) = div(x, 8)
  divide_fast(x) = x >>> 3;
end;

1
2
3
4

但是，现代编译器非常聪明，并且通常会找出在您的函数中使用的最佳指令以获得相同的结果，例如，通过在适用的情况下将整数除法 idivq 指令替换为右移 (shrq) 以加快速度。您需要自己检查汇编代码才能看到：

# Calling it with debuginfo=:none removes the comments in the armasm code
with_terminal(() -> @code_native debuginfo=:none divide_slow(UInt(1)))

1
2

  .text
  movq  %rdi, %rax
  shrq  $3, %rax
  retq
  nopl  (%rax,%rax)

1
2
3
4
5

# 分配和不变性

如前所述，主 RAM 比 CPU 缓存慢得多。然而，在主 RAM 中工作还有一个额外的缺点：操作系统 (OS) 为程序提供了一堆 RAM 以供使用。在这块内存中，程序本身需要跟踪哪些对象正在使用哪些 RAM。如果程序没有跟踪，在内存中创建一个对象可能会覆盖另一个对象，从而导致数据丢失。因此，RAM 中数据的每次创建和销毁都需要记录，这需要时间。

在 RAM 中创建新对象称为分配，销毁称为释放。实际上，（解除）分配本身并不是真正的创建或销毁，而是开始和停止跟踪内存的行为。未被跟踪的内存最终将被其他数据覆盖。分配和解除分配需要大量时间，具体取决于对象的大小，每次分配从几十纳秒到几微秒不等。

在 Julia、Python、R 和 Java 等编程语言中，释放是使用称为垃圾收集器 (GC) 的程序自动完成的。该程序会跟踪哪些对象是程序员无法访问的，并释放它们。例如，如果你这样做：

begin
  thing = [1,2,3]
  thing = nothing
end

1
2
3
4

那么就没有办法把原来的数组[1,2,3]取回来了，是不可达的。因此，它只是在浪费 RAM，什么也不做。它是垃圾. 分配和释放对象有时会导致 GC 开始扫描内存中的所有对象并释放无法访问的对象，这会导致明显的滞后。您也可以手动启动垃圾收集器：

GC.gc()

以下示例说明了在分配具有新结果的向量的函数中相对于仅修改向量而不分配任何内容的函数所花费的时间的差异：

begin
  function increment(x::Vector{<:Integer})
    y = similar(x)
    @inbounds for i in eachindex(x)
      y[i] = x[i] + 1
    end
    return y
  end

  function increment!(x::Vector{<:Integer})
    @inbounds for i in eachindex(x)
      x[i] = x[i] + 1
    end
    return x
  end
end;

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

with_terminal() do
  data = rand(UInt, 2^10)
  show(stdout, MIME"text/plain"(), @benchmark increment($data) seconds=1)
  println('\n')
  show(stdout, MIME"text/plain"(), @benchmark increment!($data) seconds=1)
end

1
2
3
4
5
6

BenchmarkTools.Trial: 4362 samples with 153 evaluations.
 Range (min … max):  664.523 ns … 35.843 μs  ┊ GC (min … max):  0.00% … 94.88%
 Time  (median):       1.137 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.457 μs ±  2.096 μs  ┊ GC (mean ± σ):  12.61% ±  8.56%

        █
  ▅▅▄▄▄██▇▅▄▄▃▃▃▃▃▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▂▂▁▂▂▂▂▁▁▁▂▂▂▂ ▃
  665 ns          Histogram: frequency by time         4.99 μs <

 Memory estimate: 8.12 KiB, allocs estimate: 1.

BenchmarkTools.Trial: 10000 samples with 980 evaluations.
 Range (min … max):  63.620 ns … 287.154 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     66.411 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   66.658 ns ±   3.548 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

          ▁  ▃  ▂  ▇  █ ▁▄                     ▁               ▁
  ▇▁█▁▁▇▁▁█▁▄█▁▆█▁▁█▆▇█▆██▅▇▆▆▅▇▇▆██▇▇▇▇▇▇▆▇█▆██▆██▆▅▄▁▄▅▄▄▄▇▇ █
  63.6 ns       Histogram: log(frequency) by time      71.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

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

在我的计算机上，分配函数平均慢 20 倍以上。还要注意在分配函数上花费的最大时间很高。这是由于代码的一些属性：

首先，分配本身需要时间
其次，分配的对象最终必须被释放，这也需要时间
第三、重复分配触发 GC 运行，造成开销
第四，更多的分配有时意味着较低的缓存使用效率，因为您使用了更多的内存

请注意，我使用了平均时间而不是中位数，因为对于这个函数，GC 大约每 30 次调用才触发一次，但它会消耗 30-40 微秒。所有这些都意味着高性能代码应该将分配保持在最低限度。

@btime 宏和其他基准测试工具会告诉您分配的数量。提供此信息是因为假设任何关心对其代码进行基准测试的程序员都会对减少分配感兴趣。

# 并非所有对象都需要分配

在 RAM 中，数据保存在栈或堆上。堆栈是一个简单的数据结构，有开始和结束，类似于 Julia 中的 Vector。堆栈只能通过从末尾添加或减去元素来修改，类似于只有两个变异操作 push! 和 pop! 的 Vector。堆栈上的这些操作非常快。然而，当我们谈论“分配”时，我们谈论的是堆上的数据。与栈不同，堆具有无限大小（好吧，它具有计算机 RAM 的大小），并且可以任意修改，删除和访问任何对象。你可以把栈想象成一个 Vector，把堆想象成一个 Dict。

从直觉上看，所有对象都需要放在 RAM 中，必须能够随时被程序检索和删除，这似乎很明显，因此需要在堆上分配。对于某些语言，例如 Python，这是正确的。然而，这在 Julia 和其他高效的编译语言中并非如此。例如，整数通常可以放在栈中。

为什么有些对象需要堆分配，而有些对象可以栈分配？要进行栈分配，编译器需要确定：

该对象是一个相当小的尺寸，所以它适合栈。由于技术原因，栈的大小不能只有数百兆。
编译器可以准确预测何时需要创建和销毁对象，因此它可以通过简单地弹出栈及时销毁它（类似于在 Vector 上调用 pop!）。这通常是编译语言中的局部变量的情况。

Julia 对栈分配的对象有更多限制。

该对象应具有在编译时已知的固定大小。
编译器必须知道对象永远不会改变。 CPU 可以自由复制堆栈分配的对象，对于不可变对象，无法区分副本和原始对象。这值得重申：对于不可变对象，无法区分副本和原始对象。这给了编译器和 CPU 在操作时一定的自由。后一点也是为什么对象在 Julia 中默认是不可变的，并导致另一个性能提示：尽可能使用不可变对象。

这在实践中意味着什么？在 Julia 中，这意味着如果您想要快速的栈分配对象：

您的对象必须在完全编译的函数中创建、使用和销毁，以便编译器确定何时需要创建、使用和销毁对象。如果对象返回供以后使用（而不是立即返回到另一个完全编译的函数），我们说该对象逃逸，并且必须被分配。
您的类型必须有大小限制。我不知道它到底有多大，但 100 字节就可以了。
编译器必须知道您类型的确切内存布局（几乎总是如此）。

注意：本文的早期版本提到 Julia 中的栈分配类型也必须是所谓的位类型。此限制在 Julia 1.5 中取消

begin
  abstract type AllocatedInteger end

  struct StackAllocated <: AllocatedInteger
    x::Int
  end

  mutable struct HeapAllocated <: AllocatedInteger
    x::Int
  end
end

1
2
3
4
5
6
7
8
9
10
11

我们可以使用实例化 StackAllocated 对象所需的代码检查实例化 HeapAllocated 对象所需的代码：

with_terminal(() -> @code_native HeapAllocated(1))

  .text
; ┌ @ pluto.jl#==#2a7c1fc6-8af1-11eb-2909-554597aa2949:9 within `HeapAllocated'
  pushq  %rbx
  movq  %rdi, %rbx
  movq  %fs:0, %rdi
  addq  $-32768, %rdi                   # imm = 0x8000
  movabsq  $139915596420025, %rax          # imm = 0x7F40A36C53B9
  movl  $1400, %esi                     # imm = 0x578
  movl  $16, %edx
  callq  *%rax
  movabsq  $139913486204208, %rcx          # imm = 0x7F4025A4FD30
  movq  %rcx, -8(%rax)
  movq  %rbx, (%rax)
  popq  %rbx
  retq
  nopl  (%rax)
; └

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

注意 HeapAllocated 中的 callq 指令。该指令调用其他函数，这意味着实际上需要更多代码来创建显示的 HeapAllocated 对象。相比之下，StackAllocated 实际上只需要几条指令：

with_terminal(() -> @code_native StackAllocated(1))

  .text
; ┌ @ pluto.jl#==#2a7c1fc6-8af1-11eb-2909-554597aa2949:5 within `StackAllocated'
  movq  %rdi, %rax
  retq
  nopw  %cs:(%rax,%rax)
; └

1
2
3
4
5
6

因为不可变对象不需要存储在堆上并且可以自由复制，所以不可变对象内联存储在数组中。这意味着不可变对象可以直接存储在数组的内存中。可变对象在堆上具有唯一标识和位置。它们与副本区分开来，因此不能自由复制，因此数组包含对存储它们的堆上内存位置的引用。从数组访问这样的对象意味着首先访问数组以获取内存位置，然后使用该内存位置访问对象本身。除了双内存访问之外，对象在堆上的存储效率较低，这意味着需要将更多内存复制到 CPU 缓存，这意味着更多的缓存未命中。因此，即使以数组的形式存储在堆上，也可以更有效地存储不可变项。

let
  Base.:+(x::Int, y::AllocatedInteger) = x + y.x
  Base.:+(x::AllocatedInteger, y::AllocatedInteger) = x.x + y.x

  data_stack = [StackAllocated(i) for i in rand(UInt16, 1000000)]
  data_heap = [HeapAllocated(i.x) for i in data_stack]

  with_terminal() do
    @btime sum($data_stack) seconds=1
    @btime sum($data_heap) seconds=1
  end
end

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

  365.630 μs (0 allocations: 0 bytes)
  1.086 ms (0 allocations: 0 bytes)

1
2

我们可以验证，确实，data_stack 中的数组存储了 StackAllocated 对象的实际数据，而 data_heap 包含指针（即内存地址）：

with_terminal() do
  data_stack = [StackAllocated(i) for i in rand(UInt16, 1)]
  data_heap = [HeapAllocated(i.x) for i in data_stack]

  println(rpad("First object of data_stack:", 36), data_stack[1])
  println(
    rpad("First data in data_stack array:", 36),
    unsafe_load(pointer(data_stack)),
    '\n'
  )

  println(rpad("First object of data_heap:", 36), data_heap[1])
  first_data = unsafe_load(Ptr{UInt}(pointer(data_heap)))
  println(rpad("First data in data_heap array:", 36), repr(first_data))
  println(
    "Data at address ",
    repr(first_data), ": ",
    unsafe_load(Ptr{HeapAllocated}(first_data))
  )
end

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

First object of data_stack:         Main.workspace2.StackAllocated(39577)
First data in data_stack array:     Main.workspace2.StackAllocated(39577)

First object of data_heap:          Main.workspace2.HeapAllocated(39577)
First data in data_heap array:      0x00007f4032b70500
Data at address 0x00007f4032b70500: Main.workspace2.HeapAllocated(39577)

1
2
3
4
5
6

# 寄存器和 SIMD

又到了更新我们简化的计算机原理图的时候了。 CPU 仅对存在于寄存器中的数据进行操作。这些是 CPU 本身内部的小型、固定大小的插槽（例如大小为 8 个字节）。寄存器旨在保存单个数据，例如整数或浮点数。正如汇编代码部分所暗示的那样，每条指令通常指的是一两个寄存器，这些寄存器包含操作所处理的数据：

$[CPU] \leftrightarrow [REGISTERS] \leftrightarrow [CPUCACHE] \leftrightarrow [RAM] \leftrightarrow [DISKCACHE] \leftrightarrow [DISK]$

要对大于一个寄存器的数据结构进行操作，必须将数据分解为适合寄存器的较小部分。例如，在我的计算机上添加两个 128 位整数时：

with_terminal(() -> @code_native UInt128(5) + UInt128(11))

  .text
; ┌ @ int.jl:87 within `+'
  addq  %rcx, %rsi
  movq  %rdi, %rax
  adcq  %r8, %rdx
  movq  %rsi, (%rdi)
  movq  %rdx, 8(%rdi)
  retq
  nopw  %cs:(%rax,%rax)
; └

1
2
3
4
5
6
7
8
9
10

没有可以进行 128 位加法的寄存器。首先，必须使用 addq 指令添加低 64 位，并放入寄存器中。然后用 adcq 指令添加高位，该指令添加数字，但也使用上一条指令的进位位。最后，使用 movq 指令将结果一次移动 64 位。

寄存器的小尺寸成为 CPU 吞吐量的瓶颈：它一次只能操作一个整数/浮点数。为了避免这种情况，现代 CPU 包含专门的 256 位寄存器（或在旧 CPU 中为 128 位，或在全新 CPU 中为 512 位），可以同时容纳 4 个 64 位整数/浮点数，或 8 个 32 位整数等。令人困惑的是，这种宽寄存器中的数据被称为“向量”。 CPU 可以访问可以对向量执行各种 CPU 操作的指令，在一条指令中对 4 个 64 位整数进行操作。这称为“单指令、多数据”、SIMD 或矢量化。值得注意的是，4x64 位操作与 256 位操作不同，例如当您添加两个向量时，4 个 64 位整数之间没有结转。相反，256 位向量运算等效于 4 个单独的 64 位运算。

我们可以用下面的例子来说明这一点：

with_terminal() do
  # Create a single statically-sized vector of 8 32-bit integers
  # I could also have created 4 64-bit ones, etc.
  a = @SVector Int32[1,2,3,4,5,6,7,8]
  @code_native debuginfo=:none a + a
end

1
2
3
4
5
6

  .text
  vmovdqu  (%rdx), %ymm0
  movq  %rdi, %rax
  vpaddd  (%rsi), %ymm0, %ymm0
  vmovdqu  %ymm0, (%rdi)
  vzeroupper
  retq
  nopw  %cs:(%rax,%rax)

1
2
3
4
5
6
7
8

在这里，两个 8×32 位向量在一条指令中相加。可以看到 CPU 利用一条 vpaddd（vectorpacked add double）指令将 8 个 32 位整数相加，以及相应的移动指令 vmovdqu。请注意，向量 CPU 指令以 v 开头。

值得一提的是 SIMD 和对齐之间的交互：如果一系列 256 位（32 字节）SIMD 加载未对齐，那么多达一半的加载可能会跨越缓存线边界，而只是 8 字节的 1/8 负载。因此，在使用 SIMD 时，对齐是一个更严重的问题。由于数组的开头总是对齐的，这通常不是问题，但在不能保证从对齐的起点开始的情况下，例如矩阵运算，这可能会产生显着差异。在具有 512 位寄存器的全新 CPU 中，问题更严重，因为 SIMD 大小与缓存行大小相同，因此如果初始加载错位，所有加载都会错位。

SIMD 向量化，例如 64 位整数可以将吞吐量提高近 4 倍，因此它在高性能编程中非常重要。如果可以，编译器将自动矢量化操作。什么可以阻止这种自动矢量化？

# SIMD 需要固定长度的不间断迭代

因为向量化操作一次对多个数据进行操作，所以不可能在任意点中断循环。例如，如果在一个时钟周期内处理了 4 个 64 位整数，那么在处理完 3 个整数后就无法停止 SIMD 循环。假设你有一个这样的循环：

for i in 1:8
    if foo()
        break
    end
    # do stuff with my_vector[i]
end

1
2
3
4
5
6

在这里，由于 break 语句，循环可以在任何迭代中结束。因此，任何加载多个整数的 SIMD 指令都可以在循环应该中断后对数据进行操作，即永远不应该读取的数据。这将是错误的行为，因此编译器不能使用 SIMD 指令。

一个好的经验法则是 SIMD 需要：

具有预定长度的循环，因此它知道何时停止
循环中没有分支（即 if 语句）的循环

事实上，即使是边界检查，即检查您没有在向量边界之外进行索引，也会导致分支。毕竟，如果代码应该在 3 次迭代后引发边界错误，那么即使是单个 SIMD 操作也会出错！为了实现 SIMD 矢量化，必须禁用所有边界检查。我们可以用这个来演示 SIMD 的影响：

begin
  function sum_boundscheck(x::Vector)
    n = zero(eltype(x))
    for i in eachindex(x)
      n += x[i]
    end
    return n
  end

  function sum_inbounds(x::Vector)
    n = zero(eltype(x))
    # By removing the boundscheck, we allow automatic SIMD
    @inbounds for i in eachindex(x)
      n += x[i]
    end
    return n
  end
end;

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

with_terminal() do
  # Make sure the vector is small so we don't time cache misses
  data = rand(UInt64, 4096)
  @btime sum_boundscheck($data) seconds=1
  @btime sum_inbounds($data) seconds=1
end

1
2
3
4
5
6

  1.411 μs (0 allocations: 0 bytes)
  137.448 ns (0 allocations: 0 bytes)

1
2

在我的电脑上，SIMD 代码比非 SIMD 代码快 10 倍。 SIMD 仅占大约 4 倍的改进（因为我们从每次迭代 64 位移动到每次迭代 256 位）。其余的收益来自不花时间检查边界和自动循环展开（稍后解释），@inbounds 注释也使这成为可能。

# SIMD 需要一个循环，其中循环顺序无关紧要

SIMD 可以更改处理数组中元素的顺序。如果任何迭代的结果取决于任何先前的迭代，使得元素无法重新排序，则编译器通常不会进行 SIMD 向量化。通常，当循环不会自动矢量化时，这是由于数据在寄存器中移动的微妙之处，这意味着数组中的元素之间会有一些隐藏的内存依赖性。

想象一下，我们想使用 SIMD 对数组中的一些 64 位整数求和。为简单起见，假设数组有 8 个元素，A、B、C ... H。在普通的非 SIMD 循环中，相加会像这样完成：

$(((((((A+B)+C)+D)+E)+F)+G)+H)$

而当使用 SIMD 加载整数时，四个 64 位整数将加载到一个向量 <A, B, C, D> 中，而其他四个则加载到另一个 <E, F, G, H> 中。将添加两个向量：<A+E、B+F、C+G、D+H>。在循环之后，将求和结果向量中的四个整数。所以总的顺序是：

$((((A+E)+(B+F))+(C+G))+(D+H))$

也许令人惊讶的是，浮点数的加法可以根据顺序给出不同的结果（即浮点加法不是关联的）：

false:

begin
  x = eps(1.0) * 0.4
  1.0 + (x + x) == (1.0 + x) + x
end

1
2
3
4

因此，浮点加法不会自动矢量化：

with_terminal() do
  data = rand(Float64, 4096)
  @btime sum_boundscheck($data) seconds=1
  @btime sum_inbounds($data) seconds=1
end

1
2
3
4
5

  2.930 μs (0 allocations: 0 bytes)
  2.925 μs (0 allocations: 0 bytes)

1
2

然而，高性能编程语言通常会提供一个命令来告诉编译器可以对循环重新排序，即使对于非关联循环也是如此。在 Julia 中，这个命令是 @simd 宏：

function sum_simd(x::Vector)
    n = zero(eltype(x))
    # Here we add the `@simd` macro to allow SIMD of floats
    @inbounds @simd for i in eachindex(x)
        n += x[i]
    end
    return n
end;

1
2
3
4
5
6
7
8

with_terminal() do
  data = rand(Float64, 4096)
  @btime sum_inbounds($data) seconds=1
  @btime sum_simd($data) seconds=1
end

1
2
3
4
5

  2.943 μs (0 allocations: 0 bytes)
  188.739 ns (0 allocations: 0 bytes)

1
2

Julia 还提供了宏 @simd ivdep，它进一步告诉编译器在循环顺序中没有内存依赖性。然而，我强烈反对使用这个宏，除非你真的知道你在做什么。一般而言，编译器最了解循环何时具有内存依赖关系，滥用 @simd ivdep 很容易导致难以检测的错误。

# 数组结构

如果我们创建一个包含四个 AlignmentTest 对象 A、B、C 和 D 的数组，这些对象将在数组中首尾相连，如下所示：

Objects: |      A        |       B       |       C       |        D      |
Fields:  |   a   | b |c| |   a   | b |c| |   a   | b |c| |   a   | b |c| |
Byte:     1               9              17              25              33

1
2
3

再次注意字节号。 8、16、24、32 为空以保持对齐，浪费内存。现在假设您想对结构的所有 .a 字段进行操作。由于 .a 字段分散为 8 个字节，因此 SIMD 操作的效率（一次最多加载 4 个字段）比将所有 .a 字段存储在一起（其中 8 个字段可以放入 256 位寄存器）低得多. 仅使用 .a 字段时，将读入整个 64 字节缓存行，其中只有一半或 32 字节有用。这不仅会导致更多的缓存未命中，我们还需要指令从我们需要的 SIMD 寄存器中挑选出一半的数据。

我们上面的内存结构被称为“结构数组”，因为它是一个充满结构的数组。相反，我们可以将 4 个对象 A 到 D 构建为“数组结构”。从概念上讲，它可能看起来像：

struct AlignmentTestVector
    a::Vector{UInt32}
    b::Vector{UInt16}
    c::Vector{UInt8}
end

1
2
3
4
5

对齐不再是问题，填充上不会浪费任何空间。当遍历所有 a 字段时，所有缓存行都包含完整的 64 字节相关数据，因此 SIMD 操作不需要额外的操作来挑选相关数据：

function Base.rand(::Type{AlignmentTest})
  AlignmentTest(rand(UInt32), rand(UInt16), rand(UInt8))
end;

1
2
3

with_terminal() do
  N  = 1_000_000
  array_of_structs = [rand(AlignmentTest) for i in 1:N]
  struct_of_arrays = AlignmentTestVector(
    rand(UInt32, N),
    rand(UInt16, N),
    rand(UInt8, N)
  )

  @btime sum(x -> x.a, $array_of_structs) seconds=1
  @btime sum($struct_of_arrays.a) seconds=1
end

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

  434.659 μs (0 allocations: 0 bytes)
  110.889 μs (0 allocations: 0 bytes)

1
2

# 专用 CPU 指令

大多数代码仅使用一些简单的 CPU 指令，例如移动、加法、乘法、位移和、或、异或、跳转等。但是，典型的现代笔记本电脑中的 CPU 支持大量 CPU 指令。通常，如果某项操作在消费类笔记本电脑中大量使用，CPU 制造商会添加专门的指令来加速这些操作。根据指令的硬件实现，使用这些指令的速度增益可能非常显着。

Julia 只公开了一些专门的指令，包括：

整数中设置位的数量使用 popcnt 指令有效计数，通过 count_ones 函数公开。
tzcnt 指令计算位中尾随零的数量一个整数，通过trailing_zeros` 函数公开
可以使用通过 bswap 函数公开的 bswap 指令反转多字节整数中各个字节的顺序。这在必须处理字节序 (opens new window)时很有用。

以下示例说明了 count_ones 函数的手动实现与使用 popcnt 指令的内置版本之间的性能差异：

function manual_count_ones(x)
    n = 0
    while x != 0
        n += x & 1
        x >>>= 1
    end
    return n
end;

1
2
3
4
5
6
7
8

with_terminal() do
  data = rand(UInt, 10000)
  @btime sum(manual_count_ones, $data) seconds=1
  @btime sum(count_ones, $data) seconds=1
end

1
2
3
4
5

  339.350 μs (0 allocations: 0 bytes)
  1.654 μs (0 allocations: 0 bytes)

1
2

您在此处观察的时间将取决于您的编译器是否足够聪明以实现第一个函数中的计算可以表示为 popcnt 指令，因此将被编译为该指令。在我的计算机上，编译器无法进行推断，第二个函数以 100 倍以上的速度实现相同的结果。

# 调用任何 CPU 指令

Julia 使直接调用 CPU 指令成为可能。通常不建议这样做，因为并非所有用户都可以使用相同的指令访问同一个 CPU，因此您的代码会在使用不同品牌计算机的用户上崩溃。

最新的 CPU 包含用于 AES 加密和 SHA256 哈希的专门指令。如果你想调用这些指令，你可以直接调用 Julia 的后端编译器 LLVM。在下面的示例中，我创建了一个直接调用 vaesenc（一轮 AES 加密）指令的函数：

begin
  # This is a 128-bit CPU "vector" in Julia
  const __m128i = NTuple{2, VecElement{Int64}}

  # Define the function in terms of LLVM instructions
  aesenc(a, roundkey) = ccall(
    "llvm.x86.aesni.aesenc", llvmcall, __m128i,
    (__m128i, __m128i), a, roundkey
  )
end;

1
2
3
4
5
6
7
8
9
10

我们可以通过检查函数的汇编来验证它的工作原理，该函数应该只包含一条 vaesenc 指令，以及 retq（返回）和 nopw（什么都不做，用作对齐 CPU 指令在内存中的填充）指令：

with_terminal(() -> @code_native aesenc(__m128i((1, 1)), __m128i((1, 1))))

  .text
; ┌ @ pluto.jl#==#25a47c54-8af2-11eb-270a-5b58c3aafe6e:6 within `aesenc'
  vaesenc  %xmm1, %xmm0, %xmm0
  retq
  nopw  %cs:(%rax,%rax)
; └

1
2
3
4
5
6

使用专门指令的算法可以非常快。在一篇博客文章 (opens new window)中，视频游戏公司 Molly Rocket 推出了一种使用 AES 指令的新非加密哈希函数，该函数达到了前所未有的速度。

# 内联

考虑这个函数的汇编：

f() = error();

with_terminal(() -> @code_native f())

  .text
; ┌ @ pluto.jl#==#36b723fc-8ee9-11eb-1b92-451b992acc0c:1 within `f'
  pushq  %rax
  movabsq  $error, %rax
  callq  *%rax
  ud2
  nop
; └

1
2
3
4
5
6
7
8

此代码包含调用另一个函数的 callq 指令。根据函数的参数和其他东西，函数调用会带来一些开销。虽然花费在函数调用上的时间以纳秒为单位进行测量，但如果调用的函数处于紧密循环中，则它可以累加。

但是，如果我们显示此函数的汇编：

call_plus(x) = x + 1;

with_terminal(() -> @code_native debuginfo=:none call_plus(1))

  .text
  leaq  1(%rdi), %rax
  retq
  nopw  %cs:(%rax,%rax)

1
2
3
4

函数 call_plus 调用 +，并编译为单个 leaq 指令（以及一些填充符 retq 和 nopw）。但是 + 是一个普通的 Julia 函数，所以 call_plus 是一个普通 Julia 函数调用另一个函数的例子。为什么没有 callq 指令来调用 +？

编译器选择将函数 + 内联到 call_plus 中。这意味着它没有调用 +，而是将 + 的内容直接复制到 call_plus 中。这样做的好处是：

函数调用没有开销
不需要构造一个元组来保存 + 函数的参数
+ 中发生的任何计算都与 call_plus 一起编译，允许编译器使用另一个中的信息并可能简化一些计算。

那么为什么不是所有的函数都内联呢？内联复制代码，增加它的大小并消耗 RAM。此外，CPU 指令本身也需要适合 CPU 缓存（尽管 CPU 指令有自己的缓存）才能有效地检索。如果一切都被内联，程序将是巨大的并且会停止。只有内联函数很小的时候，内联才只是一种优化。

相反，编译器使用启发式（经验法则）来确定函数何时足够小以进行内联以提高性能。这些启发式方法不是万无一失的，因此 Julia 提供了宏 @noinline，它可以防止小函数的内联（对于引发错误的函数很有用，必须假设它们很少被调用）和 @inline，它不会强制编译器内联，但强烈建议编译器应该内联该函数。

如果代码包含对时间敏感的部分，例如内部循环，则查看汇编代码以验证循环中的小函数是否已内联非常重要。例如，在作者的 kmer 哈希代码的这一行 (opens new window)中，如果删除此 @inline 注释，整体 minhashing 性能将下降两倍。

因此可以证明内联和不内联之间的极端差异：

begin
  @noinline noninline_poly(x) = x^3 - 4x^2 + 9x - 11
  inline_poly(x) = x^3 - 4x^2 + 9x - 11

  function time_function(F, x::AbstractVector)
    n = 0
    for i in x
      n += F(i)
    end
    return n
  end
end;

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

with_terminal() do
  data = rand(UInt, 1024)
  @btime time_function(noninline_poly, $data) seconds=1
  @btime time_function(inline_poly, $data) seconds=1
end

1
2
3
4
5

  1.240 μs (0 allocations: 0 bytes)
  625.444 ns (0 allocations: 0 bytes)

1
2

# 展开

考虑一个对 64 位整数向量求和的函数。如果向量的数据的内存偏移量存储在寄存器 %r9 中，向量的长度存储在寄存器 %r8 中，向量的当前索引存储在 %rcx 中，运行总数存储在 %rax 中，则内循环的汇编代码看起来像这样：

L1:
    ; add the integer at location %r9 + %rcx * 8 to %rax
    addq   (%r9,%rcx,8), %rax

    ; increment index by 1
    addq   $1, %rcx

    ; compare index to length of vector
    cmpq   %r8, %rcx

    ; repeat loop if index is smaller
    jb     L1

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

向量的每个元素总共有 4 条指令。 Julia 生成的实际代码与此类似，但还包括与边界检查相关的额外指令，这些指令与此处无关（当然会包括不同的注释）。

但是，如果函数是这样写的：

function sum_vector(v::Vector{Int})
    n = 0
    i = 1
    for chunk in 1:div(length(v), 4)
        n += v[i + 0]
        n += v[i + 1]
        n += v[i + 2]
        n += v[i + 3]
        i += 4
    end
    return n
end

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

如果我们假设向量的长度可以被 4 整除，结果显然是一样的。如果长度不能被 4 整除，我们可以简单地使用上面的函数对第一个元素求和，然后在另一个循环中添加最后几个元素。尽管结果功能相同，但循环的汇编不同，可能看起来像：

L1:
    addq   (%r9,%rcx,8), %rax
    addq   8(%r9,%rcx,8), %rax
    addq   16(%r9,%rcx,8), %rax
    addq   24(%r9,%rcx,8), %rax
    addq   $4, %rcx
    cmpq   %r8, %rcx
    jb     L1

1
2
3
4
5
6
7
8

每 4 次添加总共 7 条指令，或每次添加 1.75 条指令。这不到每个整数指令数的一半！速度增益来自于在循环结束时减少检查频率。我们称这个过程为展开循环，这里是四倍。当然，只有在我们事先知道迭代次数的情况下才能展开展开，因此我们不会“过冲”迭代次数。通常，编译器会自动展开循环以获得额外的性能，但值得查看程序集。例如，这是在我的计算机上为 sum([1]) 生成的最内层循环的程序集：

L144:
    vpaddq  16(%rcx,%rax,8), %ymm0, %ymm0
    vpaddq  48(%rcx,%rax,8), %ymm1, %ymm1
    vpaddq  80(%rcx,%rax,8), %ymm2, %ymm2
    vpaddq  112(%rcx,%rax,8), %ymm3, %ymm3
    addq    $16, %rax
    cmpq    %rax, %rdi
    jne L144

1
2
3
4
5
6
7
8

您可以看到它展开了四倍，并使用了 256 位 SIMD 指令，总共 128 个字节，每次迭代添加 16 个整数，或每个整数 0.44 条指令。

另请注意，编译器选择使用 4 个不同的 ymm SIMD 寄存器，ymm0 到 ymm3，而在我的示例汇编代码中，我只使用了一个寄存器 rax。这是因为，如果您使用 4 个独立的寄存器，那么在 CPU 可以开始下一个之前，您不需要等待一个 vpaddq 完成（请记住，它有大约 3 个时钟周期的延迟）。

展开与 SIMD 的相似：编译器只会展开循环，前提是它可以确定不会超过迭代次数。例如，比较以下两个函数：

with_terminal() do
  data = fill(false, 2^20)

  # any: Stops as soon as it finds a `true`
  @btime any($data) seconds=1

  # reduce: Loops over all values in the array
  @btime reduce(|, $data) seconds=1

  data[1] = true
  @btime any($data) seconds=1
  @btime reduce(|, $data) seconds=1
end

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

  258.567 μs (0 allocations: 0 bytes)
  79.610 μs (0 allocations: 0 bytes)
  1.703 ns (0 allocations: 0 bytes)
  79.870 μs (0 allocations: 0 bytes)

1
2
3
4

第一个函数在找到 true 值后立即停止并返回 - 但循环中断会禁用 SIMD 和展开。第二个函数在整个数组中继续运行，即使第一个值为 true。虽然这可以启用 SIMD 和展开，但如果一开始就看到正确的 true，这显然是浪费。因此，当我们期望在数组的 1/4 左右之前看到第一个 true 时，第一个更好，否则后者更好。

我们可以通过手动展开来达成妥协。在下面的函数中，check128 使用 inbounds 检查 128 个条目，而不会停下来检查它是否为 true，从而展开和 SIMDd。 unroll_compromise 然后使用 check128，但一旦发现 true 就会跳出循环。

begin
  @inline function check128(data, i)
      n = false
      @inbounds for j in 0:127
        n |= data[i+j]
      end
      n
    end

  function unroll_compromise(data)
    found = false
    i = 1
    while !found & (i ≤ length(data))
      check128(data, i) && return true
      i += 128
    end
    return false
  end
end;

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

with_terminal() do
  data = fill(false, 2^20)
  @btime reduce(|, $data) seconds=1
  @btime unroll_compromise($data) seconds=1

  data[1] = true
  @btime any($data) seconds=1
  @btime unroll_compromise($data) seconds=1
end

1
2
3
4
5
6
7
8
9

  80.562 μs (0 allocations: 0 bytes)
  11.161 μs (0 allocations: 0 bytes)
  1.462 ns (0 allocations: 0 bytes)
  1.693 ns (0 allocations: 0 bytes)

1
2
3
4

我们看到两个没有 true 值的数组和一开始有 true 值的数组都有出色的性能。

不幸的是，我不知道有什么方法可以自动生成这种展开，您希望在展开较小的块和在每个块之间包含分支之间进行折衷。或许在未来，这个愿望可以传达给编译器，从而自动生成最优代码。

# 避免不可预测的分支

如前所述，CPU 指令需要多个周期才能完成，但可能会在前一条指令完成计算之前排队进入 CPU。那么当 CPU 遇到分支（即跳转指令）时会发生什么？它无法知道接下来要排队的指令，因为这取决于它刚刚放入队列的指令以及尚未执行的指令。

现代 CPU 使用分支预测。 CPU 有一个分支预测器电路，它根据最近采用的分支来猜测正确的分支。从本质上讲，分支预测器试图学习在代码运行时在代码中采用分支的简单模式。将分支排队后，CPU 立即将来自分支预测器预测的任何分支的指令排队。稍后在执行排队分支时验证猜测的正确性。如果猜测是正确的，那太好了，CPU 通过猜测节省了时间。如果没有，CPU 必须清空管道并丢弃自初始猜测以来的所有计算，然后重新开始。此过程会导致几纳秒的延迟。

对于程序员来说，这意味着 if 语句的速度取决于猜测的难易程度。如果很容易猜到，分支预测器几乎一直都是正确的，并且 if 语句不会超过一条简单的指令，通常是 1 个时钟周期。在分支是随机的情况下，大约有 50% 的时间会出错，并且每次错误预测可能会花费许多时钟周期。

由循环引起的分支是最容易猜测的。如果您有一个包含 1000 个元素的循环，则代码将循环回 999 次并仅跳出一次循环。因此，分支预测器总是可以简单地预测“环回”，并获得 99.9% 的准确率。

我们可以用一个简单的函数来演示分支错误预测的性能：

# Copy all odd numbers from src to dst.
function copy_odds_branches!(dst::Vector{T}, src::Vector{T}) where {T <: Integer}
    write_index = 1
    @inbounds for i in eachindex(src) # <--- this branch is trivially easy to predict
        v = src[i]
        if isodd(v)  # <--- this is the branch we want to predict
            dst[write_index] = v
            write_index += 1
        end
    end
    return dst
end;

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

with_terminal() do
  dst = rand(UInt32, 2^18)
  src_random = rand(UInt32, 2^18)
  src_all_odd = [(2*i+1) % UInt32 for i in src_random]
  @btime copy_odds_branches!($dst, $src_random) seconds=1
  @btime copy_odds_branches!($dst, $src_all_odd) seconds=1
end

1
2
3
4
5
6
7

  900.416 μs (0 allocations: 0 bytes)
  111.029 μs (0 allocations: 0 bytes)

1
2

在第一种情况下，整数是随机的，大约一半的分支将被错误预测导致延迟。在第二种情况下，分支总是被采用，分支预测器能够快速获取模式并将达到接近 100% 的正确预测。结果，在我的电脑上，后者快了大约 8 倍。

请注意，如果您使用较小的向量并多次重复计算，就像 @btime 宏所做的那样，分支预测器能够记住小随机向量的模式，并且会比随机预测好得多。这在最现代的 CPU（尤其是 AMD 出售的 CPU，我听说）中尤其明显，其中分支预测器变得更好。这种“用心学习”是基准测试过程中循环的产物。您不会期望在真实数据上重复运行完全相同的计算：

with_terminal() do
  src_random = rand(UInt32, 128)
  dst = similar(src_random)
  src_all_odd = [(2i+1) % UInt32 for i in src_random]

  @btime copy_odds_branches!($dst, $src_random) seconds=1
  @btime copy_odds_branches!($dst, $src_all_odd) seconds=1
end

1
2
3
4
5
6
7
8

  84.795 ns (0 allocations: 0 bytes)
  62.230 ns (0 allocations: 0 bytes)

1
2

因为如果正确预测分支，它们会非常快，所以由错误检查引起的高度可预测的分支对性能没有太大影响，假设代码基本上从不出错。因此，像边界检查这样的分支非常快。如果绝对最大性能至关重要，或者如果边界检查发生在会进行 SIMD 向量化的循环中，您应该删除边界检查。

如果分支不容易预测，通常可以重新编写函数以避免分支全部在一起。例如，在 copy_odds! 上面的例子，我们可以这样写：

function copy_odds_branchless!(dst::Vector{T}, src::Vector{T}) where {T <: Integer}
    write_index = 1
    @inbounds for i in eachindex(src)
        v = src[i]
        dst[write_index] = v
        write_index += isodd(v)
    end
    return dst
end;

1
2
3
4
5
6
7
8
9

with_terminal() do
  dst = rand(UInt32, 2^18)
  src_random = rand(UInt32, 2^18)
  src_all_odd = [(2*i+1) % UInt32 for i in src_random]
  @btime copy_odds_branchless!($dst, $src_random) seconds=1
  @btime copy_odds_branchless!($dst, $src_all_odd) seconds=1
end

1
2
3
4
5
6
7

  117.822 μs (0 allocations: 0 bytes)
  119.685 μs (0 allocations: 0 bytes)

1
2

除了由循环本身引起的分支（这很容易预测）之外，它不包含其他分支，并且导致速度比完美预测的略差，但对于随机数据要好得多。

当可以使用其他指令完成相同的计算时，编译器通常会删除代码中的分支。当编译器无法这样做时，Julia 会提供 ifelse 函数，该函数有时可以帮助消除分支。

# 注意内存依赖

更深入地思考一下，为什么上面完美预测的示例比完全避免在那里有额外分支的解决方案更快？

让我们看一下汇编代码。在这里，我省略了循环的汇编（因为几乎所有的时间都花在了那里）

对于分支完整版本，我们有：、

1 L48:
2   incq    %rsi
3   cmpq    %rsi, %r9
4   je  L75
5 L56:
6   movq    (%rdx,%rsi,8), %rcx
7   testb   $1, %cl
8   je  L48
9   movq    %rcx, -8(%r8,%rdi,8)
10  incq    %rdi
11  jmp L48

1
2
3
4
5
6
7
8
9
10
11

对于无分支的，我们有：

1 L48:
2   movq    (%r9,%rcx,8), %rdx
3   incq    %rcx
4   movq    %rdx, -8(%rsi,%rdi,8)
5   andl    $1, %edx
6   addq    %rdx, %rdi
7   cmpq    %rcx, %r8
8   jne L48

1
2
3
4
5
6
7
8

全分支每次迭代执行 9 条指令（记住，所有迭代都有奇数），而无分支只执行 7 条指令。查看指令花费多长时间的表格，您会发现所有这些指令都很快。

要了解发生了什么，我们需要更深入地了解 CPU。事实上，CPU 并不像汇编代码让您相信的那样以线性方式执行 CPU 指令。相反，更准确（但仍然简化）的图片如下：

CPU 读入 CPU 指令。然后它会即时将这些 CPU 指令转换为一组更低级别的指令，称为微操作或 µops。µops 和 CPU 指令之间的重要区别在于，虽然指令只能引用几个不同的寄存器，但实际的处理器有更多的寄存器，可以通过 µops 寻址。用 µops 编写的代码称为微码。
这个微码被加载到一个称为重新排序缓冲区的内部数组中进行存储。一个 CPU 一次可以在重新排序缓冲区中保存 200 多条指令。这种存储的目的是允许以高度并行的方式执行微代码。然后将代码批量发送给执行。
最后，来自重新排序缓冲区的结果以正确的顺序发送到内存中。

重新排序缓冲区的存在对您应该如何思考代码有两个重要的影响（我知道）：

首先，您的代码通常以大块并行执行，不一定按照加载顺序。因此，具有更多、更慢 CPU 指令的程序可能比具有更少、更快指令的程序更快，如果前者程序设法更多的并行执行。

其次，分支预测（如上一节所述）不仅仅针对即将到来的分支发生，而是同时针对大量未来分支发生。

当可视化上面的小 copy_odds_branches! 循环的代码是如何执行时，你可以想象分支预测器预测所有分支，比如循环的 6 次迭代，将所有 6 次未来迭代的微码加载到重新排序缓冲区中，执行它们全部并行，然后验证 - 仍然并行 - 其分支是否被正确猜测。

顺便说一句，这种批量处理是分支错误预测对性能如此不利的原因-如果结果证明分支被错误预测，则必须取消重新排序缓冲区中的所有工作，并且 CPU 必须重新开始获取新指令，编译微码，填充缓冲区等等。

让我们考虑一下重新排序缓冲区的含义。除了创建难以预测的分支之外，还有哪些代码可以重写，从而弄乱 CPU 的工作流程？

如果我们这样做呢？

function read_indices(dst::Vector{T}, src::Vector{T}) where {T <: Integer}
  i = 1
  while i ≤ lastindex(src) - 1
    i = src[i] + 1
    dst[i] = i
  end
  return dst
end;

1
2
3
4
5
6
7
8

with_terminal() do
  dst = rand(UInt32, 2^18)
  src = UInt32.(eachindex(dst))
  @btime read_indices($dst, $src) seconds=1
end

1
2
3
4
5

  379.334 μs (0 allocations: 0 bytes)

如果你仔细想想， read_indices 的工作量比任何 copy_odds 函数都要少。它甚至不检查它复制的数字是否是奇数。然而，它比 copy_odds_branches 慢三倍多！

不同之处在于内存依赖性。我们人类，看到输入数据只是一个数字范围，可以准确地告诉函数在每次迭代中应该做什么：简单地复制下一个数字。但是编译器无法预测它加载的下一个数字是什么，因此它需要存储加载的数字。我们说代码对它从 src 加载的数字有内存依赖性。

在这种情况下，重新排序缓冲区没有用。所有指令都被加载，但只是在重新排序缓冲区中保持空闲，因为它们在“轮到他们”之前无法执行。

回到最初的例子，这就是完美预测的 copy_odds_branches! 比 code_odds_branchless! 表现更好的原因。尽管后者的指令较少，但它具有内存依赖性：存储奇数的 dst 索引取决于最后一次循环迭代。因此，与前一个函数相比，一次可以执行更少的指令，其中分支预测器提前预测多次迭代并允许多次迭代的并行计算。

# 可变时钟速度

为低功耗优化的现代笔记本电脑 CPU 在小如邮票（比人的头发还细）的芯片上消耗大约 25 瓦的功率。如果没有适当的冷却，这将导致 CPU 的温度飙升并熔化芯片的塑料，从而破坏它。通常，CPU 的最高工作温度约为 100 摄氏度。功耗以及因此产生的热量取决于时钟速度的许多因素，较高的时钟速度会产生更多热量。

现代 CPU 能够根据 CPU 温度调整其时钟速度，以防止芯片自毁。通常，CPU 温度将成为 CPU 运行速度的限制因素。在这些情况下，更好的计算机物理冷却直接转化为更快的 CPU。旧电脑通常只需清除内部灰尘，更换散热风扇和 CPU 散热膏 (opens new window)即可焕然一新！

作为一名程序员，您无法将 CPU 温度考虑在内，但了解这一点还是有好处的。特别是，CPU 温度的变化通常可以解释观察到的性能差异：

CPU 通常在工作负载开始时工作得最快，然后在达到最高温度时性能下降
SIMD 指令通常比普通指令需要更多的功率，产生更多的热量，并降低时钟频率。这可以抵消 SIMD 的一些性能提升，但 SIMD 在适用时几乎总是更有效。一个例外是相对较新的 512 位 SIMD 指令。在当前（2021 年）CPU 中，这些指令消耗了如此多的功率，以至于由此产生的时钟频率降低实际上会导致某些工作负载的整体性能下降。这个问题可能会在不久的将来得到解决，或者通过降低功耗，通过消费芯片放弃 512 位 SIMD，或者通过编译器拒绝编译这些指令。

# 多线程

在过去的糟糕日子里，随着新处理器投放市场，CPU 时钟速度每年都会增加。部分原因是由于发热，一旦 CPU 达到 3 GHz 大关，这种加速就会变慢。现在我们只看到每代处理器都有微小的时钟速度增量。重点已经转移到每个时钟周期完成更多的计算上，而不是原始的执行速度。 CPU 缓存、CPU 流水线（即整个重新排序缓冲区“工作流程”）、分支预测和 SIMD 指令都是该领域的重要贡献，此处已全部介绍。

CPU 改进的另一个重要领域是数量：几乎所有 CPU 芯片都包含多个较小的 CPU，或者说其中的内核。每个内核都有自己的小型 CPU 缓存，并行执行计算。此外，许多 CPU 具有称为超线程的功能，其中两个线程（即指令流）能够在每个内核上运行。这个想法是，每当一个进程停止时（例如，因为它遇到缓存未命中或分支预测错误），另一个进程可以在同一个内核上继续。 CPU“假装”有两倍数量的处理器。

仅当您的线程有时无法工作时，超线程才真正重要。除了像缓存未命中这样的 CPU 内部原因外，线程也可以暂停，因为它正在等待外部资源（如网络服务器）或来自磁盘的数据。如果您正在编写一个程序，其中一些线程花费大量时间空闲，则其他线程可以使用内核，并且超线程可以显示其价值。

让我们看看我们第一个运行的并行程序。首先，我们需要确保 Julia 实际上是以正确数量的线程启动的。为此，请使用 -t 选项启动 Julia - 例如 -t 8 表示 8 个线程。我已将 Julia 设置为 4 个线程：

Threads.nthreads()

# Spend about half the time waiting, half time computing
function half_asleep(start::Bool)
    a, b = 1, 0
    for iteration in 1:5
        start && sleep(0.06)
        for i in 1:100000000
            a, b = a + b, a
        end
        start || sleep(0.06)
    end
    return a
end;

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

function parallel_sleep(n_jobs)
    jobs = []
    for job in 1:n_jobs
        push!(jobs, Threads.@spawn half_asleep(isodd(job)))
    end
    return sum(fetch, jobs)
end;

1
2
3
4
5
6
7

with_terminal() do
  parallel_sleep(1); # run once to compile it
  for njobs in (1, 4, 8, 16, 32)
    @time parallel_sleep(njobs);
  end
end

1
2
3
4
5
6

with_terminal() do
  parallel_sleep(1); # run once to compile it
  for njobs in (1, 4, 8, 16, 32)
    @time parallel_sleep(njobs);
  end
end

1
2
3
4
5
6

您可以看到，通过此任务，我的计算机可以并行运行 8 个作业，几乎与运行 1 个作业的速度一样快但是 16 个作业需要更长的时间。这是因为 4 个可以同时运行，另外 4 个可以休眠，总共 8 个并发作业。

对于 CPU 受限的程序，核心只忙于一个线程，作为程序员利用超线程并没有什么可做的。实际上，对于最优化的程序，禁用超线程通常会带来更好的性能。然而，大多数工作负载并没有那么优化，并且可以真正受益于超线程。

# 可并行性

多线程比任何其他优化都更难，应该是程序员最后使用的工具之一。然而，这也是一个十分有影响力的优化。科学计算集群通常包含许多（例如数百个或数千个）CPU，每个 CPU 有数十个 CPU 内核，提供了巨大的潜在速度提升，可供挑选。

有效使用多线程的先决条件是您的计算能够分解为多个可以独立工作的块。幸运的是，大多数计算量大的任务（至少在我的工作领域，生物信息学中）都包含令人尴尬地并行的子问题。这意味着有一种自然而简单的方法可以将其分解为可以独立处理的子问题。例如，如果对 100 个基因需要进行某种独立的计算，那么每个基因使用一个线程是很自然的。问题的大小也很重要。产生（创建）线程和从线程的计算中获取结果涉及少量开销。因此，为了获得回报，每个线程都应该有一个至少需要几微秒才能完成的任务。

让我们举一个令人尴尬的并行小问题的例子。我们想要构建一个 Julia set (opens new window)。 Julia set 以 Gaston Julia 的名字命名，与 Julia 语言无关。 Julia set（通常）是复数的分形集。通过将集合成员的实数和复数分量映射到屏幕的 X 和 Y 像素值，可以生成与分形相关的 LSD 迷幻图像。

我在下面创建的 Julia 集是这样定义的：我们定义了一个函数， $f(z) = z^2 + C$ ，其中 $C$ 是一些常量。然后我们记录了可以应用于任何给定复数之前的次数。迭代次数对应图像中一个像素的亮度。我们简单地对网格中的一系列实部和虚部值重复此操作以创建图像。

首先，让我们看一个非并行解决方案：

begin
  const SHIFT = Complex{Float32}(-0.221, -0.713)

  f(z::Complex) = z^2 + SHIFT

  "Set the brightness of a particular pixel represented by a complex number"
  function mandel(z)
    n = 0
    while ((abs2(z) < 4) & (n < 255))
      n += 1
      z = f(z)
    end
    return n
  end

  "Set brightness of pixels in one column of pixels"
  function fill_column!(M::Matrix, x, real)
    for (y, im) in enumerate(range(-1.0f0, 1.0f0, length=size(M, 1)))
      M[y, x] = mandel(Complex{Float32}(real, im))
    end
  end

  "Create a Julia fractal image"
  function julia_single_threaded()
    M = Matrix{UInt8}(undef, 5000, 5000)
    for (x, real) in enumerate(range(-1.0f0, 1.0f0, length=size(M, 2)))
      fill_column!(M, x, real)
    end
    return M
  end
end;

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

with_terminal(() -> @time julia_single_threaded())

  2.067212 seconds (2 allocations: 23.842 MiB)

这在我的电脑上花了大约 2 秒钟。现在是一个并行的：

begin
  function recursive_fill_columns!(M::Matrix, cols::UnitRange)
    F, L = first(cols), last(cols)
    # If only one column, fill it using fill_column!
    if F == L
      r = range(-1.0f0,1.0f0,length=size(M, 1))[F]
      fill_column!(M, F, r)
    # Else divide the range of columns in two, spawning a new task for each half
    else
      mid = div(L+F,2)
      p = Threads.@spawn recursive_fill_columns!(M, F:mid)
      recursive_fill_columns!(M, mid+1:L)
      wait(p)
    end
  end

  function julia_multi_threaded()
    M = Matrix{UInt8}(undef, 5000, 5000)
    recursive_fill_columns!(M, 1:size(M, 2))
    return M
  end
end;

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

with_terminal(() -> @time julia_multi_threaded())

0.537189 seconds (31.45 k allocations: 26.214 MiB, 0.00% compilation time)

这几乎是快 4 倍！使用 4 个线程，这接近于最佳情况，仅适用于近乎完美的令人尴尬的并行任务。

尽管有巨大收益的潜力，但在我看来，多线程应该是提高性能的最后手段之一，原因有以下三个：

在很多情况下，实现多线程比其他优化方法更难。在所示的示例中，这非常容易。在复杂的工作流程中，它很快就会变得混乱。
多线程会导致难以诊断和不稳定的错误。这些几乎总是与多个线程读取和写入同一内存有关。例如，如果两个线程同时递增一个值为 N 的整数，则两个线程都会从内存中读取 N 并将 N+1 写回内存，其中两次递增的正确结果应该是 N+2！令人气愤的是，这些错误的出现和消失都是不可预测的，因为它们是由不吉利的时间造成的。这些错误当然有解决方案，但这是超出本文档范围的棘手主题。
最后，通过使用多线程来实现性能实际上是通过消耗更多资源来实现性能，而不是从无到有。通常，您需要为使用更多线程而付费，无论是在购买云计算时间时，还是在支付多个 CPU 内核的电力消耗增加时，或者比喻为声称拥有更多用户的 CPU 资源，他们可以在其他地方使用。相比之下，更高效的计算不需要任何成本。

# GPUs

到目前为止，我们只介绍了最重要的一种计算芯片，即 CPU。但是还有许多其他类型的芯片。最常见的替代芯片是图形处理单元或者说 GPU。

如上面使用 Julia 集的示例所示，创建计算机图像的任务通常非常并行且具有极高的并行性。在极限中，每个像素的值是一个独立的任务。这需要具有大量内核的芯片才能有效地完成。由于生成图形是计算机工作的基本组成部分，因此几乎所有商用计算机都包含 GPU。通常，它是集成到主板中的较小芯片（集成显卡，在小型笔记本电脑中很流行）。其他时候，它是一张大而笨重的卡片。

GPU 牺牲了本文档中涵盖的 CPU 的许多花里胡哨的功能，例如专用指令、SIMD 和分支预测。它们的运行频率通常也低于 CPU。这意味着它们的原始计算能力比 CPU 慢很多倍。为了弥补这一点，它们拥有大量内核。例如，高端游戏 GPU NVIDIA RTX 2080Ti 拥有 4,352 个内核。因此，某些任务可以使用 GPU 体验 10 倍甚至 100 倍的加速。最值得注意的是，对于科学应用，矩阵和向量运算是高度可并行化的。

不幸的是，我正在编写本文档的笔记本电脑只有集成显卡，而且还没有一种使用 Julia 与集成显卡接口的稳定方法，因此我无法展示示例。

还有更深奥的芯片，如 TPU（明确设计用于深度学习中常见的低精度张量运算）和 ASIC（用于单一应用的高度专业化芯片的总称）。在撰写本文时，这些芯片不常见、价格昂贵、支持不佳且用途有限，因此对非计算机科学研究人员没有任何兴趣。

jokervTv 爬坑历程

Choose mode