并行计算基础：围绕矩阵运算理解

这部分整理《并行计算与工业软件》第 3 页。原笔记中的关键词包括 CPU、GPU、cache、NUMA、SIMD、MIMD、数据搬运和任务依赖；下面将它们放回一个简单矩阵计算中理解。

1. 并行计算的对象是什么

对矩阵乘法 \(C=AB\)，每个元素满足

\[ C_{ij}=\sum_{k=1}^{n}A_{ik}B_{kj}. \]

这里有两层并行性：

不同的 \(C_{ij}\) 可以彼此独立地计算；
同一个 \(C_{ij}\) 内部的求和可以做并行归约。

但并行并不等于把循环简单分给更多处理器。程序还必须读取 \(A\)、\(B\)，写回 \(C\)，协调任务，并合并局部结果。真正耗时可以粗略分为

\[ T_{\mathrm{total}} =T_{\mathrm{compute}} +T_{\mathrm{memory}} +T_{\mathrm{communication}} +T_{\mathrm{synchronization}} +T_{\mathrm{overhead}}. \]

这正是原笔记“高效的算法不等于高效的程序”所指的区别。

2. 从外存到计算单元

数据通常沿以下层次移动：

SSD/磁盘 -> DRAM -> L3/L2/L1 cache -> 寄存器 -> 算术单元

越靠近算术单元，容量越小、延迟越低。矩阵乘法的经典分块算法会把一个小块反复留在 cache 中，以提高数据复用率，而不是每次乘加都从主存重新读取。

两个重要指标是：

延迟：一次数据访问从发起到完成所需的时间；
带宽：单位时间最多能传输多少数据。

对大量连续数据，带宽通常更重要；对频繁随机访问，延迟更明显。因而缓存友好的连续访问，往往比跳跃访问更快。

3. CPU、GPU 与加速器

CPU 的少量强核心适合复杂控制流、分支和低延迟任务。GPU 拥有大量相对简单的执行单元，适合对许多数据执行相同运算，例如矩阵乘法、卷积和有限元单元积分。

把计算移到 GPU 时还要考虑主机和设备之间的数据搬运：

\[ T_{\mathrm{GPU}} =T_{\mathrm{H2D}}+T_{\mathrm{kernel}}+T_{\mathrm{D2H}}. \]

若矩阵很小，传输和启动 kernel 的开销可能大于计算收益。只有当数据复用充分、计算量足够大，或数据能长期留在 GPU 上时，加速才更可能出现。

4. SIMD、MIMD 与 Flynn 分类

Flynn 分类按照指令流和数据流区分体系结构：

类型	含义	直观例子
SISD	单指令、单数据	单核上的普通串行程序
SIMD	单指令、多数据	CPU 向量指令、GPU 的成组执行
MISD	多指令、单数据	实际通用计算中较少见
MIMD	多指令、多数据	多核 CPU、计算节点集群

GPU 常被概括为 SIMD，但现代 GPU 更准确地说采用 SIMT：程序写成许多线程，硬件把线程按 warp/wavefront 成组执行。如果同一组线程走不同分支，就会出现分支发散，使部分执行单元暂时闲置。

5. 共享内存、分布式内存与 NUMA

5.1 共享内存

多个核心访问同一地址空间，线程之间传递数据方便，但必须处理竞争、锁、原子操作和缓存一致性。OpenMP 和线程库常用于这一模型。

5.2 NUMA

NUMA（Non-Uniform Memory Access）机器虽然提供统一地址空间，但某个 CPU 访问本地内存通常比访问另一个 CPU 插槽连接的远端内存更快。因此，“数据在哪个插槽首次分配”“线程在哪个核心运行”会影响性能。

5.3 分布式内存

集群中每个节点有自己的内存，节点间通过网络发送消息。MPI 是典型编程接口。此时算法不仅要考虑浮点运算次数，还要考虑消息数量和通信数据量。

6. 数据依赖决定能否并行

若任务 \(P_1\) 读取或修改的数据区域为 \(\Omega_1\)，任务 \(P_2\) 对应 \(\Omega_2\)，当

\[ \Omega_1\cap\Omega_2=\varnothing \]

时，两者通常可以独立执行。若有交叠，就要进一步区分：

只读-只读：一般可安全并行；
读-写：可能读到不同步的数据；
写-写：会产生数据竞争，需要重新划分或同步。

BSP（Bulk Synchronous Parallel）把程序组织为“局部计算 - 通信 - 屏障同步”的多个超级步。它易于推理，但频繁屏障会让快任务等待慢任务，所以负载均衡同样重要。

7. 一个适合笔记本的小实验

先手写三层循环：

def matmul_ijk(a, b):
    n, p, m = len(a), len(b), len(b[0])
    c = [[0.0] * m for _ in range(n)]
    for i in range(n):
        for j in range(m):
            for k in range(p):
                c[i][j] += a[i][k] * b[k][j]
    return c

再把循环次序改成 i-k-j：

def matmul_ikj(a, b):
    n, p, m = len(a), len(b), len(b[0])
    c = [[0.0] * m for _ in range(n)]
    for i in range(n):
        for k in range(p):
            aik = a[i][k]
            for j in range(m):
                c[i][j] += aik * b[k][j]
    return c

两者的浮点运算量都约为 \(2npm\)，但第二种写法更容易连续访问 b[k][j] 和 c[i][j]。这可以直观看到：算法复杂度相同，不代表真实运行时间相同。

最后与 NumPy 的 a @ b 对照。NumPy 通常调用经过分块、向量化和多线程优化的 BLAS 库，因此它代表“成熟实现”，而不是 Python 循环自然并行后的结果。

建议记录：

矩阵规模	`ijk` 时间	`ikj` 时间	NumPy 时间	结果是否一致
\(n=16\)
\(n=64\)
\(n=256\)

小矩阵上不要急于讨论“加速比”，因为解释器、计时器和函数调用开销会占主导。随着规模增大，缓存局部性和底层向量化的差距才会逐渐显现。

8. 如何判断程序受什么限制

算术强度定义为

\[ I=\frac{\text{floating-point operations}}{\text{bytes moved}}. \]

算术强度低的程序更容易受内存带宽限制；算术强度高的程序更可能受计算峰值限制。Roofline 模型用

\[ P\le \min(P_{\mathrm{peak}},\ B_{\mathrm{memory}}I) \]

表达这一上界，其中 \(P\) 是可达到的计算性能，\(P_{\mathrm{peak}}\) 是峰值算力，\(B_{\mathrm{memory}}\) 是内存带宽。

矩阵乘法通过分块反复利用数据，能够提高算术强度；稀疏矩阵向量乘法则常因数据复用少、索引访问不规则而受带宽限制。这也解释了为什么有限元求解中“装配”和“求解”会呈现不同的性能瓶颈。