沈 瑜，孫廣中

1.中國科學(xué)技術(shù)大學(xué) 超級(jí)計(jì)算中心，合肥 230026

2.中國科學(xué)技術(shù)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，合肥 230026

1 前言

在很多需要數(shù)值計(jì)算的領(lǐng)域，如何高效求解大型矩陣的本征值和本征向量是常見的問題。例如在電子結(jié)構(gòu)理論中，人們常常要把求解的Kohn-Sham方程轉(zhuǎn)換為矩陣形式，然后求出本征值[1]。很多時(shí)候系統(tǒng)需要通過多次自洽迭代達(dá)到穩(wěn)態(tài)，而在每次迭代中都需要求解該方程，尤其在幾百個(gè)原子以上的較大尺度下，該方程的求解成為最耗費(fèi)時(shí)間的步驟[2]，因此如何求解Kohn-Sham方程成為影響計(jì)算速度的關(guān)鍵因素之一。

一般情況下，Kohn-Sham方程是一個(gè)廣義本征值方程，具有如下形式：

其中，A和B為已知的N×N的矩陣；c是N×N的本征向量矩陣；λ是以本征值為對(duì)角元的矩陣。很多情況下，A和B都是對(duì)稱矩陣（如果為實(shí)數(shù)矩陣）或者厄米矩陣（如果為復(fù)數(shù)矩陣），同時(shí)B還是正定且非奇異矩陣。這種情況下，可以把廣義本征值方程轉(zhuǎn)化為一般本征值方程，然后進(jìn)行求解。

本征值的求解是數(shù)值計(jì)算的一個(gè)基本問題，直接對(duì)所有本征值和本征向量進(jìn)行代數(shù)求解的計(jì)算復(fù)雜度為O(N3)，在需求解矩陣規(guī)模增長時(shí)，所需要的計(jì)算能力增長得非?？?。但如果問題有一定特殊性，例如只需要一小部分本征值和本征向量，或者本征函數(shù)有一些約束條件等，可以采用迭代等方法來降低計(jì)算復(fù)雜度[3]。在某些情況下，甚至可以繞過求解本征值的步驟，從而使得問題復(fù)雜度降低到O(N)[4]。另一方面，也可以采用大規(guī)模并行的方法來加速本征值的計(jì)算。尤其是近些年超級(jí)計(jì)算機(jī)的飛速發(fā)展，不僅頂級(jí)超級(jí)計(jì)算機(jī)的規(guī)模迅速擴(kuò)大，而且千核乃至萬核級(jí)服務(wù)器也日趨普及（http://www.top500.org）。目前，有很多工作是研究大規(guī)模并行本征值求解，例如ScaLAPACK（scalable linear algebra package）[5]、PLAPACK（parallel linear algebra package）[6-7]、eigen_s和eigen_sx[8]、Elemental[9]、EleMRRR[10]等。同時(shí)，基于眾核的 PLASMA（parallel linear algebra for scalable multicore architectures）和基于GPU的MAGMA（matrix algebra for GPU and multicore architectures）也隨著近年來異構(gòu)計(jì)算的發(fā)展而開發(fā)出來[11-14]。

2011年德國馬克斯普朗克學(xué)會(huì)的科學(xué)家們?yōu)榱私鉀Q在FHI-aims全電子結(jié)構(gòu)軟件中的針對(duì)上千個(gè)原子體系中需要求解的大規(guī)模對(duì)稱和厄米矩陣的本征值問題，發(fā)展了面向P級(jí)應(yīng)用的本征值求解庫（eigenvalue solver for petaflop-applications，ELPA）[15]。ELPA具有相對(duì)傳統(tǒng)ScaLAPACK更好的計(jì)算速度和并行效率，除了在FHI-aims[16]里得了應(yīng)用，也為很多專業(yè)軟件采用，例如針對(duì)納米尺度材料模擬軟件OpenMX（open source package for material explorer）[17]、量子化學(xué)和固體物理軟件CP2K[18]等。

