光子传输的蒙特卡罗方法并行实现.docx

《光子传输的蒙特卡罗方法并行实现.docx》由会员分享，可在线阅读，更多相关《光子传输的蒙特卡罗方法并行实现.docx（21页珍藏版）》请在三一办公上搜索。

1、光子传输的蒙特卡罗方法在GPU上的实现L aszl o Szirmay-Kalos 等翻译：吴涛这篇文章提出一个快速的蒙特卡罗方法来解决激光和gamma射线光子在不均匀 媒介中辐射传输方程问题。激光传输在计算机图表算法中起着重大作用，比如高 能gamma射线就在医疗与物理中扮演重要角色。实时图像的应用在速度上受到 限制，因为我们通常想得到每秒至少20张图来获取运动的动态幻像。快速模拟 在人机交互系统中也十分重要，比如像放射疗法或者物理实验设计，这些地方用 户可以放置一个源，并且希望得到辐射分配的迅速反馈。多重散射模拟的速度在 迭代断层摄像术重建方面也很重要，此时散射辐射要从实际的重建数据中，且

2、要 移除测量部分来估计，而重建仍然继续针对那些假设只代表不散射的部分剩余成 分。17.1光子传输的物理特性计算多重散射以及提出一种在不均匀媒介的实际方法更是困难重重。最精确的逼 近方法是依据蒙特卡罗方法的积分运算以及通过跟踪随机光子的物理模拟现象， 来计算他们的总贡献。光子与媒介粒子在该点碰撞的条件概率密度由消光系数b七（v,分定义，如果 粒子密度不均匀，它则由光子频率v以及位置三来确定。消光系数通常可以用材 料密度C任）和一个仅与频率有关的因子的乘积来表达。光子在碰撞之外还会被反射与吸收。反射的概率称为反照率，用a（v ）表示。 随机反射方向用两个球坐标轴描述：散射角0和方位角甲（图17.1

3、）。源点与新的方向的散射角0的概率密度可能与频率有关故用物理模型来描述 （我们用瑞利散射【12】模型来模拟激光光子，用Klein-Nishina公式【11】【13】 来计算gamma光子）。方位角甲的概率密度鉴于轴向对称分布故而是均衡的。图 17.1光子在媒介中的散射。能量为%的入射光子与材料的粒子碰撞。因为碰撞，光子的方向和 能量将改变。新的方向由散射角9和方位角中指定。新的光子能量E由入射光子能量和散 射角确定。光子能量E = hv正比于它的频率，比值是普朗克常量h。能量（也就是稍高 能量的gamma光子的频率）会因为散射而改变。由康普顿定律，散射角与光子 能量的相对变化有独特的对应关系：

4、E = E-E1 + （1-cos9 ）m c 2e其中，E是散射光子能量，E0是入射光子能量，me是电子静止质量，c是 光速。检查康普顿定律公式，我们得出当E0关于电子能量（mc2）不可忽略时， 能量在散射时的变化很大。但是对于可见光而言，就不再适用 了，此时E0 = h m c ；因而，它们的散射能量将与入射光相似。结果，我们能够按不 同频率分解低能量（激光）辐射模拟，但在高能量（gamma射线），频率将成为 每次模拟光子的特征属性。在计算机模拟中，光子根据目前提到的基本定律来跟踪。模拟过程的输入数 据是媒介的密度，其由三维标量场密度。任）定义，还有仅与频率有关的反照率 函数，以及散射角的

5、概率密度场密度假定是离散数据像三维组织立体像素一样。 仅与频率有关的反照率函数，以及散射角的概率密度函数即可以像表格一样存放 在组织中也可以用代数公式表示。为了模拟光子传输过程，碰撞，吸收和散射方向需要根据物理模型定义的概 率密度来生成。一种方法是根据给定的概率密度生产随机变量，称为反演法。反 演法首先计算累计概率分配（CDF），作为概率密度的积分，然后通过反演单位 间隔均匀分布的变量CDF来生成离散抽样。在形式上，如果一个随机变量的CDF 是P3），则随机抽样x可以通过解方程来产生：P （x） r这里r是一个单位间隔内均匀分布的随机变量，通常可以调用随机数产生器 来获得。17.2光子传输在G

