本博客的目的是(作者：Phil Thierry 和 Leo Borges —  英特尔 SSG 能源工程团队成员，石油天然气技术领域)初步讨论如何研究三维地震波传播，具体地说，就是在全混合模式下研究一个、两个或多个英特尔® 至强融合™ 处理器的逆时偏移(RTM)算法，也就是充分利用 CPU 内核。

首先我们围绕主题提出一些问题，这些问题对于我们所要研究的内容以及我们如何从软件角度寻找解决方案十分重要。   除了理论探讨，我们还将分享对于有限差分模板的观点，并给出一些代码示例。

 

引言

有限差分模板这一课题十分庞大，实施方案几乎是无限的。随着更多的观众关注这些主题以及新课程的提出，我们将继续撰写博客，同时我们还希望这一讨论能够扩展到具体架构之外的主题。   我们思考的问题列述如下。您对哪些感兴趣？

1.     波方程近似法对程序行为有何影响？包括各向同性的恒定密度、可变密度、垂直横向各向同性、倾斜横向各向同性，弹性波等等。

2.     我们为快照考虑哪种 i/o 方案？我们真的需要显式 i/o 吗？

3.     如果我们考虑在英特尔至强融合协处理器上向前/向后传播，我们如何处理快照？内存中的检查点未必是一个解决方案，但是如何实现快速异步 i/o？我们需要与 PCIe 的本地存储连接，还是 PCIe/IB 上的全局存储？

4.     如何选择时间或频率域：从英特尔® 至强融合™ 协处理器的角度，关于数值法应该考虑什么？

5.     我们应当使用域分解吗？当我们考虑弹性建模或频率内容增加时，每个节点还能够处理一个或多个指令(shot)吗？这对数值法有何影响？

6.     我们有足够的内存向英特尔至强融合协处理器发送一个或多个指令(shot)以计算整个逆时偏移(RTM)吗？或者，我们需要沿用传统卸载模式吗？

7.     从计算的角度看，模板长度应该有多长？我们能够自动找到可同时满足地理位置和计算机标准的最佳模板顺序吗？

8.     为何要及时执行高阶模板？

9.     除了标准 C、FORTRAN、MPI 和 OMP 之外，还有必要考虑 Cilk Plus 和 TBB 等新编程模式吗？有必要将这里的相似性和对比分析加以扩展吗？

10.  这一全混合异构端口一定需要动态负载均衡以同时处理设备性能与指令(shot)大小的差异吗？

11.  在 CPU 内核上观察到的表现是否会在协处理器上重现？

12.  有（无）矢量化与数据调整对协处理器的影响和对相应 CPU 的影响有何不同？采用什么技术以及在什么场景下影响最大？

13.  协处理器上的先进功能需要（或有必要）硬件支持吗？

14.  我们如何使用特征曲线描述并建立模型研究 FMA 指令的影响？

15.  关于如何最高效地重叠操作计算和数据传送，已知的最佳技术是什么？如何最大程度地利用数据一致性以及协处理器到协处理器的传送？  

16.  我们还需要针对其它地震算法撰写帖子吗？如波场插值偏移、基尔霍夫概率地震需求模型(PsDM)、体层摄影术和水库模拟。

17.  电源管理是个现实问题。TTI 3DFD，第 8 个顺序显示为 1000+ pJ/Ops，而对于每秒百亿亿次浮点运算来说目标是 20 pJ/Ops。当然所需的 60x 因数很大程度上来自硬件，但是让软件尽可能接近可用峰值也很重要。那么我们该怎么办？无内在编程的标准实施可达到全负载节点 CPU 峰值的 40%，在对应协处理器上能达到多少呢？如果要获得这一高性能并保持高级别编程模式的易用性，仅仅依靠编译器行吗？

18.   ….

毫无疑问，我们应当利用本博客向列表添加更多内容。

我们请读者朋友们一起切入正题：

由于 RTM 的时间域实施，我们建议在没有交换或界限条件下启动 FD 内核。

我们可以考虑使用卸载模式模拟每个时间步的空间导数的卸载，或者我们可以使用本机模式模拟“每个协处理器的一个或多个指令(shot)”。

