本节书摘来自华章计算机《高性能科学与工程计算》一书中的第3章，第3.3节,作者:(德）Georg Hager Gerhard Wellein 更多章节内容可以访问云栖社区“华章计算机”公众号查看。
3.3 案例分析：Jacobi算法
Jacobi算法是数值分析和模拟中许多基于stencil循环方法的原型。在其最简单的形式中，可以用来求解一个标量函数Φ (r→, t)的扩散方程：
求解一个长方形区域的狄利克雷边界条件。当使用有限差分法时，其微分算子是离散的（在不丧失一般性的前提下，这里我们限定为二维方程。但请分析习题3.4中，二维和三维性能的区别）。
在每一个时间步长，Φ在坐标(xi, yi)处的修正值δΦ，由式（3-7）使用其四个邻居点的原值计算获得。当然，Φ的新值必须写回第二个数组中。当所有点都完成更新后，算法进入下一层迭代。代码清单3-1是该内核的一个可能实现。该实现解决了稳定状态的问题但缺乏收敛条件，当然这不是我们考虑的重点（注意交换t0和t1区域不需要逐个元素交换。同该算法的初步实现相比，仅通过交换第三个数组的索引，我们就获得了两倍的性能提升）。
许多优化方法都可能加速这段代码。我们首先使用平衡分析预测这段代码的性能，然后和实测值相比较。图3-5说明了一个二维Jacobi算法的五点更新过程。每次更新需要四次数据加载和一次数据存储操作，但是其“下行邻居”（downstream）phi(i+1, k, t0)在两次迭代之后一定会再次用到（从cache中），所以在计算代码平衡值时，只有三次数据读取操作（共有四次）：Bc = 1.0 W/F（如包含写分配，则代码平衡值为1.25 W/F）。然而，如图3-5所示，按行遍历会将k坐标最大（即phi(i,k+1,t0)）的元素第一次加载到cache中。这是在内存传输中不可避免的。如果cache足以容纳超过两行的数据元素，那么在cache中三次连续的行遍历会使用这个元素。在这种情况下，我们可以假设加载k和k-1行没有时间消耗，所以代码平衡值降为Bc = 0.5 W/F（如果包含写分配，其值为0.75 W/F）。如果内层维度逐渐变大（如图3-5所示，列变多），必然要增加一个，最终三个额外的数据加载操作导致代码平衡值回到令人不愉快的Bc = 1.0(1.25) W/F。
代码清单3-1 二维Jacobi算法的简单实现
根据这些平衡因子，我们可以计算出给定硬件架构上代码的lightspeed。3.1.2节已经给出了在英特尔Xeon 5160平台上的STREAM结果。当cache足以容纳两个连续行时，数据传输特点和STREAM COPY以及SCALE相匹配：一次数据加载、一次数据存储，再加上一次强制性写分配。因为Jacobi内核只包含一次乘法操作，而不是三次加法操作，所以理论的机器平衡值Bm = 0.111 W/F不得不进行修改。因此我们使用
" >
基于这个理论值，并假定写分配不能避免，因此：
" >
其中，修正的理论峰值性能为P+max = 12•4/6GFlop/s = 8 GFlop/s，则预测性能为1.78GFlop/s。然而，基于表3-3列出的STREAM COPY值，该预测性能应再乘以0.38，实际预测性能为675 MFlop/s。随着内层维度的增大，cache已经不能够同时容纳两个甚至一个连续行，代码平衡值第一次上升为Bc = 1.0 W/F，最终为Bc = 1.25 W/F。图3-6 显示了附加各种限制和预测的、不同内层维度的性能对比图（为节省内存，当N较大时（即kmax<图3-6同时也引入一个新的性能指标：更新的区域数目（LUP） /秒。因为它强调“工作完成”而不是MFlop/s，所以更适合stencil算法。在我们的案例分析中，Flop和LUP间有一个简单的1 : 4的对应关系。但在一般情况下，MFlop/s会随着算法优化方法的使用而变化，例如使用不同编译器重新排列代码序列而导致每次更新所要执行的浮点数运算数目不同。然而，用户最感兴趣的是在一定时间内可以完成多少实际工作。LUP/s使得当底层问题相同时，不管采用何种优化方法，所有的性能都具有可比性。例如，许多处理器提供了Fused Multiply-Add（FMA）机器指令执行r = a + b•c操作：一次可执行两次浮点数运算。在这种情况下，因为每次浮点数运算延迟的减少，FMA可大幅提高性能。将代码清单3-1 Jacobi算法代码重写如下：
这个版本用7次浮点数运算代替了4次浮点数运算；当算法是访存受限时，MLUP/s性能指标并没有发生变化（这留给读者使用平衡分析方法证明），但是MFlop/s却发生了变化。
3451
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