6、PU上的实现光子在空间的传输是相互独立的，因而可以用并行模拟，因此显示出其固有 的物理模型的并行性（自然是大规模并行的机器）。模拟单个粒子路径的算法是 嵌套循环。外循环的主体是负责寻找下一个散射位置，于是抽样吸收，如果还有 残存就抽样新的散射方向。这层循环在光子被吸收或者离开模拟区域时结束。不 同粒子的散射事件发生数目将会在一个大的时间间隔内波动，因而外层循环也会 以同样的次数执行。内层循环一直访问体元直到下一个散射点确定；这个过程取 决于拉伸的随机数同样与沿着光子路径的媒介密度有关。随机数与密度变化，因 此对于不同的抽样，访问的体元数目可能也是不同的。（图17.2）图 17.2在光子传输模拟

7、中控制路径的原因。光子路径有随机数目散射点而且在两个散射点间又有不 同数目的体元。在并行线程中执行不同数目的迭代的问题在GPU中变得起决定性作用。如 果我们将一个光子分配到GPU线程的技术单元，那么在warp中的最长的光子路 径与所有路径中最长的射线将会限制计算的速度。为了保持warp的连续性，线 程执行的长度必须相似，所有的线程优先执行相同的指令序列。为了克服这个困 难，我们根据GPU执行效率的需求调整传输过程的模型；也就是说，我们发展 了能在最少条件下执行的算法。为了解决不同数目的散射事件，我们重复光子步骤，也就是计算线程中的射 线。每一步包括标示另一个散射点，做终止决定，如果光子在碰撞后

8、仍存活则继 续抽样新的方向，反之在源区生成新的光子。这种策略能够保证线程都在忙碌运 行中。经过调整，线程仍可能分歧。确定一个光子是否吸收或者抽样一个新的散射 方向要根据当前的媒介密度；因而，这些任务需要解一个仅含少量变量的方程， 这个可以从实际位置确定。即使因为CDF不能计算或者用符号表示反演而变得 复杂。在这种情况下，数值方法可以得到应用。自由程抽样甚至更差因为它需要 探查媒介的大部分区域，相当于大部分当前的GPU内存通道变得相对缓慢。迭 代数据的根和访问不同密度的区域体元导致线程分歧最终降低GPU的利用率。因为不同光子的跟踪是随机的，它们访问不同模拟体元的部分时，导致数据 存储在不同并行通

9、道缓冲中。这种方式的内存使用使得高速缓存的利用率递降。 为了解决这个问题，也许有人会建议将类似路径的光子分配到一个warp中。因 为随机路径是从均匀分布的随机数中产生的，这可能意味着随机数向量的分类。 一个完整的路径可能联合了许多随机数（比如，超过20个），所有分类操作需要 在一个很高的维度工作，这在实际中是不可行的。所有光子的存储是依据第一个 随机数并且分配到已经分好序列的warp中。这一步提供很少的初始化数据，但 是早晚所有的线程必然能同时访问体元的不同部分。为了保证线程在媒介中追踪射线的连续性，我们发展了一个无条件说明的方 案来抽样吸收和散射方向。为了解决自由程问题同样为了限制数据通道的

10、分歧， 我们提出一项技术，即使用低分辨率的体元数组在除原始的高分辨率的体元数组 的基础上。这项技术显著的减少了到达高分辨率数组的组织数以及发现新的散射 点需要迭代的数目。17.2.1自由程抽样模拟一个光子的当前位置尤，方向3，以及频率V，就是抽样自由程S (就是 个光子在不碰撞下可以移动的距离)，也就是决定沿着射线以s)=x+s的下一个 交互作用点。光子与粒子发生碰撞的概率密度就是消光系数。从这里，自由程 CDF可以用指数衰减模式来表达；因而自由程抽样等于下面关于标量r即单位间 隔均匀分布的量的解：r -)1 - exp ib (V, p(ly)dl = rI )消光系数的积分称为光学深度，从

