创新计算机体系结构设计的 FMM 算法分析

1. 创新计算机体系结构设计的FMM算法分析 • 吕超 • 上海交通大学 • 软件学院 • 2010.6

2. 内容提要 • 课题背景 • 前期工作 • N-body 问题简介 • FMM 算法分析 • 为 FMM 优化的体系结构配置策略 • 结论

3. 课题背景 • 课题来源：新概念高效能计算机体系结构及系统研究开发 • 国家863 计划重点项目（2009AA012201） • 上海市科委重大科技攻关项目（08dz501600） • 项目内容：新型体系结构设计的应用分析及前端设计 • 前期应用分析 • 体系结构的前端设计 • 编译器/软件平台设计 • 应用优化 • 主要目标：设计针对高性能计算的可重构专用处理机体系结构

4. 前期工作 • 高性能计算应用分析 • CT 和 MRI 的图像重建 • 基于SURF算法的图像局部特征提取与匹配 • 应用的模拟与优化 • 基于多核 CPU 并行的 SURF 算法优化与分析 • 基于 GPGPU 的 SURF 算法实现（CUDA-SURF） • 基于 CPU 和 GPU 异构平台的 SURF 算法优化

5. N-body 问题简介 • 引入目的 • 挑选体系结构的典型应用并加以分析 • 给出针对 N-body 特定应用的体系结构优化策略 • N-body 问题 • 又称多体问题，是天体物理学、流体力学以及分子动力学的基本问题之一 • 用来模拟一个系统中相互作用的粒子的运动规律 • 高性能计算的典型应用 • 数学意义：一组已知初始值的常微分方程 （1）

6. N-body 问题简介（续） • 常见算法 • PP（Particle to Particle）算法 • 应用公式直接计算 • 时间复杂度 O(N2) • PM（Particle Mesh Method）算法 • 利用粒子网格，将多个点的作用看作整体（计算网格的势能） • 时间复杂度 O(NlogN)，近似算法，只适用于长距离间粒子的相互作用 • TM（Tree Method）算法 • 利用树形结构划分远程区域与近程区域，再分别计算 • BH（ Barnes-Hut ）算法，时间复杂度 O(NlogN)，误差在 1% 的数量级 • FMM （Fast Multipole Method）算法，BH 的改进算法，时间复杂度 O(N)，精度可控

7. FMM 算法分析 • FMM 算法的数学意义 • 相关数据结构 • 算法流程分析

8. FMM 算法的数学意义 • 适用于 FMM 算法的公式都满足以下形式，如公式（1）： • 使用树形结构（四叉树、八叉树）对空间进行划分，从而得到： • 对 Kfar函数做展开并截断，得到一近似结果： （2） （3） （4）

9. FMM 算法的数学意义（续） • 展开式 • 粒子的作用可以汇聚起来，并由两种形式表示，即 Multipole Expansion（ME）和 Local Expansion（LE） • 它们都是形如公式（4）的泰勒级数，但是 ME 和 LE 收敛于空间中的不同区域 • 计算方法 • 只需要分步骤的计算 ME、LE 和它们之间的转换 • 每个步骤的时间复杂度都是O(N) • 如图所示：

10. 第1层 第2层 第3层 相关数据结构 • 树形结构 • 平均划分各层，划分数（l表示层数， d 表示维度，二维平面划分如左图）： 2ld • 使用四叉树或者八叉树结构存储划分（如右图，以二叉树表示） • 顺序存储各个节点，每个节点包括其中心坐标、ME 和 LE

11. 相关数据结构（续） • 邻居列表与本地列表 • 邻居列表（a）：同处一层并与盒子 b 共享边界的所有盒子叫做盒子 b 的邻居 • 本地列表（b）：盒子 b 的邻居盒子加上盒子 b 自身 • 交互列表 • 对于位于第i层的盒子b，它的交互列表是指其父节点（第i-1层）所有邻居节点（第i-1层）的子节点（第I 层）所对应的盒子，在除去盒子b 所有邻居之后，剩下的所有盒子。