本文將介紹一種自行開發(fā)的基于ELPA的廣義本征值求解器GenELPA，該求解器以開源形式公開發(fā)布（https://git.ustclug.org/yshen/GenELPA_C），它實(shí)現(xiàn)了轉(zhuǎn)換廣義本征值方程為一般本征值方程的方法，然后調(diào)用ELPA進(jìn)行計(jì)算。為了方便使用，采用了類似ScaLAPACK的接口和自動(dòng)選擇最佳轉(zhuǎn)化的方法，并提供了轉(zhuǎn)化和求解分步調(diào)用的函數(shù)。本文將該求解器應(yīng)用到由中國科學(xué)技術(shù)大學(xué)量子信息重點(diǎn)實(shí)驗(yàn)室開發(fā)的第一性原理軟件ABACUS[19]（atomic-orbital based ab-initio computation at USTC，http://abacus.ustc.edu.cn）中，獲得了理想的效果。

2 原理

求解廣義本征值問題可以先將廣義本征值方程轉(zhuǎn)化為一般本征值方程，然后求解轉(zhuǎn)化后的一般本征值問題。第一種轉(zhuǎn)化方法是將B矩陣做Cholesky分解為上三角矩陣或者下三角矩陣，這里以上三角矩陣為例：

其中，U為上三角矩陣；UT是它的轉(zhuǎn)置矩陣。

然后令：

那么方程（1）就變?yōu)闃?biāo)準(zhǔn)本征值問題：

方程（5）可以通過ELPA來求解，如果需要本征向量的話，可以通過下式求得：

第二種方法是將B對(duì)角化后，求得B-1/2。即先用ELPA求解關(guān)于B的本征值方程，得到本征值和本征向量e={ei,i=1,2,…}和q={qij,i,j=1,2,…},然后有：

其中，B-1/2ij為矩陣B-1/2的矩陣元；qik和qkj為本征向量矩陣q的矩陣元；ek為第k個(gè)本征值。令：

則方程（1）可轉(zhuǎn)變?yōu)橐话惚菊髦敌问剑?/p>

同樣，上式可以通過ELPA來求解，如果需要本征向量的話，可以通過下式求得：

3 功能與使用

GenELPA主要功能是自動(dòng)進(jìn)行上述過程，用戶只需要輸入矩陣A和B以及進(jìn)程的相關(guān)信息，就可以直接調(diào)用GenELPA進(jìn)行計(jì)算。用戶可自行選擇轉(zhuǎn)化方法以及是否計(jì)算本征向量。

在實(shí)際的電子結(jié)構(gòu)計(jì)算中，有時(shí)候在進(jìn)行迭代計(jì)算時(shí)，B矩陣保持不變而A矩陣會(huì)隨著迭代的進(jìn)行而不斷更新，因此有必要將B矩陣的分解和求解最終結(jié)果分開，以避免重復(fù)計(jì)算。在GenELPA中，對(duì)這兩步計(jì)算可以單獨(dú)調(diào)用函數(shù)進(jìn)行計(jì)算：第一步將B矩陣做Cholesky分解為上三角矩陣U,然后計(jì)算U-1,或者將B矩陣對(duì)角化后計(jì)算B-1/2；第二步再利用分解后得到的U-1或者B-1/2和A矩陣計(jì)算本征值和本征向量（如果需要的話）。

在第一步中，由于ELPA提供的Cholesky分解函數(shù)在某些較小的系統(tǒng)下可能會(huì)不穩(wěn)定（https://doxygen.cp2k.org/d9/d87/classcp_fm_diag.html#a89cb3ce325957a111071b0d2200bbf17），還提供了使用Sca-LAPACK中的分解函數(shù)作為替代選擇。因此加上對(duì)矩陣B進(jìn)行對(duì)角化的方法，總共有3種方法對(duì)矩陣B進(jìn)行分解。在GenELPA中，使用者可以指定使用某種方法，也可以由程序自動(dòng)選擇。