11、而可以导出下面的方程：T (0, s) = ib (V, P(l)dl = log(1 r)0自由程抽样，也就是这个方程的解，当媒介均匀时形式很简单；换句话说， 当消光系数在媒介中处处相同时：log(1 r) s =b t (V)当媒介不均匀时，消光系数不是常数，但是可以用一个数组表示；这种情况 下，通常是方法是射线形变程即沿着射线每次行进一小步，发现步数n在黎 曼求和逼近光学深度且比log(1 r)稍大：云b (V, p(iAs)As log(1 r) Low-densityrtensiTyConstantdensityregionregion图 17.3虚拟粒子运用。左图为原始实际粒子的自

12、由程抽样(颜色为浅灰)。右图描述了整个体元包 括虚拟粒子(黑色)和实际粒子的情况。与虚拟粒子碰撞不改变光子方向。1. 产生试验性的路径长s和试验性的碰撞位置讯s)，利用上界消光系数函数b (v,分。这需要抽样最大消光的方程式的解：maxT (0, s) = fb (v, P(l )dl = - log(1- r)02. 当一个试验性的碰撞点确定时，我们用概率气(v,分/ajv, X)随机决定散射 发生在实际粒子还是虚拟粒子上。要是虚拟散射发生，那这个粒子的散射方 向不会改变，之后相同的抽样步骤将会从原点再次重复。当一个实际粒子散 射发生时，这层循环结束。注意到作为实际算法的抽样表明，通过原始消

13、光系数CDF产生的，无论抽 样中建立了什么样的上界消光系数。然而，抽样的效率受合适选择的上界系数影 响巨大。一方面，上界系数需要一个针对方程17.2简单的解。另一方面，如果 这个系数与实际系数，即虚拟粒子密度相比，差异大，虚拟散射发生的概率也会 很大，因此需要许多额外的抽样步骤直到实际散射被发现。因此，一个最佳上界 系数给出了方程17.2的简单解并且是与实际消光系数紧接的。要找到一个简单且合理准确的上边界，我们采取分段常数方程，即由低分辨 率如【8】的宏单元网格(图17.4)定义，并且分配原始体元中最大的消光系数 到每个宏单元。要是一个宏单元包括整个体元。那这种方法就和Woodcocktrac

14、king【10】相同。为了获取方程17.2的解，我们执行在宏单元网格三维数字微分分析器的访 问体元方式历经所有被射线影响的宏单元。三维DDA算法以重构为基础，即细 胞的边界在三个并行平面采集，它们垂直于各自的坐标轴。算法保持三个射线参 数代表与平面的交互作用点。其中最小值代表当前出射点。为了到达下一个元点， 稍微增加其值。增加量与三平面采集量有关，其对于给定射线是常数并且只能计 算一次。在遍历间，当射线穿过发现的宏单元边界时，射线参量s,与点一致。我们一 个个访问宏单元直到当前宏单元保持方程17.2的根。(图17.4) 选择宏单元n且包括散射点的不等式为：习6 (v )As - log(1-

15、r) 2L6 (v )Asi=0i=0Real or virtual scattering paint where 4smiri, s) =-log(1-fnd0)其值b max(v)是增加的宏单元i中常数消光系数。Optical depth of the max. densities图 17.4虚拟粒子参与下的自由程抽样。我们在宏单元网格执行三维DDA并且通过分析混合有实虚 粒子的媒介的消光系数b 来获得光学深度t。 max射线行进法与改进的逼近法的重要区别是步长As不是常数，但是通过射线 i与宏单元的作用长度获得的。当步长n在不等式中第一次满足，宏单元的散射点就确定了。由于在确定的宏单元中

16、bmax(v)是常数，它在方程17.2的积分是线性的；因而，自由程等于下 面线性方程的解：log(1- r) + Z1 b i (v )As maxiS = n-1 +-1 LVmax提出的自由程抽样方法比射线行进和Woodcock追踪能显著减少warp的分 歧。在射线行进中，迭代的轮数与自由程长度成正比，变化范围从1到数组线性 解。(如512)在Woodcock追踪中，平均步长与整体最大消光系数成反比，所以 这种算法要经历很多步，并且只能在低分辨率区域找到虚拟粒子碰撞。改进的算 法中，实粒子碰撞的概率密度等于实际系数与当前宏单元最大值之比，所以相关 于分段常数的正确性。增加宏网格的辨析率，通

17、过付出更多的宏单元步骤可以控 制其正确性。举例说明，图17.9中头部数据达到1283个体元，有83个宏网格分 辨率，平均的在找到实粒子碰撞之前寻找的虚粒子碰撞为0.95,然而平均的宏粒 子个数为1.7.如果宏分辨率增加到163，则平均虚粒子碰撞减少为0.92，但是平 均访问的宏粒子增加到3.7.通常情况下，在实际碰撞点很少（至多一到两个虚拟 粒子碰撞）发生，宏单元的可能所用步数远远小于射线行进算法的步数。17.2.2吸收抽样当光子与一个实际粒子碰撞时，我们应该确定吸收是否发生，以及碰撞之后还存 在的粒子的新方向。存活的概率由反照率定义；因而，这个随机方向通过检查反 照率是否大于一个单位间隔随机

18、分布的数：r a （y ）在我们的实现中，与频率相关的反照率函数存储在一维数组里，频率通过发 射光子的最大值y 归一化。应用线性插值方法，反照率能够以插值替换典型的 max频率。对于光线粒子，反照率存储在与对应于频率的红绿蓝光。17.2.3散射方向如果光子碰撞后存活，一个新的继续方向需要抽样，它是由散射角。和方位角甲 决定的。方位角甲是均匀分布的；因此，它用甲=2兀r生成，其中r是一个单位 间隔均匀的随机数。散射角通过解下面的散射方程得到：P （y ,cos 0） = rcos02其中，P （y,cos0）是给定频率散射角余弦的累积概率分配，而r是另一个 cos02均匀分配的随机数。散射角概率

19、密度可以通过物理模型和合适的物理表达式来描 述（我们用瑞利散射模拟光粒子，用Klein-Nishina不同穿越区域用作gamma射 线，其他相函数能够讨论掌握），通过这个方法我们可以积分的求的累积概率分 配。不幸地是，对于大多数模型，积分不能用闭区间表示，抽样方程也不能用解 析解的方式得到。另一方面，数值解，像中点或者错误的位置方法，将执行一个 循环并且速度变慢。要克服这个困难，我们规则的防着抽样点，如r=0,1/N,2/N,1并且v=0， ymax /M,2ymax /M，,ymax,解决这些抽样点（我们假设N和M的值为128）。结果 放在二维数组中，通过y / y max和匕寻址。组织里存

20、储的信息是散射角余弦，相对 于给定频率和随机数（图17.5）。模拟过程中，我们只能通过寻址组织地址获得 散射方向余弦，即通过当前变量均匀随机数和光子频率。因为光子，散射角不受 频率约束。我们能用一维组织存放抽样结果。图 17.5这个表格用来抽样散射角假定是康普顿散射和Klein-Nishina相函数。表格由（0,1）随机变 量分配参数化，并且相对能量是参照电子能量（512kev）的值，提供散射角余弦值。17.2.4并行随机数的生成所有讨论的抽样步骤需要单位间隔均匀分配的随机数。自由程抽样需要与发现第 一个实际散射点所检索的虚拟粒子个数相同的随机数。吸收抽样仅需要一个额外 随机数，但是方向抽样需

21、要两个。因为不同步骤在统计上是独立的，每一步都需 要用一套新的随机数。此外，不同的光子表现的独立，因此，每一个光子的模拟 应该根据一个新的随机数序列。很显然，如果我们在并行线程中用相同的随机数 生成器模拟不同的光子，因而我们没必要多次重复相同的计算。这个问题的标准解是为每个并行线程分配一个私有的随机数生成器。要保证 他们的序列是独立的，私有随机数生成用随机种子初始化。私有随机数生成器在 GPU上运行，并且分配到并行线程中，因此他们的初始化种子用C库函数中的 随机数生成数在CPU上生成。我们注意到，数值分析和构建灵活的随机数生成器是大难题，这是在这张内容之 外的。有兴趣的读者可以参考文献，如CU

22、DA实现【6】，这也在我们的程序中 得到应用；以及【7】的数学背景。17.3完整的系统上诉讨论的方法已经作为一个整体放入粒子成像技术应用，可以分解为发射和搜 集状态发射阶段，多重散射在体元内计算，记录碰撞点，光子能量和入射方向。 光子能量在内部碰撞积累，存在三维数组中，叫做解释缓冲器。搜集状态用标准 alpha方式直观化。模拟系统的输入为各向同性源的位置，虚拟相机的位置以及三维组织密度 场。发射状态迭代的更新三个数组：（图17.6）1. 一维种子缓冲器，储存了并行随机数生成器的种子。这个缓冲器在CPU初始 化，并且在新的随机数生成时由各自的线程更新。2. 一维粒子缓冲器代表当前模拟光子的位置，

23、方向以及能量。3. 三维解释缓冲器。图 17.6模拟程序发射状态体系结构。一个线程由三个函数构成：初始化，发射和积聚。种子缓冲器， 储存了并行随机数生成器的种子。这个缓冲器在CPU初始化，并且在新的随机数生成时由 各自的线程更新，初始化函数更新粒子缓冲器，执行两条指令，决定是否粒子退出还是存活。 它抽样已经退出的和存活的粒子的方向，并且将初始位置和新生的方向替代死亡的粒子。发 射函数，找到下次碰撞位置并决定结束。一个线程分配有自己的种种缓冲器和粒子缓冲器， 所以这里不需要同步。积聚函数将光子能量相加到由光子位置选出的解释缓冲器的一个胞元 中。粒子缓冲器与种子缓冲器由同样的大小，都等于并发的模拟

24、光子数，即，线 程数目，线程与这些缓冲器间有一对一的比对关系。粒子缓冲器由三个线程函数管理，称为循环方式：1. 初始化函数输出新发射粒子到粒子缓冲器中已经死亡的或未初始化的粒子， 并且抽样现有或存活粒子的散射粒子的方向。在死亡粒子的位置，新的粒子 从源位置生成，有随机的方向与初始能量。如果粒子存活且仍活跃，则它的 位置不改变，但是一个新的散射方向需要抽样，而且能量也要根据开普敦法 则改变。生成一个新的粒子，和抽样散射方向在概念上看上去不同，但是其 过程由相似的控制途径。所有形式的随机方向抽样需要两个随机数。如果散 射均匀，则唯一的不同是新生粒子设置到原始位置。通常情况下，初始方向 以及散射方向

25、的不同的公式会引起warp的小分歧。2. 发射函数模拟粒子存在粒子缓冲器中的一步中，它计算下一次散射事件，它 本身执行自由程抽样并且决定吸收是否发生。如果粒子离开体元没有经过碰 撞或者吸收检查告诉我们粒子被吸收了，那它的状态就被调整到死亡，来得 到下一次初始化操作，来在它的位置生成一个新的发射光子。线程最后输出 改变的位置和光子状态。3. 积聚函数搜集入射光子的能力，并在共享存储器中相加。解释缓冲器是三维 的网格，各项代表到达各自胞元的总的发光。量子化光子实际位置的坐标轴 量作为一种指数在积聚缓冲器中，这里入射能力累加。不仅散射光子而且吸 收光子都考虑进去，因为他们都到达该点。注意到这是我们算

26、法中唯一一处 不同线程要写进同一内存的区域，所以同步问题将发生。正确的实现需要信 号或微粒子增加指令来保证没有光子能力丢失。如果我们移除共有的被排除 的地方那一部分粒子将不被加到缓冲器中。而存活的线程远比数据项数目小， 这种粒子丢失可能性很小。然而，模拟过程要跟着上百万粒子，所以极少数的丢失的影响可以忽略。为了看见结果，我们需执行两种方法。第一种观点分配伪彩色给解释缓冲器 每一项，这些决定幅照度，这种可见的方法在医学系统，这里我们关心身体组织 的辐射剂量。辐照度I任)由入射流计算【5】： EI任)=iiAV其中Et是第i个光子到达这个胞元的能量，AV是胞元的体积。为了不仅看 见总能量同样看见光

