2026a

# 矩阵加法的并行化

以矩阵加法为示例讨论 Julia 中的串行和多种并行计算方式。

阅读提示

这是经典范式中的第一个综合示例。它分别包含了多线程、多进程以及 CUDA 编程的示例，每个示例相对独立，你可以选择性地跳过你不感兴趣的内容。

当你在阅读具体示例之前，建议先阅读快速入门(多线程快速入门、多进程快速入门、或 GPU 计算快速入门 )来了解相关的内容背景。

多线程环境

以下示例需要配置多线程环境，演示在 4 线程环境下运行，注意需要通过 using Base.Threads 加载多线程模块。

本例中需要加载的包：

using Base.Threads
using BenchmarkTools
using Distributed
using SharedArrays

# 背景介绍

当需要将两个矩阵相加时，需要确保这两个矩阵的大小相同。然后，可以通过将这两个矩阵中对应位置的元素相加来得到结果矩阵。

具体来说，如果有两个矩阵 A 和 B,它们的和 C 可以通过以下方式计算：

其中和的范围是 1 到 A 的行数和列数。

这个计算相对简单，且每个迭代计算相互独立，所以作为第一个示例。

接下来会从 Julia 的串行计算开始，测试标准库和朴素串行算法的性能，并在此基础上进行一些性能优化，然后使用多线程，多进程和 GPU 计算进行矩阵加法。

# Julia 串行实现

在 Julia 标准库中有矩阵加法简单的实现：

n = 1024;
A = rand(Int, (n, n));
B = rand(Int, (n, n));
C = A + B;

测试标准库计算使用时间：

function matrixadd_lib(A, B)
    return A + B
end
@btime matrixadd_lib($A, $B);

  1.315 ms (2 allocations: 8.00 MiB)

可以自行写一个朴素的串行函数来计算 n 阶矩阵加法，测试这个朴素实现的性能：

function matrixadd_serial(A, B)
    if size(A) != size(B)# 判断矩阵尺寸
        throw(ArgumentError("A, B尺寸不匹配"))
    end
    C = zeros(eltype(A), size(A))# 预分配结果矩阵
    for i in axes(C, 1)
        for j in axes(C, 2)
            C[i, j] = A[i, j] + B[i, j]
        end
    end
    return C
end

与内置的加法函数进行数值一致性比较：

matrixadd_serial(A, B) ≈ A + B

true

测试串行朴素计算时间：

@btime matrixadd_serial($A, $B);

  16.204 ms (2 allocations: 8.00 MiB)

这个计算比标准库要慢是正常的，先在此基础上做一些串行的性能优化。

由于 C 的每个元素都会被赋值，所以可以使用undef来预分配一个未初始化值的数组，这样性能上会有一定提升。

另外考虑到 Julia 的内存分布是列优先的，所以交换循环顺序。让数据沿内存读取顺序迭代。另外由于做了矩阵尺寸检验，考虑用@inbounds跳过内存边界检查来提升性能：

function matrixadd_serial_opt(A, B)
    if size(A) != size(B)# 判断矩阵尺寸
        throw(ArgumentError("A, B尺寸不匹配"))
    end
    C = Matrix{eltype(A)}(undef,size(A))# undef定义预分配数组
    for j in axes(C, 2)
        for i in axes(C, 1)
            # 添加@inbounds宏
            @inbounds C[i, j] = A[i, j] + B[i, j]
        end
    end
    return C
end
@btime matrixadd_serial_opt($A, $B);

  1.182 ms (2 allocations: 8.00 MiB)

接下来考虑在此基础上进行多线程，多进程，GPU 的计算。

# 多线程计算

考虑在单机上用共享内存多线程并行化矩阵加法。

需要以多线程的方式启动 Julia，用 julia -t 4 或 julia --threads 4 启动 4 个线程的 Julia。Julia 程序不使用多线程计算的话会默认运行在主线程上。

对于 for 循环，Threads 库提供了简单的宏@threads来并行化它：

function matrixadd_threads(A, B)
    if size(A) != size(B)
        throw(ArgumentError("A, B尺寸不匹配"))
    end
    C = Matrix{eltype(A)}(undef, size(A))
    @threads for j in axes(C, 2)
        for i in axes(C, 1)
            @inbounds C[i, j] = A[i, j] + B[i, j]
        end
    end
    return C
end
@btime matrixadd_threads($A, $B);

  488.200 μs (22 allocations: 8.00 MiB)

