理论计算机科学电子笔记253（2010）89-104www.elsevier.com/locate/entcsSLAMM -数值算法John M. 丹尼斯1，2国家大气研究Boulder，CO 80305，美国伊丽莎白河Jessup杰瑟普3，4科罗拉多大学计算机科学系Boulder，CO 80309，USAWilliam M. Waite5科罗拉多大学电气与计算机工程系Boulder，CO 80309，USA摘要内存效率正在取代浮点运算的数量，成为数值算法性能的决定因素。从一开始就将内存效率集成到算法中变得更加容易通过计算工具可以量化它的记忆传输。稀疏线性代数内存模型（SLAMM）由一个源到源翻译器实现，该翻译器接受MATLAB指定的算法，并添加代码来预测内存传输。我们对许多小内核和解决稀疏线性系统的算法的完整实现进行的测试表明，SLAMM可以准确地预测从内存层次结构加载的数据量，L1高速缓存在三个不同计算平台上误差在20%以内。SLAMM允许我们在迭代算法的设计阶段快速评估特定选择的内存效率，并且它提供了一种自动化的机制来调整迭代实现。它将执行先验记忆分析的时间从长达几天减少到20分钟。关键词：内存分析，稀疏线性代数，源到源翻译，MATLAB1作者的工作得到了国家科学基金会合作基金NSF 01的支持，该基金为国家大气研究中心（NCAR）提供资金，并通过赠款：#OCI-0749206和#OCE-0825754。额外资金通过能源部CCPP计划赠款#DE-PS 02 - 07 ER 07 -06提供。2电子邮件：dennis@ucar.edu3作者的工作得到了美国国家科学基金会的资助，资助号为CCF-0430646和CCF 0830458。4电子邮件：Elizabeth. Colorado.edu5电子邮件：William. Colorado.edu1571-0661 © 2010 Elsevier B.V. 在CC BY-NC-ND许可下开放访问。doi：10.1016/j.entcs.2010.08.03490J.M. Dennis等人/理论计算机科学电子笔记253（2010）891引言传统上，数值算法的设计是为了在给定的浮点运算次数下达到最佳的数值精度[26，31，30，1，25]。然而，这种方法是基于这样的假设，即执行浮点运算的时间占总计算成本的主导地位。这种假设一般不再成立，因为与浮点运算的成本相比，内存访问的成本已经上升微处理器的性能以每年60%的速度提高，而动态随机存取存储器的性能仅以10%的速度提高。此外，算法开发的进步意味着解决许多问题所需的浮点运算的总数正在减少[19]。因此，求解时间不能再被认为是单独执行的浮点运算次数的函数，而必须被描述为浮点运算和内存访问成本的组合[19]。创建新的算法，证明数值和内存效率[35，7，2，5，6，37，14，20，17]是一项艰巨的任务，尚未得到广泛或系统的解决。最常见的方法是在实现完成之前忽略不幸的是，为了提高内存效率而对代码进行追溯性的重新工具化是一项令人生畏的任务.最好从一开始就将内存效率集成到算法中手动分析涉及数据移动的分析表达式的推导，需要计算机体系结构和软件工程的广泛知识它是费力的，容易出错的，并且过于复杂，无法定期执行。幸运的是，我们可以通过使用计算工具来量化内存传输来自动化这个过程。静态分析通过计算代码本身中的变量引用来估计内存使用情况;动态分析在代码运行时执行，因此反映实际的数据访问。一个通用的设计策略是使用不同的内存布局在MATLAB中建立算法的原型，并通过对原型进行测量来确定最终版本的性能。使用适当的知识收集代码创建有用的原型需要专门的知识;如果我们不打包这些知识，这种方法就不会广泛使用。因此，我们开发了稀疏线性代数记忆模型（SLAMM）处理器，这是一种源到源翻译器，它接受描述候选算法的MATLAB代码，并添加代码块来预测数据移动[10]。源到源翻译器在数值算法的开发中有着悠久的历史1970年初，Marian Gabriel报道了一个翻译器，它接受了一个FORTRAN公式，并产生了一个计算输入公式导数的FORTRAN公式[16]。该翻译器是一个更大的最小二乘拟合系统的一部分[15]。加布里埃尔曲线拟合程序…需要函数的偏导数进行拟合，因此用户必须为这些导数提供FORTRAN表达式J.M. Dennis等人/理论计算机科学电子笔记253（2010）8991SLAMMBody的内存分析： BLGMRES总计：存储需求兆字节（SR）：10.21总计：从L2-> L1 Mbytes（WSL）加载：574.64 +-6.16DGEMM兆字节：0.00 +-0.00稀疏操作MB：79.50 +-6.16Fig. 1. 记忆分析提维斯为一个复杂的函数找到它们，充其量是一件麻烦事。在最坏的情况下，这是一个常见的错误来源多年来，已经开发了一些工具来自动化从声明性规范中生成分析和翻译算法[22];使用它们，现在可以更容易地构建像SLAMM这样的处理器。我曾以《易经》为题，一个公共域编译器生成器，具有对诸如名称分析等标准任务的复杂支持，以生成SLAMM处理器。第2节解释了SLAMM处理器如何分析原型并将其转换为MATLAB程序。在第3节中，我们总结了我们的经验，在应用SLAMM算法开发。2 SLAMM处理器SLAMM的目标是预测一个算法移动的数据量，该算法将在编译语言中实现，给出该算法的MATLAB原型。图1显示了SLAMM为名为BLGMRES的算法提供给用户的信息，BLGMRES是25，228阶矩阵的稀疏线性系统求解器• 所有BLGMRES• 工作集加载（WSL）大小是从内存加载的数据量到处理器，对于BLGMRES是574.64 Mbytes。由于SLAMM在预测某些编程结构的数据移动方面的困难，这个数字可能会有高达6.16 Mbytes的误差。• 没有密集矩阵-矩阵运算可以作为对基本线性代数子程序（BLAS）DGEMM的调用来实现[12]。• 稀疏矩阵-向量操作占总负载的79.50 MB和所有可能的错误。为了获得这些信息，需要对原型进行静态和动态分析。SLAMM直接执行静态分析，并且还产生已经用附加MAT-LAB代码块增强的原型的副本。MATLAB解释器执行转换后的代码以提供动态分析。SLAMM处理器接受包含SLAMM指令的MATLAB原型。所有的SLAMM指令都由前缀%SLM、一个命令和一个终止符组成。命令（表1）是一个关键字，后跟一个或多个参数，参数之间用逗号分隔。92J.M. Dennis等人/理论计算机科学电子笔记253（2010）89表1SLAMM命令开始名称标记MATLAB代码主体的开始，进行分析，并给它命名。endname标记MATLAB代码体的结束名名字FuncStart名称将MATLAB函数的开头标记为分析，并给它命名。FuncEndname标记MATLAB函数名的结束。printname请求为name打印内存分析。函数识别，. . . 将一个或多个标识符分类为函数名对于这些数据移动将被忽略。IncFuncident，. . . 将一个或多个标识符分类为函数名包含大量数据移动。在2.1节中，我们解释了原型中变量引用对应的内存访问次数与它发生的区域之间的关系不幸的是，并不是每个变量引用都具有相同的效果，如第2.2.基于原型的静态分析的翻译在2.3节中描述。2.1区域和变量访问在最简单的情况下，原型中每个变量的出现都代表最终算法对该变量的内存位置的访问。但是，访问内存位置的次数取决于执行路径和包含变量出现的区域。作为一个例子，考虑图2所示的原型。条件保证B2或B3中的一个，但不是两个都执行。类似地，出现在由迭代控制的区域中的变量将被多次访问SLAMM指令还可以用于划分要执行内存分析的区域。在图2中，使用指令定义区域B1并将其命名为foo。print命令会导致类似于图1所描述的输出。名称分析是用于发现和记录程序区域内标识符出现之间的所有关系的过程。这个过程是很好理解的，可以用几个简单的概念来描述[23]：范围：程序文本的区域图2中的区域B 0-B3中的每一个构成范围。对于每个范围，SLAMM计算：总存储需求（SR），从内存加载的数据量，工作集加载（WSL）大小，以及存储到内存的数据量，工作集存储（WSS）大小。定义发生：一个标识符的发生，即使它是合法的，J.M. Dennis等人/理论计算机科学电子笔记253（2010）8993B0B2B3B1图二、计数器关联的作用域这是唯一一次出现这种身份证。SLAMM中的所有已识别事件均为定义事件。对于每个定义事件，SLAMM通过将slm前置到标识符并捕获类型，范围和所需的内存存储来创建新的数据结构。绑定：标识符和唯一实体之间的关系。作用域：绑定有效的程序文本区域。在SLAMM中，绑定在包含标识符的定义性出现的整个范围内保持。同一标识符在嵌套范围中的定义性出现将创建一个新绑定。在新绑定的作用域内，旧绑定无效。因此，新绑定的作用域在旧绑定的作用域中构成了一个名称空间：范围、出现次数、绑定和作用域的集合。SLAMM有两个命名空间：变量命名空间和计数器命名空间。每个标识符出现在每个名称空间中都有一个绑定变量名称空间中唯一相关的范围是那些构成文件和函数的范围。在图2中，区域B0是变量名称空间中的唯一范围。变量名空间中绑定到标识符的实体表示由该标识符命名的原型变量所有范围都与计数器名称空间相关。绑定到计数器名称空间中的标识符的实体表示原型%SLM起始块;r=1;%SLM打印foo;%SLM结束块;结束新%SLM end foo;return-3;其他return [r，b];%SLM start foo;b=rand（1，1）如果b.594J.M. Dennis等人/理论计算机科学电子笔记253（2010）89由该标识符命名的变量在该绑定的范围内访问2.2正在更正访问计数在原型中出现标识符并不一定表示必须从内存层次结构中加载相应的变量。幸运的是，可以通过一系列校正来提高基本标识符计数的准确性。这些校正捕获从MATLAB到编译语言的翻译的属性。表达式n =size（A，1）;的一个高效实现，它使用内置的MATLAB函数size将变量n设置为矩阵A的行维度，不需要访问整个矩阵A。通过访问单个整数变量，可以在编译语言中提供必要的的函数调用校正通过递减出现在函数调用参数列表中的那些标识符的计数来解决计数不匹配。 我们注意不要减少参数列表中表达式中出现的标识符的计数。例如，函数调用校正不适用于m =size（r_m1'*r）;中的向量标识符r m1和r。我们将在下一节中描述函数调用的内存分析方法。点积alpha=在这种情况下，标识符r在一个表达式中出现两次。然而，向量r仅从存储器层级加载一次在一个表达式中重复出现r表示在内存分析计算中不能忽略的缓存重用。为了解决缓存重用问题，我们减少了在单个语句中多次出现的任何标识符的标识符计数。为简单起见，此重复更正不解决多个语句之间缓存重用的可能性表达式r_m1 =r;将向量r复制到向量rm1。在MATLAB中，内存复制对于重命名是必要的，因为MATLAB缺少指针构造。在编译语言中，单个变量赋值可以作为内存复制或指针赋值来实现。我们假设，对于高效实现的算法，具有大存储需求的变量（例如，向量）使用指针赋值。存储需求较小的变量（如单个浮点值）使用内存副本。因为小变量的内存复制成本复制校正递减标识符计数以正确地处理指针分配。到目前为止讨论的校正涉及递减计数，以正确计算特定标识符的发生。 在r =A*w中需要不同形式的校正。 这里SLAMM需要关于由A和w表示的实体类型的附加信息。例如，A可以是稀疏矩阵或密集矩阵。类型信息仅对MATLAB解释器可用我们通过将某些操作符转换为函数调用来解决类型信息的缺乏 例如，为了应用特殊算子校正，我们将表达式r=A*w变换为等价的r=mtimes（A，w）。mtimes的一个profilled版本计算数据J.M. Dennis等人/理论计算机科学电子笔记253（2010）8995B0B2B3B1图3.第三章。图2中代码的翻译在运行时基于类型的确定的移动2.3将SLAMM转换为MATLAB图3示出了SLAMM处理器对图2 的图3中的椭圆框表示被引入以执行动态分析的附加代码块。请注意，SLAMM指令是MATLAB注释，可以作为文档复制到输出标题包含所用MATLAB结构的定义和初始化 为每个范围和标识符积累信息SizeOf代码块在对变量进行赋值后立即插入。它通常包含对MATLABwhos函数的单个调用。下面是图3的范围B2中的assignment和SizeOf块：return[r，b];[slm new] = whos这将在为用户变量命名的生成变量（slm new结构的字段字节%SLM start foo;b=rand（1，1）SizeOf（b）如果b.5new =[r，b，3];SizeOf（new）独占内存分析（B2）其他new =-3;SizeOf（new）独占内存分析（B3）端新%SLM end foo;%SLM起始块;报头r=1;SizeOf（r）包容性记忆分析%SLM打印foo;打印内存分析（foo）%SLM结束块;96J.M. Dennis等人/理论计算机科学电子笔记253（2010）89whos（'new'）调用创建的变量指定变量占用的存储量。新的.独占内存分析块的目的是从特定作用域累积每个范围都需要一个独占内存分析代码块。B1的独占内存分析块（如图3所示）包含（以及其他内容）：slm_foo.wsl = slm_foo.wsl + slm_new.bytes + slm_b.bytes;slm_foo.wss = slm_foo.wss + slm_b.bytes;这通过if表达式和new的获取来计算从存储器hierarchy加载的数据的量，将结果累积在结构slm foo的字段wsl中。类似地，存储到存储器层次结构的数据量（在块开始处分配给b）记录在字段wss中。这些赋值中的每一个都有两个组成部分：在B1中执行代码之前加载或存储的总量slmfoo.wss）和在该执行期间加载或存储的量（例如，slmnew.bytes）。前者在Header中初始化为零，后者是由whos确定的常量。请注意，此计算仅考虑B1绑定范围内的访问。new和b都在B2中有定义性的出现，因此B2在B1绑定的范围中为这些标识符构成了一个洞。与B2中的分配相关联的存储器访问信息由排他存储器分析代码块B2中的代码累积。一般来说，定义变量访问的每个术语都乘以第2.2节中描述的出现次数。对于图3中的所有变量访问，该计数都是1，因此SLAMM省略了它。包含内存分析块的目的是从函数或程序中的所有独占内存分析块收集数据它在函数或程序的文本末尾生成，在那里它将被执行一次。图3的包含性内存分析块包含（除其他外）：slm_foo.wsl = slm_foo.wsl + slm_elseLn7.wsl;slm_foo.wsl = slm_foo.wsl + slm_ifLn5.wsl;这里ifLn5和elseLn7分别是B2和B3的SLAMM打印内存分析（foo）由一个调用库函数组成，SlmPrtAnalysis，并以适当的结构作为其参数。SLAMM区分对应于函数的标识符和对应于变量的标识符。函数被进一步分类为对数据移动有重大贡献的函数和没有贡献的函数。我们将那些对数据移动有重要贡献的函数称为profiled函数。SLAMM在手动构建的表中为常用的内置MATLAB函数集合维护正确的分类。例如，MATLAB函数大小决定变量的范围able通常对数据移动没有贡献，但是函数qr计算输入矩阵的QR分解，当然可以。J.M. Dennis等人/理论计算机科学电子笔记253（2010）8997两种类型的代码转换是profiled函数所必需的。第一个涉及改变函数调用，而第二个涉及改变函数本身。SLAMM通过将字符串SLM作为其名称的前缀并添加输出参数，将函数调用转换为已配置的内置MATLAB函数。额外的输出参数是一个MATLAB结构，其中包含SLAMM计算的系统分析。profiled函数对于返回多个输出参数的profiled函数，代码转换只需要添加一个额外的输出参数。对于返回单个输出参数的profiled函数，可能需要额外的代码转换例如，使用profiled函数的输出参数作为操作数的单个表达式需要生成多个子表达式。所有子表达式都由临时变量链接。生成的表达式t=sin（r）+cos（z）的MATLAB，其中r，z，t∈Rn为：[slm[slm_L1C5，slm_sin L1C5] = SLMsin（r）; [slm_L1C14，slm_cos L1C14] =SLMcos（z）;t= slm_L1C5 +slm_L1C14，在这里，单个原始表达式t=sin（r）+cos（z）被分成三个独立的语句。临时变量slm L1C5和slm L1C14分别是sin和cos函数的原始输出自变量，并用于计算预期结果t。结构slmsinL1C5和slm cos L1C14分别包含sin和cos函数的存储器分析的结果。函数调用转换需要生成新版本的profiled函数。SLAMM为所有必要的内置MATLAB函数提供了一系列profiled函数。每一个都包括对原始函数的调用和内存分析结构的分配。对于用户提供的MATLAB函数，SLAMM语言处理器生成上述各种内存分析代码块，并适当地更改名称和输出参数。3 SLAMM的应用自动内存分析提供了在设计阶段快速评估特定设计选择的内存效率的能力，以及提高预先存在的求解器的内存效率的能力我们通过使用SLAMM来减少并行海洋计划（POP）[32]的执行时间来说明这两种应用，POP是洛斯阿拉莫斯国家实验室开发的全球海洋模型。POP被广泛用作共同体气候系统模型[8]的海洋部分，它使用预处理共轭梯度求解器来更新正压部分的表面压力。分布式内存计算机上的消息传递是通过消息传递接口（MPI）[34]标准支持的我们使用POP版本2.0.1 [29]来检查预处理中的数据移动98J.M. Dennis等人/理论计算机科学电子笔记253（2010）89表2POP中预处理共轭梯度求解器单次迭代的数据移动测试网格 WSLP和WSLM的值以千字节为求解器执行WSLP版本Ultra IIPOWER4R14KWSLM误差WSLM误差WSLM误差PCG2+2D4902v1v2516349055.3%百分之零点一50684865百分之三点四-0.7%57284854百分之十六点九-1.0%PCG2+1D32183164-1.7%3335百分之三点七34737.9%减少百分之三十四百分之三十九百分之三十四百分之三十九使用测试和gx1v3网格的共轭梯度求解器测试网格是随源代码提供的粗略网格，以便于将POP移植到其他计算平台。使用gx1v3网格的POP，在赤道的网格点之间有一个度的间隔，比测试网格的分辨率更高，占国家大气研究中心（NCAR）总计算周期的20%[4]。POP使用三维计算网格。水平维度被分解为逻辑上的矩形二维（2D）块[21]。通过在每个处理器上放置一个或多个2D块，计算网格分布在多个处理器2D数据结构的主要优点是它为矩阵-向量乘法提供了常规的跨一访问。共轭梯度求解器内的2D数据结构的缺点是它包括潜在的大量代表陆地的网格点。因此，将多个明确存储的零添加到矩阵。使用压缩稀疏行存储的替代1D数据结构将避免包括着陆点，但将引入间接寻址。有关数据结构变化的其他详细信息见[11]。我们使用SLAMM来评估基于1D和2D数据结构的共轭梯度求解器所需的数据移动我们描述了两种类型的数据移动，预测的数据移动（WSLP）和测量的数据移动（WSLM）。WSLP由SLAMM预测，WSLM使用硬件性能计数器测量，并引用从L2加载到L1缓存的数据。对于一维数据结构，我们使用MATLAB提供的PCG求解程序。我们在MATLAB中编写了一个新的PCG求解器，它使用与POP中相同的2D数据结构。使用第2节中描述的SLAMM指令，我们为每个算法的一次迭代创建了一个区域。表2的第二列提供了具有2D数据结构（PCG2+2D）的算法的WSLP和测试网格的1D数据结构版本（PCG2+1D）。我们测试了一个天真的PCG2+2D（v1）和优化的一个（v2）的实现，稍后将在此科. 两者的SLAMM预测是相同的。SLAMM预测，与现有的2D数据结构相比，1D数据结构的使用将数据移动减少了34%。基于数据移动减少34%的预期，我们在POP中实现了1D数据结构。如果SLAMM预测了数据移动的轻微减少，甚至数据移动的增加，1D数据结构将永远不会实现。像这样的先验分析提供了一个信心，即实现的编程时间J.M. Dennis等人/理论计算机科学电子笔记253（2010）8999表3微处理器计算平台及其缓存配置的描述CPU公司MhzUltra IISUN400POWER4IBM1300R14KSGI500L1数据缓存32KB32KB32KBL2高速缓存4 MB1440 KB8 MBL3高速缓存–32 MB–！==！代码块：求解器v1！==do iblock=1，nblocksP（：，：，iblock）= Z（：，：，iblock）+P（：，：，iblock）*betaQ = operator（P，iblock）WORK0（：，：，iblock）=Q（：，：，iblock）*P（：，：，iblock）enddodelta=global_sum（WORK0，LMASK）！==！代码块：求解器v2！==delta_local=0.d0do iblock=1，nblocksP（：，：，iblock）= Z（：，：，iblock）+P（：，：，iblock）*betaQ = operator（P，iblock）工作0=Q（：，：，iblock）*P（：，：，iblock）delta_local = delta_local + local_sum（WORK0，LMASK（：，：，iblock））enddodelta=gsum（delta_local）图四、实现版本v1和v2的一部分PCG算法的代码块代码修改或新算法不会被浪费。为了评估我们实现的质量并检查基于SLAMM的预测的准确性，我们使用本地开发的性能分析库（Htrace）来检测所有版本的求解器，该库基于PAPI [27]硬件性能计数器API。HTrace通过跟踪在内存层次结构的不同组件中移动的高速缓存行的数量来计算数据移动我们关注三个主要的微处理器计算平台，它们为从存储器层次结构加载到L1缓存的缓存行提供计数器：Sun Ultra II [24]，IBM POWER 4 [3]和MIPS R14K [36]。表3中提供了每个计算平台的缓存配置说明。测量数据移动（WSLM）也在表2中列出，每个计算平台上的2D数据结构实现（PCG2+2D v1和v2）和1D版本（PCG2+1D）。我们报告了十次迭代中的平均数据移动。虽然对于Ultra II和POWER4平台，PCS 2 +2D v1求解器的测量和预测WSL之间的差异极小，但R14 K的测量值5728 KB比预测值4902 KB大17%。三个计算平台之间的数据移动的差异PCG2+2D v1求解器从内存层次加载的数据比R14 K上所需的多17%，这表明可以提高实现的质量虽然SLAMM提供了潜在性能问题的指示，100J.M. Dennis等人/理论计算机科学电子笔记253（2010）89并不能隔离问题的根源。定位和解决问题仍然是软件开发人员的责任。对PCG2+2D v1求解器源代码的检查表明，对点积计算的微小更改减少了数据移动。图4中提供了求解器版本v1和v2的代码块。函数运算符应用局部矩阵乘积，并且数组LMASK是屏蔽出对应于着陆点的点的数组。在do循环的v1版本中，创建了一个临时数组WORK0，它包含两个向量Q和P的逐点乘积。在do循环之外，LMASK和WORK0数组的点积通过全局求和计算。如果在do循环中访问的数据的大小大于L1缓存，则do循环末尾的WORK0数组的一部分不再位于L1缓存中，必须重新加载才能完成计算。在do循环的优化版本v2中，我们采用标量临时增量局部来累积每个块然后，我们 使 用 一 个 函 数localsum 来 应 用 land mask 来 完 成 点 积 。 最 后 ， 我 们 将 函 数globalsum替换为对gsum的调用。子例程gsum在单个处理器上执行时是delta本地到delta的赋值。因为代码块的版本v2在do循环之外不访问WORK0，所以它潜在地减少了数据移动。表2中R14 K上PCG2+2D求解器v1和v2版本的WSLM请注意，Ultra II和POWER4平台上的数据移动也减少了，但程度较小表2还提供了求解器的预测和测量数据移动之间的相对误差。SLAMM预测PCG2+2D v2求解器的数据移动平均误差在0.6%以内，PCG2+1D求解器的数据移动平均误差在4.4%以内。(The质量较差的PCG2+2D v1的偏差更大。）表2表明，数据移动的实际减少百分比与预测减少非常相似。表2中的结果清楚地表明，自动内存分析可以准确地预测数据移动量，从而在编译语言中实现之前提供特定设计选择的内存效率的先验知识请注意，完全可以在没有SLAMM的情况下手动执行相同的分析。虽然2D情况的计算将是简单的，但1D情况将更加复杂。具体地，1D数据结构的数据移动量的准确预测将需要关于着陆点的位置的详细信息。通过比较与MATLAB原型相关的代码行和内存分析代码，可以估计手动内存分析的难度表4显示了MATLAB中1D PCG原型的非注释行的数量。对于1D PCG原型，原型从41条线增加到总数48行，增加了7%的SLM指令。SLAMM生成的输出代码（原型和内存分析代码）为290行。在这种情况下，内存分析所需的代码行数大约是J.M. Dennis等人/理论计算机科学电子笔记253（2010）89101算法本身。这说明手动计算通常是繁重的并且容易出错，并且因此可能不定期地（或根本不）执行接下来，我们研究减少数据移动对执行时间的影响POP的时间步包括斜压分量和正压分量。正压分量是由一个单一的线性解算器的表面压力。我们在每个计算平台的单个处理器上使用测试网格执行POP，总共执行了20个时间步。这种配置平均需要69次PCG2迭代每时间步长的算法。表5给出了使用求解器的三种实现的正压执行时间（以秒为单位）。我们包括使用2D数据结构的求解器的初始实现的执行时间，以准确地反映自动内存分析对执行时间的总体影响。请注意，PCG1+1D求解器的执行时间始终低于基于2D数据结构的求解器。表5中的最后一行包含PCG2+2D v1与PCG2+1D求解器正压执行时间的减少百分比。表5表明，通过自动内存分析评估或识别的代码修改平均减少了46%的执行时间奇怪的是，表2和表5的比较表明，执行时间减少的百分比差异的原因尚不清楚。接下来，我们将使用gx1v3网格来研究1D数据结构对POWER4平台上并行执行时间的影响。 我们专注于POP在64个POWER4处理器上的执行时间。这是POP的典型配置，NCAR每年消耗约240万CPU小时。我们写一个通用的并行分散例程，在MPI下提供并行执行，解算器的1D版本。最近，一种使用单点积[9]（PCG1）的替代预处理共轭梯度算法被添加到POP中。PCG1算法为并行执行提供了性能优势，因为它消除了一个分布式点积计算。我们确认POP在64个IBM POWER4处理器上总共执行200个时间步，平均每个时间步151次迭代请注意，虽然求解器的成本很重要，表4MATLAB原型的源代码行。开发人员生成源极线生成的SLAMM原型和分析代码原始可持续土地管理百分比417290表5在单个处理器上使用测试求解器执行Ultra IIPOWER4R14KPCG2+2D v121.174.578.58PCG2+2D v220.494.017.97PCG2+1D12.742.114.61减少百分之三十九百分之五十四百分之四十六102J.M. Dennis等人/理论计算机科学电子笔记253（2010）89不支配时间步的总成本表中提供了四种求解器实现的时间步进循环的执行时间（以秒为单位）第六章这些结果表明，使用PCG1+1D求解器与使用PCG1+2D表6在64个POWER4处理器上使用gx1v3网格执行200个时间步的总执行时间求解器实现PCG2+2D v1PCG1+2D v1PCG2+1DPCG1+1D总时间（秒）86.281.578.873.9solver将POP的总执行时间减少了9%。POP的总执行时间减少9%是非常重要的，因为该模型已经被广泛研究和优化[21，33];额外的9%转化为NCAR每年减少216，000个四、结论数据移动是当今计算机上因此，在进行昂贵的实施之前，应探讨设计决策对数据移动的影响。数据移动的预测是可能的，通过添加仪器的MATLAB原型的算法。 就像加布里埃尔近40年前描述的偏导数计算问题一样，这个过程充其量是一个麻烦，最坏的情况是一个常见的错误来源SLAMM自动化修改MATLAB原型的过程，将测试成本从几天减少到几分钟。它封装了有关变量访问特性的知识以及有关在原型运行期间积累信息的策略如果随着时间的推移，有了更详细的知识和更好的战略自动化测试过程的大部分困难在于分析原型文本以选择适当的测试代码的常见任务。由于编译器技术的进步，这些分析任务可以通过声明性规范来描述，代码可以自动生成。这样做可以大大降低开发工具的成本比如SLAMM我们已经表明，SLAMM可以用来指导目前使用的大型代码的性能改善它提供了另一个源代码到源代码翻译器接管软件开发的繁琐和容易出错的方面的实用程序的例子引用[1] E. 安德森Z。拜角，加-地Bischof，J. Demmel，J. J. Dongarra，J. DuCroz，A.格林鲍姆，S. 哈马林 LAPACK用户 SIAM，1992年。[2] A. Baker，J. Dennis，and E. R. 杰瑟普向内存高效的线性求解器。在J.M.L. M 帕尔马，J. 栋加拉河谷Hernandez和A.A. Sousa，编辑，VECPARJ.M. Dennis等人/理论计算机科学电子笔记253（2010）89103《计算科学的高性能计算：论文选编和特邀演讲》，计算机科学讲义第2565卷，第315-327页。施普林格，柏林，2003年。[3] S.贝林河作者：Bell，P. Holtho，F. O'Connell和W.威尔POWER4处理器介绍和调优指南。IBM红皮书，2001年11月。 http://www.redbooks.ibm.com。[4] T. 贝奇国家大气研究中心，12月。2005年个人沟通。[5] J. Bilmes，K. Asanovic，C. W. 阿成和杰·德梅尔使用PHiPAC优化矩阵乘法：一种可移植的、高性能的ANSI C编码方法。在ICSPress.[6] L. Carter，J. Ferrante，and S. F.哈梅尔 分层平铺，提高超标量性能。在IPPSIEEE计算机协会。[7] A. T. Chronopoulos和C. W.齿轮. 对称线性方程组的s步迭代法。 杂志计算与应用数学，25：153 -16 8 ， 19 8 9 。[8] W. D. 柯林斯角M. 比茨湾L. 布莱克蒙湾B. 博南角S. Bretherton，J.A. Carton，P.张，S. C.作者声明：J. B.亨德森，J.T. Kiehl，W. G. Large，D. S.麦肯纳湾D.桑特，和R.D. 史密斯社区气候系统模型：CCSM3。气候杂志：CCSM特刊，11（6），2005年。[9] E. F. D'Azevedo，V. L. Eijkhout和C. H.罗曼共轭梯度算法与分布式内存多处理器上减少同步开销。技术报告56，Lapack工作说明，2002年8月。[10] J. M.丹尼斯自动内存分析：改进迭代算法的设计和实现。博士论文，科罗拉多大学，博尔德，2005年7月[11]J. M. Dennis和E.R. 杰瑟普应用自动内存分析改进迭代算法。SIAM科学计算，29（5）：2210-2223，2007。[12] J. Dongarra，J. DuCroz，S.哈马林和我。你好。算法679：一组3级基本线性代数子程序。ACM Trans.Math. Software，16：18[13] Eli ： 一 个 用 于 编 译 器 构 造 的 集 成 工 具 集 。 文 档 、 示 例 ， 可 从 www.example.com 下 载 http://eli-project.sourceforge.net/。[14] 工程和科学子程序库（ESSL）和并行ESSL。http://www.ibm.com/systems/p/software/essl.html，2007年。[15] M.的加百列用于最小二乘拟合的专用语言。技术报告ANL-7495，应用数学部，阿贡国家实验室，1968年9月[16] M.的加百列一种用于最小二乘拟合的专用语言的符号导数提取器。技术报告ANL-7628，应用数学部，阿贡国家实验室，1970年2月[17] K.后藤河范·德·盖恩。3级BLAS的高性能实施。技术报告TR-2006-23，德克萨斯大学奥斯汀分校，计算机科学系，2006年。[18] R. W.格雷，S。 P. Levi，V. P. Heuring，A. M. Sloane和W. M. 等等 Eli：一个完整的，可扩展的编译器构造系统。Communications of the ACM，35（2）：121http：project.sourceforge.net/。[19] W. D. Gropp，D. K. Kaushik，D. E. Keyes，and B. F.史密斯隐式计算流体动力学代码的实际性能界限。以. 埃塞尔等人，编辑，Proceedings of Parallel CFDElsevier，1999年。[20] 英特尔数学内核图书馆http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm，2007年出版。[21] P. W.琼斯，P.H. Worley，Y.吉田，J.B. III White，and J. Levesque.在并行海洋计划（POP）中的实际性能移植。并发计算练习Exper，17：1317-1327，2005.[22] 联合Kastens，A. M. Sloane和W. M.等等从规范生成软件。琼斯和巴特利特，萨德伯里，MA，2007年。[23] 联合Kastens和W. M.等等属性文法中的模块性与可重用性。Acta Informatica，31：601[24] 孙微系统。Ultra2架构：技术白皮书，2005年。http://pennsun.essc.psu.edu/customerweb/WhitePapers/的网站。104J.M. Dennis等人/理论计算机科学电子笔记253（2010）89[25] A. A.米林河H.科恩湾C.柯蒂斯，W。P. Dannevik，A. M. Dimits，M. A. Duchaineau，D. E. Eliason，D. R. Schikore，S. E. Anderson，D. H.波特，公共关系。伍德沃德湖J. Shieh和S. W.白.基于ibm-sp系统的可压缩湍流的高分辨率模拟。在SupercomputingPress.[26] D. P·奥利里块共轭梯度算法及其相关方法。Linear Algebra and its Applications，29：293[27] 性能应用程序编程接口：用户http：//icl.cs.utk.edu/papi，2005年。[28] D. Patterson，T.安德森，N。卡德威尔河Fromm，K.基顿角科济拉基斯河托马斯