27、谱，一种辐照度变量在不同频率范围计算，一个光子加到它 适合的频率辐照度上。第二种方法计算辐射到达虚拟相机，并且可以使用，举例，在全局解释器着 色。为了计算辐射，辐照度用材料反照度相乘，以及入射散射光子的概率密度反 映了虚拟相机的方向。17.4结果与评价上述讨论方法的实现是基于CUDA且完全在GPU上运行的。我们坚持这个系统 在NVIDIA GeForce 260 GTX图像处理芯片。当前版本假设单一的各项同源点并 且在一个GPU上运行，但是扩充到多个通常源和多GPU芯片开发以及GPU群 也是直截了当的问题。在我们的实现中，辐射源能用交互式替代，更加适应GPU 组织存储器。通过我们编写的光子传输

28、模拟，交互式的等级也可以在商业上可行的应用， 我们广泛使用蒙特卡罗工具箱，比如像Geant4【2】。 GATE(http:/www.opengatecollaboration.org)， MCNP (http:/mcnp-green.lanl.gov/), SimSet (http:/depts.washington.edu/simset)，或者 PeneloPET【3】，但是这个问题 相同的计算可能会花去数个小时，甚至几天。一种最优化多线程CPU(Intel Q8200)程序能够每秒追踪0.4百万，使用 Woodcock追踪方法依据【14】。我们的实现要模拟36百万和166百万个粒子， 分别

