1. 简介
近年来，材料计算与模拟在新材料研发和材料物性分析中发挥越来越重要的作用，特别是基于量子力学的第一原理(ab initio)计算，构成材料计算与模拟的微观基础 [1] - [8]，由于实验技术和条件的限制，也是对计算需求特别大的阶段。描述微观粒子运动需要求解偏微分方程，对计算能力提出了巨大的挑战。随着高性能计算的能力提升，材料计算软件的并行性能也不断提高，出现了越来越多的适应高性能计算的第一原理并行计算软件，如VASP (Vienna Ab-initio Simulation Package) [9]、QE (Quantum Espresso) [10] 和PWmat [11] 等。考虑到构成宏观物质的微观粒子的数量高达1023量级，现有的高性能计算能力面对材料设计的计算需求，还有不小的距离。
第一原理材料计算软件是基于密度泛函理论(Density Functional Theory) [12] [13] 和周期势条件下电子运动的偏微分方程求解器。计算过程中主要涉及的高性能计算包括：1) 周期势条件下的平面波基组和FFT变换所需的网格分布；2) 矩阵对角化求解DFT基态方程本征值的迭代优化算法。以VASP为例，在网格划分方面，采取了“双网格技术(dual grid)” [14]，即在倒空间用粗网格计算基础上，将网格加密，然后实施FFT到正空间(密网格)，计算有关物理量后，变换回到倒空间(密网格)，最后得到粗网格上的物理量，确保计算精度。而矩阵迭代优化算法则采用了Davison方法、共轭梯度(CG)和RMM-DIIS [15] [16] 等多种优化算法实现。
对称特征值问题作为数值计算的重要问题，存在于很多学科领域中，尤其是材料计算第一性原理，该问题的求解中，矩阵三对角化和三对角矩阵特征值的求解是两个重要的步骤，鉴于问题规模问题，本文中使用分治策略和OpenMPI进行矩阵并行求解，并测试了在分布式模式下的并行求解效率。
本文按以下顺序组织：第二部分介绍HPC并行环境，OpenMPI和LSF集群管理软件，第三部分介绍VASP在HPC作业下的调度、跨节点并行测试和结果进行分析，分别从VASP并行、MPI的体系结构和作业管理系统的资源分配讨论；第四部分为矩阵对角化程序并行算法及测试效果；最后为本文结论。
2. 使用须知
作为以理论计算为基础的计算材料科学来说，计算研究是第一要务。通过模拟计算，可以在不实际制备的前提下深入理解材料的微观结构性质，从而节约研发成本；同时还可以模拟特殊/极限环境下的实验结果，比如对人体健康有害，或者处在压、超低温、强磁场等某些极端条件下，实验测量很难实现或者耗费巨大，可以使用模拟计算来代替或指导实验。
随着HPC能力飞速提高，材料科学应用领域大幅度扩展，更大尺度的体系、更长时间的动力学演化、更精确的理论计算描述成为可能。高性能计算在材料科学中的应用也更加广泛。高性能计算给材料科学研究带来极大便利的同时，材料科学的发展也推动了高性能计算的进步。
2.1. OpenMPI的体系结构
适用于机群、MPP等分布式内存结构的并行编程环境，通常由“并行虚拟机”(Parallel Virtual Machine, PVM)或“消息传递接口”(Message Passing Interface, MPI)来实现。利用PVM工具，可以把互连的各种计算机虚拟为一台并行机，从而为编程人员提供了一个便于管理和使用的编程环境，而由PVM的编译库对程序进行转换，将程序的计算任务分解为若干子任务后合理分配到各个节点机进行并行处理。MPI是一种基于消息传递的并行计算规范，消息(Message)一般包括数据、指令或其它各种控制信号等，MPI提供了一套消息传递库，基于消息传递的并行编程实际上就是通过调用MPI的消息传递库函数实现节点机之间的数据交换，并提供并行处理任务之间的同步等。目前，基于PVM和MPI并行编程环境，都可以支持C、C++和FORTRAN等的并行编程。
Open MPI的基于构件的体系结构不仅为第三方研究提供了稳定的平台，也使得独立软件附件能够在运行时组合。Open MPI的运行时环境提供了启动和管理并行应用的基本服务。Open MPI的设计以MPI构件架构(MPI Component Architecture, MCA)为中心，为所有其它层次提供管理服务的基础组件结构。对Open MPI的架构分析表明，OpenMPI并不会开启额外进程，因此跨节点时的新增进程的唯一来源是作业管理系统对计算资源的分配与监控过程。
2.2. 作业管理系统的资源分配与监控
高性能计算是以大型集群和计算队列为硬件基础。作业管理系统包括集群管理和作业管理两个方面，集群管理系统将一组独立的计算机组合成资源池，共同协作以完成任务。集群广泛应用于高性能计算领域，提供低成本，可扩展及高性能的计算能力，作业管理系统是构成集群的重要软件系统，它的主要任务是对集群的资源进行集中的监控和管理，为用户提交的任务分配可用的计算资源，并监控和管理作业的执行及结果的返回；同时，还提供系统容错和错误恢复能力，可以在异常或故障发生时最大程度的减少任务中断的代价。
随着数据中心在规模和复杂性上的快速增加，使得对集群工作负载和应用的管理更加困难，同时也难以确保计算硬件和软件等资源的正常使用。用户希望在任何地方都能灵活使用应用程序，并自动操控数据流，管理者希望能够监督集群的资源和负载，管理软件使用权，确定瓶颈问题，监督面向用户的服务协议是否满足，并计划系统的扩容数量。
北京市计算中心采用IBM Spectrum LSF (Load Sharing Facility)平台系列软件管理集群和作业。对于要求高的分布式任务型高性能计算环境，该软件能够通过综合的基于智能的，策略驱动的调度策略高效地管理负载，方便客户使用所有计算基础设施资源，保障最佳的应用性能。软件将集群内服务器分为两类，分别是管理节点和计算节点。管理节点负责监控集群中所有节点状态，并分配任务到计算节点，计算节点负责运行用户提交的任务，图1是LSF平台在集群中的系统环境下的任务提交执行过程。
作业管理系统展示了任务提交的完整过程：
1) 作业提交
用户通过LSF客户端，或者可以执行bsub命令的服务器上提交一个作业，当提交这份作业时，如果不指定哪个队列，这份作业就会被提交到系统默认的队列中，作业在队列中等待安排，这些作业处于等待状态。
Figure 1. The submit progress when using LSF in Cluster environment
图1. LSF平台在集群中的系统环境下的任务提交执行过程
2) 作业调度
后台的主进程mbatchd将处理队列中的作业，在一个预定的时间间隔里将这些作业按设定的计划，传递给主调度进程mbschd。
主调度进程mbschd评估这份工作时，根据作业的优先权制定调度决策、调度机制和可利用资源。主调度进程选择可以运行作业的最佳计算节点，并将它的决策返回给后台主进程mbatchd。主负载信息管理进程(lim)收集资源信息，主lim与mbatchd主进程交流这些信息，反过来mbatchd主进程使用之前交流信息支持调度决定。
3) 作业分配
mbatchd主进程收到mbschd发过来的决定后，立即分配作业到计算节点。
4) 作业运行
从批处理进程(sbatchd)，从mbatchd主进程接到要求，为作业创建子sbatchd和一个执行环境，通过使用一个远程执行服务器开始这个作业。对于使用MPI执行的作业，LSF会在MPI和应用中间创建一个TaskStarter进程，用于管理和监控每个作业的执行。
5) 返回输出
当一个作业完成时，如果这个作业没有出现任何问题，它处于一个完成状态。如果有错误作业无法完成，这份作业处于退出状态。sbatchd传达作业信息，包括错误提示和给mbatchd的输出信息。
6) 信息反馈
mbatchd通过预定方式给提交者反馈作业输出信息、作业错误、提示信息、作业信息。
通过上述作业提交和管理、监控过程，不难看出，在跨节点计算任务中，作业管理系统将在MPI通信启动时，在计算节点上开启两类进程，一类用于监督作业运行状况(节点数−1)*CPU数，另一类是管理进程，即检查跨节点计算资源的运行和稳定状况。因此是(节点数−1)。如果计算中的资源不跨节点，作业管理系统启动MPI的进程同时充当监督和守护进程，因此作业管理系统没有启动额外的进程。
3. VASP并行效果测试
VASP是维也纳大学开发PAW方法 [17] [18] 的第一原理材料物性模拟的主要软件，VASP的MPI并行效率较高。在北京市计算中心的服务器的计算环境为Intel® Xeon (R) CPU E5-2680 V3 2.50 GHz，内存DDR4 64 G，2133 MHz，数据网络连接为Infiniband 56 Gb/s。每个节点有24核，每核4线程，通过作业管理系统，分别分配的计算资源包括单节点8核、16核、24核和跨节点24核(N = 2)、48核(N = 3)、72核(N = 4)等几种情况，VASP正常结束后，作业管理系统将统计并输出的最大进程数和线程数。由此确定VASP并行过程中节点、进程等资源的使用情况，为进一步优化计算资源作准备。
3.1. VASP的并行实现
本文实验的运行环境中，VASP的并行计算，是由Open MPI支持的。Open MPI是在LAM/MPI，LA-MPI和FT-MPI的基础上的一种全新的基于构件概念的MPI实现，VASP的MPI实现只关注高性能计算的部分方面，即倒空间布点的分配和能带并行快速求解的问题。Open MPI并不是LAM/MPI，LA-MPI和FT-MPI的简单组合，而是一种全新的MPI实现，完全实现了MPI-1.2和MPI-2规约，并且完全支持并发和多线程应用(也就是MPI_THREAD_MULTIPLE)。
VASP运行过程中，并行计算部分主要集中在主要倒空间网格点的并行分配、k点能带本征值求解(矩阵对角化)和DFT总能量计算这几部分。VASP的主程序将计算资源(分配到的CPU)按控制参数NPAR和NCORE分解(要求满足NPAR*NCORE = NCPU)，在实际并行计算过程中，各部分对计算资源的利用方式不尽相同：
1) 倒空间网格点并行分配和FFT变换，在该阶段，NCPU个计算资源是作为整理对待的，即空间划分的网格被均分到NCPU个进程上，形成网格点和NCPU个进程的拓扑映射，在此基础上，根据NCPU个进程的网络关系完成MPI支持的FFT变换，实现最高效的并发式FFT变换。
2) k点能带本征值求解时，全部CPU分解为NPAR个组，每组NCORE个CPU上实现MPI通信和数据共享(通过共享内存)来求解一个能带的本征值，每次可同时求解NPAR个能带本征值。
3) DFT总能量计算，与k点能带本征值求解时相似，将NPAR组计算的结果汇总，得到体系的总能量。只是简单的MPI进程。
3.2. VASP并行测试
VASP并行计算的节点、进程统计结果列于表1，计算结果表明，单节点(N = 1)时，进程数即CPU数目，当出现跨节点(N > 1)计算时，进程数显著增加，并且进程数随着每个节点上的CPU数目增加呈现出有规律的变化，线程数则始终保持CPU数目的4倍递增。为讨论跨节点计算额外新增进程的可能起源。我们从VASP软件的并行实现、MPI的体系结构和作业管理系统的资源分配和监控三个方面讨论计算过程中启动的进程数目。
从图2中可以看出总的计算时间(紫色线，单位小时)随着并行核数的增加和而降低。当全部核被占据时，并行效率明显降低(计算时间增加)。VASP程序的并行执行主要是自洽迭代(SCF)中矩阵对角化，可以看出，每次自洽迭代的对角化过程RM-DIIS的耗时(蓝色线，单位秒)也是类似规律。在单机最大节点24核的情况下，并行核数达到20核时，RM-DIIS达到并行峰值。
Table 1. The number of threads when using different number of nodes and CPUs
表1. VASP并行计算时不同节点产生进程和线程
跨节点并行时(图3)：满负荷(24核)并行时，随着并行节点数的增加，单次自洽迭代的矩阵对角化时间逐渐降低，但降低的梯度逐渐变小(计算效率有所下降)。半负荷(12核)并行时，计算效率比满负荷并行效率显著提高，由此可以推断，LSF作业调度系统下，每个阶段满负荷运行计算效率不高，考虑到单节点并行的情况，建议20核作为跨节点并行的推荐值。
综上所述，VASP计算过程中，无论实际计算中是否需要跨节点，VASP将NCPU个计算资源作为总体对待，虽然有一定的控制参数将计算资源进行划分，但是VASP运行中并未开启额外的进程，总进程保持为NCPU，因此额外的进程可能由MPI跨节点并行和作业管理系统对作业和计算资源监控产生的。
Figure 2. VASP parallel test single-node
图2. 单节点VASP并行测试
Figure 3. VASP parallel test inter-nodes
图3. 跨节点的VASP并行测试
根据上述讨论，我们可以得到以下结论：在高性能集群中的VASP计算，从作业管理的角度考虑，当指定的计算资源是跨节点计算时，除了VASP作业自身需开启的流程，作业管理系统会开启一部分进程，实现对作业和计算资源的管理、监控和维护。这些额外的进程通过主进程mbatchd监控。一般情况下，这些进程占用的系统资源不大。因此，对于VASP这样的计算作业，由于作业管理产生的进程，对计算资源和计算效率的影响可以忽略不计。但是，由于这些进程伴随作业任务的全过程，因此Wall-time的影响，则必须考虑。
4. 矩阵并行对角化及实验效果
VASP软件中的矩阵对角化计算采用两个阶段，分别采用EDDIAG和RMM-DIIS算法，为了降低问题复杂度并体现MPI并行对矩阵对角化的效率提升，我们采用分治法对矩阵进行分块，每个块独立求解后再用Thomas算法进行归约计算。
4.1. 算法描述
对角矩阵是指只有主对角线上含有非零元素的矩阵，即，已知一个n × n矩阵M，如果对于i ≠ j，Mij = 0，则该矩阵为对角矩阵。如果存在一个矩阵A，使A−1MA的结果为对角矩阵，则称矩阵A将矩阵M对角化。对于一个矩阵来说，不一定存在将其对角化的矩阵，但是任意一个n × n矩阵如果存在n个线性不相关的特征向量，则该矩阵可被对角化。
Figure 4. Pseudo code of algorithm
图4. 矩阵并行对角化程序伪代码
算法描速：
1) 矩阵分块，假设矩阵规模为n，并行数量为p，令k = n/p，则矩阵划分为k*k的矩阵块；
2) 在每个进程中对矩阵进行独立求解；
3) 用Thomas算法对2)中的结果进行归约。
求解算法的伪代码见(图4矩阵并行对角化程序伪代码)
4.2. 实验结果
实验以串行的Thomas算法作为参考，求解规模分为1,000,000，5,000,000，10,000,000，50,000,000四种，分别在2、4、10、16、20、32核的并行度下完成计算。计算时间、加速比和并行效率如表2所示。
Table 2. Parameters of matrix diagonalization
表2. 矩阵对角化参数
实验结果表明利用分治算法，在OpenMPI和分布式并行环境下，矩阵对角化计算具有较强的扩展能力，在文中的规模上并行效率基本稳定，运行时间随着核数增加线性减少，直观效果见图5~8。
Figure 5. Processing time under 1 M, 5 M, 10 M size
图5. 1 M，5 M，10 M矩阵规模下执行时间
Figure 6. Processing time under 50 M size
图6. 50 M矩阵规模下的执行时间
Figure 7. Speedups in different core count
图7. 不同核数的并行加速比
Figure 8. System utilized in different size of problem and core count
图8. 不同矩阵规模和并行核数下系统的使用效率
5. 结论
在HPC环境中的材料计算的核心问题是计算分解、通讯和归约。应用在实际运行中需要兼顾计算资源、通讯效率等多方面问题。在指定的计算资源是跨节点计算时，除了应用软件自身需的资源，操作系统、存储系统和管理软件也会占用一定资源。虽然，这些程序占用的系统资源不大，但对于VASP之类的CPU密集型应用，由于辅助程序产生的进程，对作业任务的全过程均会产生影响，在系统满负荷运行时，会导致计算任务的不均衡问题，不如空置少量核作为预留资源。
矩阵对角化在科学计算领域占有重要地位，如何更加合理地解构问题，处理分块求解通讯和同步问题，保证加速比随着并行度提升线性增加是矩阵优化的重要研究方向。
致谢
本项目由国家重点研发计划项目“产学研用协同的高通量材料计算融合服务平台”(2018YFB0703900)，课题3高通量材料计算的工作流设计与交互图形化(2018YFB0703903)支持。