GenELPA中的函數(shù)采用類似ScaLAPACK的命名規(guī)則：第一個(gè)字母如果是p，表明是個(gè)并行的函數(shù)；第二個(gè)字母表示處理的數(shù)據(jù)類型，d為雙精度實(shí)數(shù)，z為雙精度復(fù)數(shù)。ELPA本身也有兩個(gè)求解引擎“ELPA1”和“ELPA2”，因此需要調(diào)用ELPA的函數(shù)也有兩個(gè)不同的版本，分別對(duì)應(yīng)“ELPA1”和“ELPA2”，其函數(shù)名分別用1和2結(jié)尾區(qū)別。

表1是目前GenELPA提供的函數(shù)。GenELPA的輸入?yún)?shù)類似于ScaLAPACK，此外還加入了在整個(gè)程序中只需要一次初始化的一些參數(shù)，以避免重復(fù)計(jì)算。本文以pdSolveGenEigen1為例說明輸入?yún)?shù)。

其中，各個(gè)參數(shù)意義詳見表2。

使用GenELPA時(shí)，有以下需要說明的地方：

（1）ELPA底層為FORTRAN語言編寫，因此二維數(shù)組（a、b、q、work）中的數(shù)據(jù)在內(nèi)存中的存儲(chǔ)需要按FORTRAN語言使用的列主序方式存儲(chǔ)。至于進(jìn)程矩陣，則不限定采用行主序或者列主序。

Table 1 Function list of GenELPA表1GenELPA函數(shù)列表

Table 2 Parameter list of pdSolveGenEigen1表2 pdSolveGenEigen1參數(shù)列表

（2）程序運(yùn)行中，數(shù)組a、b會(huì)被用作中間變量，因此程序結(jié)束后，變量a、b的值會(huì)被改變。如果希望保留輸入數(shù)組，請(qǐng)事先將輸入數(shù)組復(fù)制一份。

（3）除非必要，否則method方法最好選擇0作為默認(rèn)參數(shù)。當(dāng)使用參數(shù)0時(shí)，在程序結(jié)束后method會(huì)返回實(shí)際使用的方法序號(hào)。

（4）程序目前尚在開發(fā)中，會(huì)盡快采用ELPA2進(jìn)行計(jì)算和處理復(fù)數(shù)方程的部分補(bǔ)齊，同時(shí)會(huì)進(jìn)一步對(duì)性能調(diào)試優(yōu)化。

4 性能測(cè)試

本文使用了20 000×20 000大小的矩陣進(jìn)行了實(shí)際測(cè)試，測(cè)量了不同方法所需要的矩陣分解時(shí)間、本征值計(jì)算時(shí)間和總計(jì)算時(shí)間。使用的測(cè)試平臺(tái)是中國科學(xué)技術(shù)大學(xué)超級(jí)計(jì)算中心的TC4600集群[20]，該集群的每個(gè)計(jì)算節(jié)點(diǎn)有24個(gè)CPU核心，因此測(cè)試使用的CPU核心數(shù)為24的整數(shù)倍，即48、96、192、384和768個(gè)。

Fig.1 Time to calculate all eigenvalues圖1 計(jì)算所有本征值的時(shí)間

Fig.2 Time to calculate half of all eigenvalues圖2 計(jì)算一半本征值的時(shí)間

圖1和圖2分別是計(jì)算全部（20 000）和一半（10 000）本征值的時(shí)間統(tǒng)計(jì)，橫軸為使用的CPU核心數(shù)，縱軸為時(shí)間。其中圖（a）是3種方法所用整體時(shí)間的對(duì)比，黑色方塊線為方法1，紅色圓圈線為方法2，綠色三角線為方法3?？梢钥闯鲈谑褂肅PU核心數(shù)較少時(shí)方法1和方法2的速度幾乎一樣，而使用CPU核心數(shù)較多時(shí)，方法2要稍快一些。至于方法3，所需時(shí)間一直要比其他兩種方法長很多。圖（b）～（d）分別是方法1到方法3的分步驟時(shí)間統(tǒng)計(jì)，其中黑色方塊線為右側(cè)矩陣分解時(shí)間，紅色圓圈線為計(jì)算本征值時(shí)間，綠色三角線為總時(shí)間。由于方法3通過對(duì)角化右側(cè)矩陣進(jìn)行分解時(shí)需要計(jì)算右側(cè)矩陣的本征值，方法3的矩陣分解時(shí)間要接近計(jì)算本征值的時(shí)間。而方法1和方法2的矩陣分解使用Cholesky方法，所需時(shí)間要遠(yuǎn)遠(yuǎn)小于計(jì)算本征值的時(shí)間。