29、在 NVIDIA GeForce 260 GTX 和 NVidia GeForce 480 GTX 这几个 GPU 上运 行；这些加速比分别达到90倍和415倍，执行时间在主函数中的配额分别为初 始化30%，发射为60%，最后累加为10%；通过模拟结果，GPU的带宽范围很 大。我们用门户提供的包来验证我们的程序。我们选取同等大小的密度场，来模 拟光子传输，然后比较结果。当然，在蒙特卡罗方法中，我们不能期待两种方法 完全的相同；因而，我们同样计算变化的结果，并且计算这个不同是否小于标准 要求的背离，对于光粒子，还有另一种方法来确认正确性。如果媒介是可用肉眼 看见的厚，那么扩散方程就是辐射传输很好

30、的逼近，可用解析的方法来解决均匀 介质的问题。单一各项同性源【9】。我们的应用通过了各项测试。2.5 million photons70 msec5 million photons140 msec25 million photons 70Ci insec图 17.7可见的辐射剂量，由一个点源的低能量光子在各项同性媒介中引起。体元大小定义为128 3的 分辨率；宏单元网格有16 3个胞元。2.5 rnKJIon photons5 million photons25 million photons300 msecE00 msec3 sec图 17.8可见的辐射剂量，由一个点源的低能量光子在各项同性