12. 自底向上运算 自顶向下运算 ME到远程节点LE的转换 ME自底向上传递 LE自顶向下传递 计算Multipole Expansion 计算Local Expansion P2M M2M M2L L2L L2P 算法流程分析 • 五个 Kernel • 自底向上运算：P2M（Particle to Multipole Expansion）， M2M（Multipole Expansion to Multipole Expansion） • 自顶向下运算：M2L（Multipole Expansion to Local Expansion）， L2L（Local Expansion to Local Expansion）， L2P（Local Expansion to Particle）

13. 为 FMM 优化的体系结构配置策略 • PetFMM 应用 • PetFMM 应用的计算与访存 • 优化的体系结构配置策略

14. PetFMM 应用 • 涡流质点方法（Vortex Particle Method） • 在流体力学中，可以把流体看作是由大量的粒子所组成，通过分别计算各粒子的速度，来得到流体的速度[11] • 可看作是 N-body 问题 • PetFMM 应用 • 涡流质点方法计算中的一部分 • 计算在二维空间中高斯分布下的涡流粒子按照毕奥－萨伐尔定律发生的速度变化[12] • 由C++实现，以 PETSc Library 为基础设计了相关的数据结构和算法，并初步实现了CPU上的并行。

15. PetFMM 应用的计算与访存 • 计算量分析 • 主要为双精度加法与乘法 • 少量双精度除法与自然指数运算 • 各阶段计算量分析见下表： 展开式截断值 p，四叉树高度为l，粒子总数为N，底层盒子中平均粒子数NP，交互列表平均盒子数NIL，本地列表平均盒子数NLL，盒子总数为NB

16. PetFMM 应用的计算与访存（续） • 各数据结构存储量（Bytes） • 分阶段访存量（Bytes） 展开式截断值 p，四叉树高度为l，粒子总数为N，底层盒子中平均粒子数NP，交互列表平均盒子数NIL，本地列表平均盒子数NLL，盒子总数为NB

17. PetFMM 应用的计算与访存（续） • 访存计算比 • 为计算单元与访存单元的比例配置提供依据 • 计算规模在百亿级别以上，即N>=1010，为此需要选择合适的四叉树深度。这里选择截断值16和18，深度14和16，粒子数1010和1011一共八种情况进行讨论 最大比例出现在 L2P 和 P2P 阶段，约为 0.5

18. 优化的体系结构配置策略 • 设想中的可重构专用处理器体系结构 • 由一颗通用处理器和多块可重构的 SIMD 的 PE 阵列组成 • 45/65 nm 制程，512个双精浮点乘加单元，DDR2/DDR3 内存控制器

19. 优化的体系结构配置策略（续） • 计算单元配置 • 双精度浮点乘加单元 FMA 负责各 Kernel 的运算 • 通用处理器负责树结构和列表结构的初始化运算 • 加入 SFU（Special Function Unit），以加速自然指数运算等复杂运算 • 访存单元配置 • 假设理论带宽为10.667 GB/s 的 DDR3-1333 内存，以最大访存计算比 0.5 来计算，对于一个运行在 500MHz，有 512 个双精度浮点运算单元的处理器来说，需要至少 10 个独立的访存部件才能保证处理器在该阶段达到性能峰值。 • 无临时数据复用的极端情况下会成为性能瓶颈 • 片上缓存配置 • 可编程的片上缓存，适应各 Kernel 规则的访存 • 优化 FMM 的迭代次序，提高临时数据的复用率

20. 结论 • N-body问题中的 FMM 算法更适合作为典型应用 • O(N) 的时间复杂度 • 精度可控 • 针对 FMM 算法优化的体系结构 • 以双精度乘加单元为主，个别特殊运算单元 • 数据结构的初始化工作由通用处理器完成 • Kernel 规则，SIMD 方式运算 • 可编程的片上缓存，并以此对FMM的访存次序进行优化 • 访存单元在个别 Kernel 中可能成为瓶颈

21. 未来工作展望 • 综合考虑多个应用确定体系结构的最终配置 • FFT • LINPACK • 依据最终配置，对FMM算法进行一定的优化 • 优化策略的验证与模拟 • 在多核 CPU 上的模拟 • 在 GPGPU 上的模拟

22. Any question?

23. The End Thank you