測(cè)試結(jié)果表明：首先，使用Cholesky分解比對(duì)角化能夠更快地將矩陣B分解；其次，在使用CPU核心數(shù)較少時(shí)，采用ScaLAPACK和ELPA提供的Cholesky分解方法耗時(shí)幾乎一樣。

需要說明的是，該測(cè)試程序還沒有得到充分優(yōu)化，因此在大規(guī)模并行效率上還不如ELPA本身的效率，后續(xù)工作中將進(jìn)一步對(duì)其優(yōu)化。

將此求解器應(yīng)用到第一性原理軟件ABACUS，并使用一個(gè)由2 048個(gè)硅原子（Si）組成的晶體進(jìn)行了初步測(cè)試（圖3），該體系在每一個(gè)電子步的迭代中需要計(jì)算一次26 000×26 000大小的廣義本征值方程。同樣在TC4600系統(tǒng)中，采用96個(gè)CPU核心進(jìn)行了測(cè)試。作為對(duì)比，ABACUS原始使用的求解器是由中科院網(wǎng)絡(luò)信息中心開發(fā)的HPSEPS（high performance symmetric eigenproblem software）軟件包[21]。測(cè)試結(jié)果表明，采用GenELPA進(jìn)行計(jì)算后，本征值求解部分計(jì)算時(shí)間由4 300 s減少到2 000 s，程序整體運(yùn)行時(shí)間由4 800 s減少到2 500 s，程序性能取得了顯著提高，如圖4所示。

Fig.3 Asample containing 2 048 silicon atoms圖3 包含2 048個(gè)硅原子的測(cè)試算例

Fig.4 Comparison of computing time圖4 計(jì)算時(shí)間對(duì)比

5 總結(jié)

本文基于具有良好的計(jì)算速度和并行效率的本征值求解器ELPA編寫了廣義本征值求解器GenELPA，該求解器為用戶提供了3種將廣義本征值問題轉(zhuǎn)化為一般本征值問題的方法，簡化了使用方式，避免了在求解某些問題時(shí)ELPA中可能存在的問題；實(shí)際測(cè)試表明該求解器還具有較好的計(jì)算速度和并行效率。針對(duì)電子結(jié)構(gòu)計(jì)算等需要大規(guī)模并行求解廣義本征值的問題，能夠有效地簡化難度，提高效率。

致謝本工作得到了中國科學(xué)技術(shù)大學(xué)中科院量子信息重點(diǎn)實(shí)驗(yàn)室的何力新教授、任新國教授和劉曉輝博士在應(yīng)用集成和實(shí)例測(cè)試中的幫助。

[1]Martin R M.Electronic structure:basic theory and practical methods[M].New York:Cambridge University Press,2004.

[2]Li Pengfei,Liu Xiaohui,Chen Mohan,et al.Large-scale ab initio simulations based on systematically improvable atomic basis[J].Computational Materials Science,2016,112(44):503-517.

[3]Golub G H,van der Vorst H A.Eigenvalue computation in the 20th century[J].Journal of Computational and Applied Mathematics,2000,123(1):35-65.

[4]Goedecker S.Linear scaling electronic structure methods[J].Reviews of Modern Physics,1999,71(4):1085-1123.

[5]Blackford L S,Choi J,Cleary A,et al.ScaLAPACK users'guide[M].Philadelphia:SIAM,1997.

[6]Alpatov P,Baker G,Edwards C,et al.PLAPACK:parallel linear algebra package design overview[C]//Proceedings of the 1997 ACM/IEEE Conference on Supercomputing,San Jose,Nov 15-21,1997.New York:ACM,1997:1-16.