31、媒介中引起。体元大小定义为5123的 分辨率；图17.7和17.8展示了当光子传输时，我们程序的透视结果。图像展示了辐 射程度的编码，通过假定的层次性在图像的场密度上表现。通常，蒙特卡罗方法 的错误与模拟的路径数的平方根成反比，而透视图的完成时间与路径数成正比。 在我们的模拟中，噪声初始化由高架结果显示，甚至能显示视觉上更真实的图像 （相对误差小于1%）在交互时获得速度。图17.8说明方法能扩展到高分辨率模型中区。这里，密度采取人类可见的 数据调整到5123分辨率上。大小尺寸足够大，能够适合我们的GPU。512 initial photon energy 256 kV initial phot

32、on sn&rgy图 17.9由gamma射线引起的可见辐射剂量。体元大小分辨率定义为128 3。图17.9描述了 gamma射线的传输结果，假定两种不同的源出射粒子能量分 别在512KeV和256KeV的级别。注意到，gamma射线的模拟时间与光粒子时间 相同，尽管稍微有些更加复杂。17.5未来的方向当前提出的实现还有几个限制，他们有助于解决系统，变得简单而且小。比 如，瑞利散射和Klein-Nishina模型是由电路模拟转化为代码的，同样的单一各 项同性源也是这样的，他们发射512KeV的光子或者一个光子分成三份频率的， 即红绿蓝。尽管如此，这些设定也是普通的，512KeV级别的光子对比于

33、电子位 子的淹没现象，当光子分为三份时到标准区域在时，在电脑的图像处理中，更普 遍的用户交互都能够扩展应用的领域。限制中将密度场对应到GPU的内存中更加具有挑战意义。因为我们改进的 方法通过利用高分辨率数组来决定是否一个真实的碰撞会发生，否则使用低分辨 率的宏单元数组，这是自然的扩充，来存储低分辨率仅在GPU上的内存，然后 从CPU中复制实际的消光系数当需要的时候。这个系数将允许特别的武断的高 分辨率模型来进行模拟；这些模拟非常需要在医学诊断中。然而，数据通道高分 辨率的是与实际图解不相关的，对于其在GPU上的模拟也并不简单，我们考虑 它是将来更重要的工作。另一种可能性，是考虑程序上的生成模拟