关于波方程近似法，我们可以讨论各向同性情景，该情景受 CPU 带宽限制(不受缓存限制)，或者我们可以直接考虑 TTI 实施，后者不受 CPU 带宽限制(甚至不受缓存限制)。

关于性能评测结果，每秒很多点或很多格对于不同的实施来说并不能说明什么，所以我们推荐采用每秒浮点运算或峰值的百分数(%)作为性能评测指标。该指标的一大优势是只需要每格点的浮点运算数即可评测总体性能。最后合计的每秒浮点运算性能指标还可以保护每点浮点运算的数值秘密。

问答

以下是我们和一家合作伙伴公司关于我们所提问题的对话摘要。他们回答问题，我们点评：

问1：包括各向同性的恒定密度、可变密度、垂直横向各向同性、倾斜横向各向同性，弹性波等等

答1 (合作伙伴)：提高建模目标介质的复杂度(从各向同性恒定密度一直到弹性各向异性) a) 由于额外的地球模型以及更多的波场体量而需要提高内存需求(取决于公式) b) 提高每格点的计算次数 c) 更改计算/内存率（也取决于公式）。这些更改会产生一连串的效果 — 例如：提高内存需求会增加解决特定问题所需的节点数量，从而增加通信带宽需求。选择哪种介质类型最为建模对象最终取决于用户，以及工作所在的地质环境。例如：TTI 各向异性和各向同性相比，成像位置会造成 100 米的误差。

答1.2 (英特尔)：答对了！最重要的是 MUL 和 ADD 之间的平衡，它决定着这些代码的效率级别。我们已经介绍，由于 MUL/ADD 的平衡以及字节/浮点运算的要求较低，弹性(visco-elastic)和 TTI 公式的效率高于各向同性，尤其是当使用缓存限制时。内存占用一定程度上可通过域分解处理以支持生产运行。

英特尔至强融合协处理器的实施还有一个问题：“卸载与否”。一旦我们能够完全重叠操作数据传送，为主机计算的每个域卸载热点就像某些标准的“加速”过程。

问 2：我们为快照考虑哪种 i/o 方案？我们真的需要显式 i/o 吗？

答2 (合作伙伴):我不确定您说的显示 I/O 是什么意思。我想您是在比较：从卸载代码/本机执行访问文件系统，以及向主机返回传送数据并写入磁盘。这要取决于性能。如果操作方便只是速度较慢，也不失为一种好方法。这将取决于传播对 I/O 的性能。

答2.1 (英特尔) 我说的 i/o 是关于我们如何处理快照以完成逆时过程。我们需要回到如何选择最佳 i/o 方案的难题，这是我们面对标准集群时都会遇到的。现在对于一个可共享文件系统、内存容量有限并能够执行异步传送的 IA 设备来说又遇到了一个问题：存储在波场的全部磁盘或部分磁盘中、存储在界限中、或使用检查点技术保存在内存中。现在我们回到关键的问题：“卸载与否”。

问 5：我们应当使用域分解吗？当我们考虑弹性建模或频率内容增加时，每个节点还能够处理一个或多个指令(shot)吗？这对数值法有何影响？

答5 (合作伙伴):对于生产 3D 来说域分解是必需的，这是由于数据集的大小和所需的频率内容(尤其是对于 TTI)。

问 6：我们有足够的内存向英特尔至强融合协处理器发送一个或多个指令(shot)以计算整个逆时偏移(RTM)吗？或者，我们需要沿用传统卸载模式吗？

答6 (合作伙伴):参考答5.  只有低频率和/或 2D 偏移才足够小，可以在一个卡上容纳一个以上。

答5.2 和答6.2 (英特尔):我们需要将“数据集”视为叠前地震输入数据，并将“图片大小”这一 3D 表按照参数与字段的数量成倍放大(速率、密度、不均、波场，……)数据集的大小仅和指令(shot)与接收器之间的空间相关。地震跟踪的数量会不断增长，采集面也呈现这一趋势吗？无论如何，这些输入数据很快就会达到 50 PB。

为了更好地实现频率内容恢复，除了输入信号和接收器的能力，我们将看到间距减小和时间采样，这将导致全弹性传播情况下，3D 字段大小随字段数量的增加而增加。