[7]Baker G,Gunnels J,Morrow G,et al.PLAPACK:high performance through high-level abstraction[C]//Proceedings of the 1998 International Conference on Parallel Processing,Minneapolis,Aug 10-14,1998.Washington:IEEE Computer Society,1998:414-422.

[8]Imamura T,Yamada S,Machida M.Development of a highperformance eigensolver on a peta-scale next-generation supercomputer system[C]//Proceedings of the 2010 International Conference on Supercomputing in Nuclear Applications and Monte Carlo,Tokyo,Oct 17-21,2010.Tokyo:Atomic Energy Society of Japan,2011:643-650.

[9]Poulson J,Marker B,van de Geijn R A,et al.Elemental:a new framework for distributed memory dense matrix computations[J].ACM Transactions on Mathematical Software,2013,39(2):13.

[10]Petschow M,Peise E,Bientinesi P.High-performance solvers for dense hermitian eigenproblems[J].SIAM Journal on Scientific Computing,2013,35(1):C1-C22.

[11]Haidar A,Solcà R,Gates M,et al.Leading edge hybrid multi-GPU algorithms for generalized eigenproblems in electronic structure calculations[C]//LNCS 7905:Proceedings of the 28th International Supercomputing Conference,Leipzig,Jun 16-20,2013.Berlin,Heidelberg:Springer,2013:67-80.

[12]Agullo E,Demmel J,Dongarra J,et al.Numerical linear algebra on emerging architectures:the PLASMA and MAGMA projects[J].Journal of Physics:Conference Series,2009,180:012037.

[13]Dongarra J J,Tomov S.Matrix algebra for GPU and multicore architectures(MAGMA)for large petascale systems[R].Knoxville:University of Tennessee,2014.

[14]Buttari A,Dongarra J,Kurzak J,et al.The impact of multicore on math software[C]//LNCS 4699:Proceedings of the 8th International Workshop on Applied Parallel Computing:State of the Art in Scientific Computing,Ume?,Jun 18-21,2006.Berlin,Heidelberg:Springer,2007:1-10.

[15]Marek A,Blum V,Johanni R,et al.The ELPA library:scalable parallel eigenvalue solutions for electronic structure theory and computational science[J].Journal of Physics:Condensed Matter,2014,26(21):213201.

[16]Blum V,Gehrke R,Hanke F,et al.Ab initio molecular simulations with numeric atom-centered orbitals[J].Computer Physics Communications,2009,180(11):2175-2196.

[17]Ozaki T.Variationally optimized atomic orbitals for largescale electronic structures[J].Physical Review B,2003,67(15):155108.

[18]Hutter J,Iannuzzi M,Schiffmann F,et al.CP2K:atomistic simulations of condensed matter systems[J].Wiley Interdisciplinary Reviews:Computational Molecular Science,2014,4(1):15-25.

[19]Liu Xiaohui,Chen Mohan,Li Pengfei,et al.Introduction to first-principles simulation package ABACUS based on systematically improvable atomic orbitals[J].Acta Physica Sinica,2015,64(18):187104.

[20]Li Huimin.User manual of Sugon TC4600 100TFLOPS supercomputing system[EB/OL].(2016-06-15).http://scc.ustc.edu.cn/zlsc/tc4600/tc4600-userguide.pdf.

[21]Zhao Yonghua,Chi Xuebin,Wang Wu.HPSEPS software and parallel computing performance on thousand cores[J].e-Sciense Technology&Application,2010,1(4):20-28.

附中文參考文獻(xiàn):

[19]劉曉輝,陳默涵,李鵬飛，等.基于數(shù)值原子軌道基組的第一性原理計(jì)算軟件ABACUS[J].物理學(xué)報(bào),2015,64(18):187104.

[20]李會(huì)民.曙光TC4600百萬億次超級(jí)計(jì)算系統(tǒng)使用指南[EB/OL].(2016-06-15).http://scc.ustc.edu.cn/zlsc/tc4600/tc4600-userguide.pdf.

[21]趙永華,遲學(xué)斌,王武.HPSEPS軟件包及其在千核上的并行計(jì)算性能[J].科研信息化技術(shù)與應(yīng)用,2010,1(4):20-28.