34、，代表高分辨率模型，作为一种算 法来代替高分辨率的数组。对于这些模拟模型，任意有效的方法都将实行。补充章节：并行算法的讨论不同并行算法所运行的层次常用粒度来描述。所谓粒度是指各个多处理器可 独立并行执行的任务大小的量度。大粒度反映可并行性的运算量与程序量大，也 称为粗粒度。相反的，如果并行的运算量小，则为小粒度，也称作细粒度。从不同角度对并行算法分类可以有不同的分法:数值计算和非数值计算的、 同步和异步、SIMD机器上的、MIMD机器上的和VLVI结构上的，等等。数值计算是指基于代数关系运算的一类，诸如矩阵运算、多项式求值、解线 性方程组等计算问题。基本上属于数值分析范畴。非数值计算是指基于比

35、较关系运算的一类诸如排序、选择、搜索、匹配以及 图论等方面的计算问题。基本上是属于符号处理的范畴。数值计算与非数值计算 的分类是针对研究内容所属的范畴分的。同步算法是指某些进程必须等待别的进程的一类并行算法。因为一个进程的 执行依赖于输入数据和系统中断，所以所有进程均需同步在一个时钟，等待最慢 的进程。有人将运行在SIMD机器模型上的算法称为同步并行算法。异步算法是指诸如进程的执行一般不必相互等待的一类并行算法。在此情况 下，进程的通信是通过动态地读取（修改）共享存储器的全局变量来实现。有人 将运行在MIMD共享存储器模型上的算法叫做异步并行算法。同步算法和异步 算法的分类是针对进程间的同步方

36、式分的。就设计方法而言，并行算法的设计方法大体上有三种，他们各有利弊：（1）检测和开发现有的串行算法中固有并行性并将其并行化。该方法代价 小，易于开发。由于算法的思想移植于现有的串行算法，其算法的非并行部分的 正确性和可用性较有保证。但如果现有的串行算法具有内在的顺序性，该算法很 难成功。（2）修改已有的并行算法使其可求解另一个相似问题。当已有的算法与需 解决问题有相似性时，该方法将有很好的效果。但该方法有赖于特定的一类问题， 适用范围狭窄，难以普遍应用。（3）从问题本身的特点出发，从头开始设计一个全新的并行算法。该方法 难道最大，可能遇到的问题最复杂，但可以应用于未知领域的并行算法开发。下面

37、我们讨论一下并行算法性能的评价标准：并行算法的性能评价标准与串行算法有所不同。这是因为影响并行计算执行 时间的因素比串行算法要多。对于串行算法，如果运行的平台确定，保证了所需 的存储空间，基本上可以认为算法的运行时间与算法面临的工作量成正比。由此， 在串行计算中，当算法承担的工作量一定时，算法运行时间可以较为客观的评定 算法性能，但是，在并行计算中，各处理器间的通信开销和负载平衡的开销可能 会很大，它们对算法的运行时间有很大影响。在这种情况下，除了运行时间外， 还需要采用更多的标准来衡量并行算法的性能。（一）运行时间编写并行程序首先我们关心的是并行的实现到底有多快。在消息传递系统 中，在并行算

38、法的整个执行时间中必须考虑通信所耗费的时间。并行执行时间，p 是由两部分组成：一个计算部分tcom和 一个通信部分、mm，即有p comp comm其中计算时间可用类似顺序算法的方法加以推演。通信时间与通信量的大 小、互联网络的速度以及传送方式等有关。例如对PC机群来讲，通信时间就依 赖于许多因素，包括网络速度、网络结构和网络竞争，因而无法用一个简单的公 式来表示，只能给出下面的近似公式：comm starup data其中Ltarup为启动时间，实际上是发送空数据包的消息所需的时间。ata表 示发送单位数据字所需的传送时间，则表示数据字的数目。如果并行程序是运行在由互联网连接各处理器而组成的