之前完成了地震学(多亏低频率和更大的空间采样)，难道您不认为可变的网格间距甚至可以在弹性 RTM 的情况下支持节点上的 SMP 实施吗？好……这里我需要向您提供一些数字，呵呵。

处理自己指令(shot)集合的众核可出色地将样本代码运行在动态负载均衡充分的任何主机或英特尔至强融合协处理器上。为了在未来几年内实现这一点，我们能够评估所需内存大小和带宽吗？

这一问题完全不同于 FWI，在后者情况下，每个指令(shot)都有助于每次迭代的梯度升级，梯度需要对任一其它指令(shot)可用。

问 13：协处理器上的先进功能需要（或有必要）硬件支持吗？

答13 (合作伙伴)：是，TTI 存在计算限度，还有大量三角法。

答13.1 (英特尔)：对，但是英特尔至强融合协处理器目前不包含正弦/余弦，只包含指数。我们试图使用特征曲线表示“动态正弦/余弦（有或无 VML）；需要数据传送的预先计算的表；使用指数获得超越正弦/余弦”之间的平衡。TTI RTM 需要带各向异性参数(3D 表)的正弦/余弦。在 WEM 条件下，我们还需要计算大量正弦/余弦以获得相速度，后者不需要 3D 存储从而避免采用完全不同的策略。

问 14：我们如何使用特征曲线描述并建立模型研究 FMA 指令的影响？

答14 (合作伙伴)：解答这一问题最好深入了解算术管线，英特尔最适合解答这一问题。

可使用编译器 12.1 和免费可用的 SDE 进行首次评估。通过模拟运行在 Haswell上的应用，可以生成 AVX2/FMA 指令组合。然后就是对该问题进行形式化，并根据指令计数预测未来性能！

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

有限差分模板 — 第 1 部分

这里我们讨论有限差分模板并请您分享您对这一主题的研究结果与经验。

具体地说，我们希望借助本文抛砖引玉，讨论如何实施多维有限差分模板，重点是核内优化和并行计算。这里使用 ISO-3DFD 实施作为示例。

作为示例的计算是简单、恒定、对称、第 8 个顺序三维有限差分(3DFD)并且及时更新各向同性(ISO)。这一 25 点模板计算可按照图式写成时间步外层循环 t=1..nt，每个迭代更新 n1Xn2Xn3 3-dimensional array next( , , )，使用恒定系数 coeff(0), coeff(1), …, coeff(4) 以及来自 n1Xn2Xn3 arrays prev( , , ) 与 vel( , , )的值。在伪码和简单实施的情况下：

for t=1..nt                             //  nt time steps    for i3=4..n3-4                      //  Third dimension in space     for i2=4..n2-4                    //  Second dimension in space       for i1=4..n1-4                  //  First dimension in space         div = coeff(0)*prev(i1,i2,i3)         for r=1..4                     //  8th order stencil            div += coeff(r)*( prev(i1+r,i2  ,i3  ) + prev(i1−r,i2  ,i3  )    // Space first dimension                             +prev(i1  ,i2+r,i3  ) + prev(i1  ,i2−r,i3  )         // Space second dimension                             +prev(i1  ,i2  ,i3+r) + prev(i1  ,i2  ,i3-r)          // Space third dimension         done         next(i1,i2,i3) = 2*prev(i1,i2,i3) - next(i1,i2,i3) + div*vel(i1,i2,i3)    // time update       done     done   done   swap prev <--> next done

这里 n1 代表最快的维度(统一跨度)，n3 代表最慢的维度(跨度程度为 n1*n2)。在(宽松)定义问题以后，让我们选择性能指标。有限差分法的评测指标通常是每次迭代的实际时间(秒)，每秒处理的相邻单元格的数量(i1,i2,i3)(Cell/s)，或每秒浮点运算次数(Flops)。这里我们推荐使用 Flops 指标，因为它可以比较更为复杂的模板方案的不同实施情况。在 ISO-3DFD 示例中，指定时间迭代的每网格浮点运算次数为 7*R+5 = R+3 (乘法)加(6*R+2)，其中 R=4 称为模板半长。每秒的浮点运算次数和处理网格数量之间可以直接转换：Flops = (7*R+5) * Cell/s。

这里第一篇帖子的重点是并行处理工作负载。如需其它背景信息，请点击<这里> 查看示例。在这第一篇帖子中，我们概不讨论关于如何完善有限差分的技巧。

