并行廣義本征值求解器GenELPA*

沈 瑜，孫廣中

1.中國科學技術大學 超級計算中心，合肥 230026

2.中國科學技術大學 計算機科學與技術學院，合肥 230026

1 前言

在很多需要數值計算的領域，如何高效求解大型矩陣的本征值和本征向量是常見的問題。例如在電子結構理論中，人們常常要把求解的Kohn-Sham方程轉換為矩陣形式，然后求出本征值[1]。很多時候系統需要通過多次自洽迭代達到穩態，而在每次迭代中都需要求解該方程，尤其在幾百個原子以上的較大尺度下，該方程的求解成為最耗費時間的步驟[2]，因此如何求解Kohn-Sham方程成為影響計算速度的關鍵因素之一。

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

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

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

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

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

2 原理

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

其中，U為上三角矩陣；UT是它的轉置矩陣。

然后令：

那么方程（1）就變為標準本征值問題：

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

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

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

則方程（1）可轉變為一般本征值形式：

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

3 功能與使用

GenELPA主要功能是自動進行上述過程，用戶只需要輸入矩陣A和B以及進程的相關信息，就可以直接調用GenELPA進行計算。用戶可自行選擇轉化方法以及是否計算本征向量。

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

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

GenELPA中的函數采用類似ScaLAPACK的命名規則：第一個字母如果是p，表明是個并行的函數；第二個字母表示處理的數據類型，d為雙精度實數，z為雙精度復數。ELPA本身也有兩個求解引擎“ELPA1”和“ELPA2”，因此需要調用ELPA的函數也有兩個不同的版本，分別對應“ELPA1”和“ELPA2”，其函數名分別用1和2結尾區別。

表1是目前GenELPA提供的函數。GenELPA的輸入參數類似于ScaLAPACK，此外還加入了在整個程序中只需要一次初始化的一些參數，以避免重復計算。本文以pdSolveGenEigen1為例說明輸入參數。

其中，各個參數意義詳見表2。

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

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

Table 1 Function list of GenELPA表1GenELPA函數列表

Table 2 Parameter list of pdSolveGenEigen1表2 pdSolveGenEigen1參數列表

（2）程序運行中，數組a、b會被用作中間變量，因此程序結束后，變量a、b的值會被改變。如果希望保留輸入數組，請事先將輸入數組復制一份。

（3）除非必要，否則method方法最好選擇0作為默認參數。當使用參數0時，在程序結束后method會返回實際使用的方法序號。

（4）程序目前尚在開發中，會盡快采用ELPA2進行計算和處理復數方程的部分補齊，同時會進一步對性能調試優化。

4 性能測試

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

Fig.1 Time to calculate all eigenvalues圖1 計算所有本征值的時間

Fig.2 Time to calculate half of all eigenvalues圖2 計算一半本征值的時間

圖1和圖2分別是計算全部（20 000）和一半（10 000）本征值的時間統計，橫軸為使用的CPU核心數，縱軸為時間。其中圖（a）是3種方法所用整體時間的對比，黑色方塊線為方法1，紅色圓圈線為方法2，綠色三角線為方法3。可以看出在使用CPU核心數較少時方法1和方法2的速度幾乎一樣，而使用CPU核心數較多時，方法2要稍快一些。至于方法3，所需時間一直要比其他兩種方法長很多。圖（b）～（d）分別是方法1到方法3的分步驟時間統計，其中黑色方塊線為右側矩陣分解時間，紅色圓圈線為計算本征值時間，綠色三角線為總時間。由于方法3通過對角化右側矩陣進行分解時需要計算右側矩陣的本征值，方法3的矩陣分解時間要接近計算本征值的時間。而方法1和方法2的矩陣分解使用Cholesky方法，所需時間要遠遠小于計算本征值的時間。

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

需要說明的是，該測試程序還沒有得到充分優化，因此在大規模并行效率上還不如ELPA本身的效率，后續工作中將進一步對其優化。

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

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

Fig.4 Comparison of computing time圖4 計算時間對比

5 總結

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

致謝本工作得到了中國科學技術大學中科院量子信息重點實驗室的何力新教授、任新國教授和劉曉輝博士在應用集成和實例測試中的幫助。

[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.

附中文參考文獻:

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

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

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