39、并行计算机上的，我 们希望减少通信开销，因为一般来说互联网的速度比计算的速度慢的多，网络上 的传输会成为并行程序减少运行时间的瓶颈。由于各节点处理器间的通信将占用 大量时间，在保持有足够并行性的同时，尽量增大计算/通信比是非常重要的。（二）并行加速比并行加速比和并行效率是最早出现也是最常用的并行算法评价标准，它们体 现了在并行计算机上运用并行算法求解实际问题所能获得的好处与付出的代价 之间的关系。并行加速比是定量地描述在对一个串行程序实现并行化的过程中，由于运行 时间的减少而获得的性能。加速比是对并行算法的并行化程度的一种通用的度量 方式。它有两种形式：设T表示用串行机求解某个计算问题所需的时

40、间，二是用p个处理器求解 该问题所需时间，T是用单个处理器求解同一个问题所需的时间，则：1该问题所需时间，T是用单个处理器求解同一问题所需的时间，则：1其一：p个处理器的并行加速比：TS =p Tp其二：并行加速比为八 TS seqp Tp在一般情况下，人S P， 这称为超线性加速，出现超线性加速的情况通常暗示原先的串行算法不是最优化 的，或者是由于使用了某些有利于并行构造的独特的，例如在多处理器的系统中 一般每个处理器都有一定的高速缓存，这些高速缓存的存在减少了机器在I/O上 花费的时间，从而提高机器执行并行程序的速度。如果这种执行速度的提高抵消 了并行计算中通信上多消耗的时间，超线性加速就

41、很有可能会出现。另外，超线 性加速还经常出现在搜索算法中，搜索算法中如果大量存在某些搜索路径先于其 他搜索路径先找到符合搜索条件的结果的情况，则很容易导致整个并行算法的非 线性加速。（三）并行效率由于并行算法的加速比只反映了并行计算的收益而没有反映采取并行计算 所付出的代价，所以需要一个反映并行计算代价的指标。并行效率就是这样一 个指标，它反映了新增处理器对提高加速比的贡献。根据这个指标，用户可以 得知自己在并行计算上的付出与收获之间的关系。并行算法的效率可定义为算法的加速比通处理器数目之比，即SE p p p或c SE p.p p其中，p为处理器个数，当加速比S,接近p时，效率Ep接近于1。

42、若一个 并行算法的S广P，Ep =1，则称这个算法是最优的并行算法。这是理想情况， 但实际上达到S P是困难的，引起效率减小的原因在于：算法中并行度不高， p或缺乏负载平衡；通信、访存冲突以及同步的时间开销随参与计算的处理器数目 的增加而增加。（四）Amdahl定律设某个问题中只能串行执行的运算量百分比为S，其余为p个处理器可以并 行执行的运算量百分比为1 - 5。一般来说，任务量和完成时间在一定条件下是成 正比的，所用的时间比也可以用任务量的比来计算，此时设一个常量k（ k与处 理器的执行速度有关），则T 1x k11 sT = （+ s） x kp p如果我们忽略通信与同步等由并行引起的开

43、销，加速比5,此时应该等于T 1S = = PTp三 + sP这个公式就是Amdahl定律。所以当处理器个数p增加时Sp S加速比上限与处理器无关。在理想的情况下，程序的每由于上式中0 s 1Sp的上界为当问题规模固定时，一个部分都能完全并行，那时p个处理器的加速比似乎应该等于p。实际下，这是不可能的。一个并行程序的加速比能接近于p，这已经是很好的，是高度并行 的。然而加速比越接近于p，其并行性越好。从另一个角度来说，当加速比5,接 近p时，并行效率Ep接近于1。（五）可扩展性可扩展性是评价并行算法性能的重要指标之一。可扩展性的含义是：在确定 的应用背景中，算法（或程序）的性能能否随处理器数的

44、增加而按比例的提高， 它也是网络并行算法在多大程度上能够有效的利用多处理器数增加的能力的一 个量度。随着处理机的增加，如果效率曲线基本保持不变，或者略有下降，则称 该算法在所用的并行机上可扩展性好，如果曲线下降很快，则称可扩展性差。影 响一个并行算法的可扩展性因素很多，如计算方法、粒度、加速比、效率、并行 系统与通信开销等，评判的准则也不尽相同。本文就不再讨论。需要注意的是， 对于一个大规模的并行计算而言，并行算法的效率对于其可扩展性是至关重要 的。参考文献I K. Assie, B. Breton, et al., Monte Carlo simulation in PET and SPECT instrumentation using GATE, Nucl.Instrum. Methods Phy. Res. 527 (1