采用并行处理方法是因为我们选择了数据模块化策略：数据在系统可用的内核/线程之间进行分区 — 定义线程模块的大小 n1_Tblock x n2_Tblock x n3_Tblock。我们将 n1 X n2 X n3 域分割为索引列表，描述 n1_Tblock X n2_Tblock X n3_Tblock 数据块。该列表的长度 num_blocks=num_n1_blocks * num_n2_blocks * num_n3_blocks = ceiling((n1-2*R)/n1_Tblock) * ceiling((n2-2*R)/n2_Tblock) * ceiling((n3-2*R)/n3_Tblock)，随后进行填充，确保每个数据块描述符合高速缓存限制：

struct block_struct{   int i1_idx;   int i2_idx;   int i3_idx;   int padding [10]; // 64B cacheline padding };   __declspec(align(64)) block_struct blocking[num_blocks+1];  // one entry per cacheline int index=0; for (int i3b=4; i3b<n3-4; i3b+=n3_Tblock)   for(int i2b=4; i2b<n2-4; i2b+=n2_Tblock)     for(int i1b=4; i1b<n1-4; i1b+=n1_Tblock) {       blocking[index].i1_idx = i1b;       blocking[index].i2_idx = i2b;       blocking[index].i3_idx = i3b;       index++;     }

 

首先请注意：显式数据块列表可支持程序员任意定义线程数据块的列表顺序以及如何分配至线程 — 例如，可以更改循环的顺序。我们这里不探讨重新排序问题。其次，这一框架还可用于不同的线程模式(OpenMP, Cilk Plus, TBB, Pthreads)。这里我们使用 OpenMP:

 

#pragma omp parallel for num_threads(num_threads) schedule(dynamic)       \    firstprivate (n1, n2, n3, num_blocks, n1_Tblock, n2_Tblock, n3_Tblock) \    shared( coeff, prev, next, vel, blocking)    for (int i=0; i<num_blocks; i++) {      int i1b = blocking[i].i1_idx;      int i2b = blocking[i].i2_idx;      int i3b = blocking[i].i3_idx;      apply_stencil(next, prev, vel, coeff, i1b, i2b, i3b, n1, n2, n3, n1_Tblock, n2_Tblock, n3_Tblock);    }

我们使用动态调度确保顺利完成任务(子程序调用)，假设 n*_Tblock 数据块采用小数值时 num_blocks 总是大于 num_threads。

采用正确的数据块大小须经测试。我们这样选择是参考了关于模板循环区块的文献，如 Leopold 撰写的《严格限制三维模板代码的容量缺失》(2002) 以及 Datta 等人撰写的《基于现代微处理器的模板计算的优化与性能建模》(2009)。Leopold 的论文认为，NxNxN 域的矩形区块的形状为 (N-2) x s x (s*L/2)，其中 L 代表高速缓存大小，s 代表模块化因素。在这一区块方案的条件下，我们可以将这一思想加以扩展，并给出模块化因素的初步估计值。模板计算的限制因数是为了确保内存读取偏移量项 prev(i1,i2±r,i3±r)时的数据局部性。为了确保 KNC 架构下使用 512KB L2 高速缓存时的局部性：

n1 * s * s * (16 (每管线 SP) / 2) * 4 (基于上一阵列 SP 字节) <= 512 KB (L2 缓存)

假设统一跨度维度下顺序为 n1=1K 时的阵列大小 n1：

s^2 <= 512K/(1K * 8 * 4)，得出：s <= 4

也就是说，我们区块的初步评估值是 n1_Tblock=n1，n2_Tblock=4 和 n3_Tblock=4*16/2=32。这还意味着在同一跨度方向上没有模块化，从而保持持久连续的内存读取流。首先，在空间第二维度中它似乎显得太“狭小”模块化大小 n2_Tblock=4。不过调用该计算会加载 halo cells prev(i1,i2±r,i3), r<= R=4，所以第二维度中的有效模块化大小是 n2_Tblock+2*R。

请注意在第一版本中：1) 我们不实施任何时间模块化/时间偏移方案。不存在专门利用内部单元属性的代码，或者利用整个问题可单个 SMP 系统中解决的事实。2) 也没有循环级的模块化：区块仅发生在“线程模块化”级别 — 在这一情况下任务通过 blocking[*] 索引列表指派至线程。您也可以考虑多级别的模块化以便更好地控制数据到内核的指派：当线程内部应用实际模板计算时，线程模块化使用稍大一些的数据块以及循环区块。3) 我们未考虑主机到协处理器以及协处理器到主机的传送时间。仅考虑英特尔® 至强融合™ 协处理器上的模板计算时间。