可以看到，对于这段代码简单地添加一个宏符号@threads，就能利用多线程计算极大地提升代码速度。@threads 是一个用于并行执行 for 循环的宏，将 for 循环根据线程数进行均匀划分并行执行。

但线程的创建销毁和调度分配是有开销的，故时间并非是 4 线程应对应的 4 倍加速。

# 多进程计算

现在考虑用分布式内存多进程并行化计算矩阵加法。

分布式内存多进程方式可以用在大型机器集群中，也可以在单台计算机上使用，这样可以充分利用 CPU 的多核心，或是工作站的多路 CPU，且不需要担心多线程中因共享内存存在的数据竞争问题。

Julia 中的分布式内存多进程计算使用的是 Distribute 库。

在只含有主进程的 Julia 中使用 Distributed 库中的addprocs添加 4 个工作进程：

addprocs(4)
procs()

5-element Vector{Int64}:
 1
 2
 3
 4
 5

Julia 进程之间的通信通过 Julia REPL 主进程调度，该进程接受用户程序输入并控制所有其他进程。我们也可以在启动 Julia 时通过Julai -p n在本地主机上创建 n 个额外的进程。procs() 可以用于查看进程数，返回所有可用的 Julia 进程的 id。

看看用法简单的 @distributed，需要注意，因为 Distributed 自动发送的是数据拷贝。如果不使用SharedArray创建共享数组，@distributed的 for 循环内只会使用拷贝数据进行计算，所以对于赋值语句，更改的只是拷贝数组，最终函数只会返回初始化的 C 数组。

function matrixadd_distributed(A, B)
    if size(A) != size(B)
        throw(ArgumentError("A, B尺寸不匹配"))
    end
    C = SharedArray{eltype(A)}(size(A))
    @sync @distributed for j in axes(C, 2)
        for i in axes(C, 1)
            @inbounds C[i, j] = A[i, j] + B[i, j]
        end
    end
    return C
end

测试数值正确性：

matrixadd_distributed(A, B) ≈ A + B

true

检测性能：

@btime matrixadd_distributed($A, $B);

  25.814 ms (957 allocations: 40.31 KiB)

相比朴素实现基本没有进步，为什么会这样？

由于分布式内存的性质，多进程间的数据通信和同步是重量级的，需要很大开销。这是导致分布式计算没有太多性能改善的主要原因。

可以看到，在这个例子中多进程并行计算提升的资源利用率并不能抵消进程调度和进程间通信带来的开销。

# GPU 计算

CUDA 支持

GPU 计算需要 Nvidia GPU 硬件设备支持，用户需要提前手动安装 CUDA.jl 包并using CUDA

测试对比 GPU 上进行加法计算的性能：

A_d = CUDA.rand(Int, (n, n))
B_d = CUDA.rand(Int, (n, n))
function matrixadd_GPU(A, B)
    CUDA.@sync A + B
end
@btime matrixadd_GPU(A_d, B_d);

  38.400 μs (34 allocations: 2.05 KiB)

不用做太多的程序改写，就可以得到成倍的加速。

只有CuArray数据类型才能在 GPU 上进行计算，故上方代码使用CUDA.rand为 GPU 计算创建矩阵，而非使用rand创建。

我们也可以用CuArray进行类型转换，将数据在 CPU 上生成，然后使用CuArray发送至 GPU 上。

n = 1024
A_d = CuArray(rand(Int, (n, n)))
B_d = CuArray(rand(Int, (n, n)))

不过这种数据转移是相对于直接在 GPU 上生成矩阵性能开销是更大的。

CUDA.@sync用于同步 GPU 上的计算，如果不添加CUDA.@sync的话，测量时间并非是计算完成所需时间而是开始计算需要的时间，用户需要注意这一点。

# 对比与方案选择

我们得到以下时间对比：

`@btime` 测试	1024x1024 `Int`矩阵加法时间
`matrixadd_lib`内置函数 CPU 计算	1.315 ms
`matrixadd_GPU`内置函数 GPU 计算	38.400 μs
`matrixadd_serial`函数朴素单线程 CPU 计算	16.204 ms
`matrixadd_serial_opt`函数单线程优化后 CPU 计算	1.182 ms
`matrixadd_threads` 函数 4 线程 CPU 计算	488.200 μs
`matrixadd_distributed`函数 4 进程 CPU 计算	25.814 ms