另外，在这第一篇帖子中，我们不做大幅度的代码更改，以便研究有限差分的矢量化。未经具体优化的 C/C++ 版本的子程序的基本示例如下：

// apply 8th order ISO stencil on block [i1b..i1b+n1_Tb]X[i2b..i2b+n2_Tb]X[i3b..i3b+n3_Tb] void apply_stencil(float *ptr_next, float *ptr_prev, float *ptr_vel, float *coeff, // arrays                                             const int i1b,   const int i2b,   const int i3b,                             // block indexes                                      const int n1,    const int n2,    const int n3,                             // full domain                                        const int n1_Tb, const int n2_Tb, const int n3_Tb) {                 // block sizes                         const int n1_end = n1-4, n2_end = n2-4, n3_end = n3-4;   const float c0 = coeff[0], c1=coeff[1], c2=coeff[2], c3=coeff[3], c4=coeff[4];   const int n1n2 = n1*n2;   const int n1_2 = 2*n1, n1n2_2 = 2*n1n2;   const int n1_3 = 3*n1, n1n2_3 = 3*n1n2;   const int n1_4 = 4*n1, n1n2_4 = 4*n1n2;     const int n3b_end = MIN(i3b+n3_Tb, n3_end);   const int n2b_end = MIN(i2b+n2_Tb, n2_end);   const int n1b_end = MIN(i1b+n1_Tb, n1_end);     for (int i3=i3b; i3< n3b_end; i3++) {     float *prev  = &ptr_prev[i3*n1n2+i2b*n1];     float *next  = &ptr_next[i3*n1n2+i2b*n1];     float *vel    = &ptr_vel [i3*n1n2+i2b*n1];     for(int i2=i2b; i2< n2b_end; i2++, prev+=n1, next+=n1, vel+=n1) {   #pragma vector always #pragma ivdep       for(int i1=i1b; i1<n1b_end; i1++) {         float tmp  = c0*  prev[i1]                        + c1*( prev[i1        +1] + prev[i1        -1] +                                  prev[i1      +n1] + prev[i1       -n1] +                                  prev[i1   +n1n2] + prev[i1   -n1n2]  )                        + c2*( prev[i1        +2] + prev[i1        -2] +                                  prev[i1   +n1_2] + prev[i1    -n1_2] +                                  prev[i1+n1n2_2] + prev[i1-n1n2_2]  )                        + c3*( prev[i1        +3] + prev[i1         -3] +                                  prev[i1   +n1_3] + prev[i1    -n1_3] +                                  prev[i1+n1n2_3] + prev[i1-n1n2_3]  )                        + c4*( prev[i1        +4] + prev[i1         -4] +                                  prev[i1   +n1_4] + prev[i1    -n1_4] +                                  prev[i1+n1n2_4] + prev[i1-n1n2_4]  );         next[i1] = 2.0f*prev[i1] -next[i1] +tmp*vel[i1];       }     }   } }

所以这是一个简单实施，也可以使用 Fortran 语言。在此我们不探讨源代码级别的数据一致性属性，也不探讨以矢量化为导向的代码更改，如何提高内存访问效率等等。

在以上测试条件下，我们可以运行我们的第一个实例。我们使用包含 60 个内核的英特尔至强融合协处理器处理 n1=928 x n2=448 x n3=840 问题。每个代码包含三个线程(合计 180 个线程)并且首次测试还将使用 KMP_AFFINITY=balanced。这一点上，测试不跟踪数据传送开销：该实施将任务卸载至协处理器，但是该测试需要在至强融合卡内执行 800 次迭代算法。这样以来就可以忽略与卸载初始化及终止相关的传送时间。只有在第一篇帖子中同时使用 Flops 和 Cell/s 指标报告性能。这里 Cells/s = (n1-2*4)*(n2-2*8)*(n3-2*8)/time-for-one-iteration；Flops = (Cell/s) * (7*4+5)。

首先，在数据到线程的指派过程中，区块大小的估计值可按照以下方法测试：赋值 n1_Tblock=n1 和 n3_Tblock=32，然后应用 n2_Tblock=4 附近值，其中维度应对数据块/区块大小敏感：

 

在这一配置下，ISO-3DFD 模板峰值为 118.4GFlops (3.6GCell/s)且 n2_Tblock=1，同时性能对第二维度模块化因数的不同值相当敏感。在此类简单模板中，在最低维度模块化 n3_Tblock 情况下与第二维度情况下相比，对变量的敏感度也十分低：

 
更复杂的模板可显示更高的性能变量以区分最低维度下的模块化因数。

显然本帖子仅仅是举例说明如何使用区块同时作为模块化因数的启发。并行处理方案本身仅仅是众多方案和可能实施中的一个示例。本示例不拟探讨数据的多级分区：假设并行处理方案将较大数据块指派至每个线程，然后在模板循环计算中实施狭小区块。所有示例均不测试以提高数据局部性为目标的显示线程固定技术。

以后我们会讨论面向数据一致性和矢量化的代码更改。关于实施的带宽分析也是一个有益的话题。

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

有限差分模板 — 第 2 部分

这里我们提供样本源代码供您参考，请点击这里。

ISO-3DFD (版本 1) 原始码应当支持构建并运行卸载测试(make clean; make arch=offload; ./offload.sh) 或原生测试(make clean; make arch=mic; mic_native.sh)。

iso-3dfd_V1/src/iso-3dfd_stencil.cc 是实际模板计算。请注意：这一版本未针对矢量化和数据一致性进行特殊优化。这是模板的简单实施。

iso-3dfd_V1/src/iso-3dfd_parallel.cc 是实际任务并行处理，其中 blocking[*] 阵列被构建同时数据被指派至线程，这一点在上篇博客帖子中已经讨论。

iso-3dfd_V1/src/iso-3dfd_main.cc 是主驱动程序。为了简化代码，驱动程序不进行数据初始化 — 就像随机数据或雷克小波分析原函数。也不进行错误检查。当前版本仅支持单个 MIC 卡：在测试多卡系统时，您可以考虑使用特定的 MIC 卡 ID 编辑卸载语句。

默认值运行 n1=928 x n2=448 x n3=840 问题大小，线程模块化 n1_Tblock=n1=928 (无模块化), n2_Tblock=1 和 n3_Tblock=124。OpenMP 亲和度未在源代码内设置，所以开发人员可通过 KMP_AFFINITY 环境参数测试不同的亲和度方案。驱动程序支持以下命令行参数：

[n1]  [n2]  [n3]  [# threads]  [# iterations]  [thread block n1]  [thread block n2]  [thread block n3]

作为参考，script ./offload.sh 示例说明了如何使用 KMP_AFFINITY= "granularity=thread,balanced" 运行样本代码，以及如何使用包含 60 个内核的英特尔至强融合协处理器支持每代码的 3 个线程：

 

$ make clean; make arch=offload; ./offload.sh

using 180 threads

allocating prev, next and vel: total 3996.56 Mbytes

n1=928 n2=448 n3=840 nreps=800 num_threads=180

n1_Tblock=928 n2_Tblock=1 n3_Tblock=124

-------------------------------

With data transfer IN and OUT

time:          71.47 sec

throughput:  3770.01 MPoints/s

flops:        124.41 GFlops

 

Script ./mic_native.sh 示例说明了如何使用 KMP_AFFINITY= "granularity=thread,balanced" 运行样本代码，以及使用本机的 61 个内核支持每代码的 3 个线程：

 

$ make clean ; make arch=mic ; ./mic_native.sh

allocating prev, next and vel: total 3996.56 Mbytes

n1=928 n2=448 n3=840 nreps=20 num_threads=183

n1_Tblock=928 n2_Tblock=1 n3_Tblock=124

-------------------------------

time:           1.75 sec

throughput:  3850.36 MPoints/s

flops:        127.06 GFlops

 

 

 

期待分享您关于本主题的研究经验和实例！