面向E級計算的線性代數解法器研究綜述

何連花，徐 順，金 鐘

中國科學院 計算機網絡信息中心，北京 100190

數值線性代數算法是大多數科學算法的支柱，大量科學與工程問題最終都歸結為數值線性代數即矩陣計算問題。數值線性代數主要研究線性方程組求解、特征值問題、奇異值分解、最小二乘問題等，目前已經存在非常多的優秀算法來求解這些問題[1-5]。

數值線性代數解法器的發展可以追溯到20世紀70 年代，以LINPACK（linear system package）的產生為標志。此后，產生一系列優秀的數值線性代數解法器，為科學計算提供了有力的支撐。實現塊操作優化的LAPACK（linear algebra package）[6]在20 世紀80 年代被正式推出（http://www.netlib.org/lapack/）。1995年用于并行分布式多核系統的ScaLAPACK（scalable linear algebra package）[7]產生了。隨著多核處理器的出現，關注多核計算的PLASMA（parallel linear algebra software for multi-core architectures）[8-9]在2008年被正式推出。同時，隨著GPU的出現，GPU版本的線性代數庫MAGMA（matrix algebra on GPU and multicore architectures）[10]產生。最近幾年，隨著準E 級機及E級機的出現，出現了很多面向E級機架構設計的數值線性代數解法器，例如美國的SLATE（software for linear algebra targeting exascale）[11]、中國的xMath[12]、日本的KMATHLIB（https://www.r-ccs.riken.jp/labs/lpnctrt/projects.html）。如圖1 所示，可以明顯看出數值線性代數解法器一直隨著計算機體系結構的變化而發展。

圖1 稠密矩陣數值線性代數解法器發展Fig.1 Development of dense matrix numerical linear algebra solvers

本文總結了國外高性能數值線性代數解法器已有工作，特別關注主流稠密及稀疏線性代數解法器面向現代體系架構在功能及性能方面的優化。介紹了稀疏線性代數解法器，包括：PETSc（portable，extensible toolkit for scientific computation）[13]、Trilinos[14]、Hypre[15]、Ginkgo[16]、SuperLU[17]/STRUMPACK（structured matrix package）[18]。

1 E級架構特點和挑戰

如今，高性能計算（high-performance computing，HPC）正處于一個關鍵時刻，計算機體系結構多方面的結合，給數值線性代數算法的發展，既帶來了挑戰，又提供了機遇。本章首先介紹E級架構具有的特點，然后介紹數值線性代數算法研究趨勢。

1.1 E級架構特點

當今高性能計算機具有如下特性[19]：

（1）異構化

目前各種異構架構的CPU（如：x86、ARM、Power）不斷涌現，同時各種基于不同體系結構的計算加速部件（如：GPU、張量處理器（tensor processing unit，TPU）、FPGA（field-programmable gate array）、其他應用專用加速卡）也陸續出現。E級計算機節點數量很大，且每個節點都包含一個重量級多核處理器和許多重量級硬件加速器。大規模并行和異構體系結構給并行編程帶來巨大困難。

（2）數據傳輸成本變高

目前，超級計算機的浮點性能增長速度超過了內存訪問性能或節點間通信性能。這導致計算速度和數據傳輸速度之間存在嚴重差異，預計這種差異在未來會越來越大。可以將數據傳輸性能大概分為延遲和吞吐量。延遲的影響在AllReduce 類型通信中最為突出（如內積、范數）。

（3）能耗增加

目前Top500排名前十的超級計算機的功耗已普遍在MW 甚至10 MW 量級，而E 級計算機功耗將達到30 MW，這給大規模應用的運行帶來了極高的能耗成本。盡管很多芯片設計和制造技術能夠根據實時負載等信息降低微處理器功耗，但在數值算法設計和實現層面依然需要考慮降低能耗。

（4）擁有多種精度的浮點運算

在E級大規模科學計算中，雙精度運算不足的情況會比現在更頻繁地發生。但以四倍精度執行所有運算的成本太高。支持多種精度的GPU異構計算在超算環境已普遍應用，這也為混合精度計算提供了必要的硬件基礎。

（5）故障率增加

隨著技術的進步，大規模集成電路（large-scale integration，LSI）的工藝規則將變得更加精細，從而允許更多的組件集成到一臺機器中。盡管這將顯著提高超級計算機計算能力和效率，但也帶來了可靠性和故障率方面的挑戰。

1.2 數值線性代數算法研究趨勢

算法、軟件和硬件都至關重要。隨著E級機器的誕生，數值線性代數算法也需要不斷更新以便滿足各個領域的應用需求。本節將簡單介紹數值線性代數算法領域最近的一些研究方向（見圖2），這些研究旨在應對1.1節所提出的挑戰。

圖2 數值線性代數算法發展趨勢Fig.2 Research trends of numerical linear algebra algorithms

（1）設計具有高度并行性的算法

為了利用超級計算機的計算核心，數值線性代數算法需要具有高度的并行性。對于稠密矩陣，最近十年基于有向無環圖（directed acyclic graph，DAG）[20]的調度技術已經成為標準技術之一。針對高瘦矩陣的QR分解，提出了基于分而治之策略的TSQR算法[21]。對于稀疏矩陣，目前塊Krylov方法被廣泛使用。

（2）減少數據傳輸量

在E級機上，主內存訪問和節點間通信等數據傳輸的成本遠高于計算。因此，有必要盡量減少數據傳輸量。在內存層次結構中，較高級別的存儲速度更快，但容量有限，而較低級別的存儲容量較大，但數據傳輸速度較低。因此，通過提高數據可重用性，并在數據位于緩存等上層存儲中時對數據執行盡可能多的計算，可以最大限度地減少來自較低層存儲的數據傳輸并提高性能。目前使用各種block 方法（例如block GMRES（generalized minimal residual algorithm）方法）可有效減少數據傳輸量。

（3）降低能耗

這些新的機器特征，要求新一代數值線性代數算法應盡可能最大限度地降低能耗。國內外有很多學者努力通過改進數值線性代數算法以達到降低能耗的效果。2014年，Tan等人[22]對改善高性能數值線性代數運算的功率和能耗的技術進行概述總結。

（4）混合精度算法

目前已經有很多工作關注混合精度數值線性代數算法，Abdelfattah等人[23]對混合精度數值線性代數算法進行了全面的總結。研究表明混合精度算法的性能超過了在高精度下運行的傳統算法。對于性能受限的稠密線性代數算法，采用低精度誤差矯正的混合精度迭代更新（iterative refinement，IR）算法仍然是利用低精度計算能力的首選算法。對于稀疏線性代數，算法的內存約束性質使將內存精度與算術精度解耦的思想更具發展前景。此外，有效的預處理器精度調整，可以在不影響迭代求解器收斂的情況下節省運行時間。2022 年，Higham 等人[24]對數值線性代數中混合精度算法進一步概述總結。

（5）算法層級容錯

目前已經出現一些具有算法級容錯能力的數值線性代數算法，例如Fault-Tolerant GMRES 方法[25]。但目前這方面的研究還遠遠不夠，需要新的范式將算法容錯、應用程序容錯以及系統容錯綜合考慮。

2 稠密線性代數解法器

在稠密數值線性代數中，求解線性方程組Ax=b時，一般先將矩陣A分解為具有一定規律的矩陣的乘積，然后通過向前回代和向后回代（有時計算矩陣的逆），從分解的矩陣中輕松求解x。不同的矩陣分解適用于特定的問題類別。圖1 展示了稠密線性代數解法器的發展歷程，數值線性代數軟件庫在不斷進行優化以便適應新的體系架構，從而發揮出更好的計算性能。本章分別介紹了ScaLAPACK、PLASMA、MAGMA及SLATE的功能，且表1匯總了它們的特點。

表1 ScaLAPACK、PLASMA、MAGMA及SLATE特點匯總Table 1 Summary of ScaLAPACK,PLASMA,MAGMA and SLATE features

2.1 ScaLAPACK

ScaLAPACK（http://www.netlib.org/scalapack/）是一個高性能線性代數庫，適用于分布存儲的多指令流多數據流（multiple instruction multiple data，MIMD）并行機，例如Intel Paragon、IBM SP 系列和Cray T3 系列。此外，該軟件的設計使其可以通過網絡環境與工作站集群一起使用，并通過PVM（parallel virtual machine）或MPI（message passing interface）與異構計算環境一起使用。

ScaLAPACK 可以在任何支持PVM 或MPI 的機器上運行。ScaLAPACK 可以求解線性方程組、線性最小二乘問題、特征值問題和奇異值問題，可以處理矩陣分解或估計條件數等，可求解稠密和帶狀矩陣。ScaLAPACK 主要設計思想包括：將數據矩陣以塊狀循環分布的方式分布到各進程上；使用塊分區算法以確保高質量的數據重用，減少內存層次結構中不同級別之間的數據移動頻率。ScaLAPACK基本構建模塊是PBLAS（parallel basic linear algebra subprograms）和BLACS（basic linear algebra communication subprograms）。

2.2 PLASMA

PLASMA（https://bitbucket.org/icl/plasma）是 由美國田納西大學開發的稠密線性代數包，主要目的是解決LAPACK 和ScaLAPACK 庫在多核處理器上的性能缺陷。PLASMA目前提供了一系列用于求解線性方程組、最小二乘問題、特征值問題和奇異值問題的例程。在過去的十年中，PLASMA 已經被廣泛用于使用英特爾CPU處理器、IBM POWER處理器和ARM 處理器的各種系統中。PLASMA 的設計是為了超越LAPACK 和ScaLAPACK，主要是通過重構軟件結構，盡可能在基于多核處理器的現代計算機上實現更大的效率，同時實現新的或改進的算法。然而，目前由于功能有限，PLASMA 并不能完全取代LAPACK。PLASMA 旨在為多核處理器系統提供高性能服務，這一目標是通過結合并行算法、系統調度和軟件工程中最先進的解決方案實現的。特別是，PLASMA是圍繞以下三個概念構建的。

（1）Tile矩陣布局

PLASMA 采用基于tile 的存儲方法。矩陣被細分為相對較小的方塊，稱為tiles，每個方塊占據一個連續的內存區域。Tiles被有效地加載到緩存中。使用tiles 布局可最大程度地減少沖突cache 失效、TLB（translation lookaside buffer）失效和錯誤共享，并最大限度地提高預取（prefetch）的能力。PLASMA 包含并行和緩存高效的例程，用于在傳統的LAPACK布局和tile 布局之間進行轉換。PLASMA 目前存儲兩個版本的矩陣，因此它的內存需求比LAPACK大。

（2）Tile算法

PLASMA 引入了重新設計的新算法，從而最大限度地提高多核系統緩存級別的數據重用率。tiles被加載到高速緩存中，并在被傳回主內存之前被完全處理。對小tile的操作產生了細粒度的并行性，從而提供足夠的工作以保持大量的核心占用。

（3）動態調度

PLASMA 依賴于順序任務的并發運行時調度。運行時調度基于在任何給定的時間點上根據處理數據的可用性將工作分配給內核的想法，因此有時也稱為數據驅動調度。該概念與通過DAG表達計算的想法密切相關，并且在運行時采用DAG 的靈活性。這與fork-and-join 調度（ScaLAPACK 所采用的）直接相反，fork-and-join 調度中人工同步點暴露了代碼的串行部分，在順序執行時多個內核處于空閑狀態。PLASMA最初使用的是自制的QUARK運行時調度器。目前，PLASMA 依靠OpenMP（open multi-processing）進行動態、基于任務的調度[26]。

DPLASMA（distributed parallel linear algebra software for multicore architectures）[27]通過PaRSEC（https://icl.utk.edu/parsec）運行時系統將PLASMA 算法擴展到分布式環境，并使用CUDA（compute unified device architecture）支持GPU 加速器。在DPLASMA 中，許多算法使用參數化任務圖（parameterized task graph，PTG）來描述任務和依賴關系，其中每個任務通過state-space描述其所有后續任務。在此表示中，不需要在運行時推斷任務依賴關系，并且可以有效地執行PTG。

Chameleon（https://project.inria.fr/chameleon/）[28]以不同方式擴展了PLASMA 算法，從而可以支持各種運行時系統：QUARK、OpenMP、PaRSEC或StarPU（https://starpu.gitlabpages.inria.fr）。其中PaRSEC 和StarPU 支持分布式環境和使用CUDA 的GPU 加速器。在DPLASMA 和Chameleon 中使用高級運行時是一個成功的策略。但面臨的困難是它與運行時相關，這往往需要與運行時的開發者密切協作。例如，增加ROCm或oneAPI的支持很可能需要首先在運行時中增加這種支持。

2.3 MAGMA

MAGMA（https://icl.utk.edu/magma/）是用于異構/混合架構的線性代數軟件包，MAGMA 是由開發LAPACK 和ScaLAPACK 的團隊設計和實現的，融合了混合避免同步和通信避免算法以及動態運行時系統的最新發展。MAGMA 庫的設計在功能、數據存儲和接口方面與LAPACK 相似，從而用戶可以輕松地將他們現有的應用程序從LAPACK 移植到MAGMA，以便充分發揮混合架構的計算能力。

許多科學計算既需要解決大型數值線性代數問題（如BLAS、卷積、奇異值分解線性系統求解器），也需要批量求解可以并行求解的并發小問題。MAGMA包括完善的求解大型矩陣的能力，并支持AI 和高性能數據分析（high performance data analytics，HPDA）、高階有限元方法、圖形分析、神經科學、量子化學、信號處理和許多其他問題領域中發生的較小矩陣計算。例如，在ECP CEED（https://ceed.exascaleproject.org/）項目中，MAGMA 批處理方法用于對許多分組（例如，批處理）的小型高強度GEMM（general matrix to matrix multiplication）進行拆分計算。

MAGMA還支持各種稀疏數據格式，包括壓縮稀疏行格式（compressed sparse row，CSR）、ELL（ELLPACK）、SELL-P 等。目前該團隊的工作重心放在GPU 硬件加速上。基準測試結果所示[29]，MAGMA-sparse在性能方面可以與NVIDIA 稀疏矩陣庫cuSPARSE 8.0 相競爭。

2.4 SLATE

SLATE（https://icl.utk.edu/slate/）[11,30]為當前和未來的分布式內存系統提供基本的稠密線性代數計算能力。SLATE 將提供現有的LAPACK 和ScaLAPACK功能，包括BLAS的并行實現、線性系統求解器、最小二乘法求解器以及奇異值和特征值求解器。在這方面，SLATE 將作為LAPACK 和ScaLAPACK 的替代品，LAPACK 和ScaLAPACK 經過20 年的運行，已經無法完全適應現代的復雜架構體系。

SLATE 的目標是為ECP（exascale computing project）應用提供稠密線性代數計算功能，如EXAALT（exascale atomistics for accuracy，length，and time）、計算化學軟件NWChemEx、量子蒙特卡洛軟件QMCPACK、從頭計算量子化學軟件GAMESS，以及其他軟件庫和框架，如FBSS。SLATE主要特點有：

（1）面向現代硬件體系設計，例如已經面世的E級計算機。

（2）通過依靠標準的計算組件（如BLAS 和LAPACK的供應商實現）和標準的并行編程技術（如MPI、OpenMP）或可移植的運行時系統（如PaRSEC）來保證可移植性。

（3）依靠直觀的單程序多數據（single program/multiple data，SPMD）編程模型和一組簡單的abstractions 來表示稠密矩陣和稠密矩陣運算，從而提高工作效率。

（4）通過使用C++語言的有用工具（例如模板以及函數和運算符的重載）來確保可維護性，并最大程度地減少代碼膨脹。

3 稀疏線性代數解法器

稀疏線性代數對稀疏矩陣進行操作，即在線性方程組Ax=b中，A的矩陣元素主要為零，稀疏矩陣中零元素的比例被稱為稀疏性。因此，它可以利用壓縮數據結構來實現計算和內存效率，其中只存儲非零值及其索引。系數矩陣A的稀疏性決定了兩種類型的線性代數之間的一個基本區別，即稠密線性代數一般是受計算復雜度約束，而稀疏線性代數通常是受內存約束的。

目前有很多優秀的稀疏線性代數解法器。Anzt等人[31]介紹了ECP 項目在為E 級計算平臺提供稀疏求解器方面的挑戰、策略和進展。本章介紹了PETSc、Trilinos、Hypre、Ginkgo 及SuperLU/STRUMPACK 的功能和特點。

3.1 PETSc

PETSc（http://www.mcs.anl.gov/petsc/）是由美國阿貢國家實驗室開發的可移植可擴展科學計算工具包。PETSc是通過MPI實現進程之間的通信，通過調用BLAS 實現數值計算。PETSc 包含線性代數方程解法器、非線性代數方程與無約束優化方程解法器、偏微分方程（partial differential equation，PDE）或常微分方程（ordinary differential equation，ODE）時間依賴方程解法器等，且提供與MATLAB、ESSL、ParMETIS（https://scicomp.ethz.ch/wiki/ParMETIS）等的接口。PETSc 可用于用C、C++、Fortran 和Python 編寫的應用程序。PETSc 支持分布式內存、共享內存和異構CPU-GPU架構，以及任何混合架構。2015年和2016年“戈登·貝爾”獎（Gordon Bell Prize）的獲獎應用均使用了PETSc。

為了更好地利用基于眾核和GPU 的系統，并提供對E級系統的擴展性，PETSc團隊近幾年主要關注以下三點：（1）開發部分無矩陣的可擴展求解器；（2）減少同步算法，與帶有同步點的求解器相比，可以擴展到更大的并發性；（3）所有核心數據結構的性能和數據結構優化。

將PETSc數值庫擴展到E級系統的挑戰有：（1）用于線性、非線性和ODE 求解器的傳統基于稀疏矩陣的技術以及優化算法都受到內存帶寬的限制；（2）任何需要在所有計算單元之間進行的同步（如，內積或范數）都會極大地影響解法器的可擴展性；（3）需要支持E級系統里的各種加速器，以及支持應用程序團隊用于性能可移植性的編程模型。為了解決數值庫的可擴展性問題，PETSc實現了新的求解器和數據結構，包括：延遲使用內積和范數結果的pipeline Krylov方法，允許歸約（reduction）和其他計算的重疊；使用高階方法的部分無矩陣求解器，具有高浮點-內存訪問比率，并且具有良好的使用眾核和GPU 系統的潛力；代數多重網格所需的稀疏矩陣-矩陣乘積的節點內優化，以更好地利用眾核系統。

3.2 Trilinos

Trilinos（https://trilinos.github.io/）是Sandia 國 家實驗室開發的用來解決大規模復雜科學計算問題的科學軟件庫的集合，尤其以線性求解器、非線性求解器、瞬態求解器、優化求解器和不確定性量化（uncertainty quantification，UQ）求解器而聞名，支持并行計算機上復雜的多物理場和多尺度工程應用程序。Trilinos最初是為“flat MPI”而設計的，但近年來節點并發性隨著從多核到眾核架構（例如GPU）的廣泛采用而增長，Trilinos采用了兩級并行模型。在頂層，對象（矩陣、向量）分布在MPI rank 中。分布由映射給出，映射是Petra對象模型中的一階對象。

在Trilinos 中，求解器包括直接稀疏方法，如Amesos2軟件中的稀疏LU/Cholesky，以及Belos中包含的迭代方法（主要是Krylov 方法）[32]。預處理程序包括穩態方法和不完全因式分解（Ifpack2）、區域分解方法和（不）完全因式分解（ShyLU），以及代數多重網格（MueLu）。最近幾年，Trilinos 的主要關注點是在Belos 中開發以加速器為中心的迭代求解器；在MueLu中開發多線程和面向GPU加速的多重網格方法；在ShyLU 中開發基于Kokkos[33]的直接因式分解方法，如Basker[34]和Tacho[35]；在Ifpack2 和ShyLU 中開發面向GPU加速的區域分解方法。

3.3 Hypre

Hypre（http://www.llnl.gov/casc/hypre）是由勞倫斯利弗莫爾國家實驗室開發的軟件包，為在大規模并行計算機上解決大型稀疏線性系統提供了高性能的預處理器和求解器。Hypre 特點之一是提供概念線性系統接口，其中包括結構化、半結構化和傳統的線性代數系統接口。這些接口為用戶提供了一種更自然的方式來描述其線性系統。Hypre 庫提供非結構化和結構化多重網格求解器，在各種高性能計算機上表現出出色的可擴展性，例如Blue Gene 系統（非結構化求解器BoomerAMG 已擴展到125 萬個MPI內核，總共有450萬個硬件線程）[36]。Hypre已經被許多ECP 應用程序團隊所使用，包括ExaAM、Subsurface、ExaWind、CEED等。

非結構化代數多重網格求解器BoomerAMG，它是Hypre 最常用的預處理器，包含啟動和求解兩部分。粗網格算子是通過三個矩陣相乘得到的，這導致粗網格矩陣中每行的非零元數增加，并且相鄰進程的數量也隨之增加。以前在分布式內存機器上實現最佳收斂和性能的插值策略不適合在GPU或類似架構上實現，且計算機體系結構上的通信和計算成本之間的差異越來越大，因此開發試圖最小化通信與計算比率的算法非常重要。這對于具有最優O(N)計算復雜度的多重網格方法來說尤其困難。消除通信通常會降低收斂性，因此可能不會加快整體求解時間（或降低功耗）。但Hypre 團隊已經在許多方法上取得了一些成功，包括：通過稀疏插值技術使三維可變系數擴散問題加速了10 倍以上[37]；具有乘性多重網格收斂速率的高并發加性多重網格方法使三維非結構化Laplace問題的速度提高了2.5倍[38]；用于減少粗網格耦合和通信的非Galerkin 方法使固定系數三維Laplace問題的速度提高了15%至400%[39]。

3.4 Ginkgo

Ginkgo（https://ginkgo-project.github.io/）為迭代求解稀疏線性系統提供了求解器、預處理器以及中心基本構件。它適用于多核和眾核架構，主要是C++實現的，它提供了多種稀疏矩陣格式，例如COO（coordinate，坐標的形式）、CSR、ELL、HYBRID、SELLP。同時，它提供了不同格式之間的高性能轉換程序。

Ginkgo 提供了多種迭代求解器，例如共軛梯度方法（conjugate gradient，CG）、靈活的共軛梯度方法（flexible conjugate gradient，FCG）、雙共軛梯度方法（biconjugate gradient，BiCG）及其穩定版本（biconjugate gradient stabilized method，BiCGSTAB）、廣義最小殘差法（GMRES）等Krylov 子空間方法，還提供更通用的方法，如迭代更新方法（IR）。Ginkgo 還支持帶不完全因式分解預條件的直接和迭代三角求解。基于內存訪問器，Ginkgo團隊部署了壓縮基GMRES（compressed basis GMRES，CB-GMRES）求解器，該求解器通過以較低的精度存儲Krylov基矢從而加速內存訪問，且性能優于GMRES[40]。

Ginkgo 具有一些最先進的通用預處理器，如支持通過動態調整內存精度減少內存帶寬壓力的塊Jacobi預條件[41]，該版本通過動態調整內存精度來減少對內存帶寬的壓力，以滿足數值要求，這種方法已被證明對具有塊結構的問題非常有效[42]。Ginkgo 還具有高度并行的不完全因式分解預條件，如ParILU和ParILUT 預處理[43]、近似逆預條件（如不完全稀疏近似逆預處理）以及多重網格預條件。

由于稀疏算法在幾乎所有的現代架構上都受到內存的限制，人們探索了各種策略來減輕對內存帶寬的壓力，如分層結構的稀疏數據格式壓縮技術和混合精度算法。Ginkgo 解決了內存瓶頸問題，它采用模塊化精度生態系統，將內存精度與算術精度格式解耦，只為算術操作保留全精度，并根據算法要求動態調整內存精度。Ginkgo還包含了最新的混合精度算法[41]，這使得該庫適用于各種科學應用。Ginkgo是超大規模軟件開發工具包（xSDK，https://xsdk.info/）的一部分，是超大規模科學軟件棧（E4S，https://e4s-project.github.io/）的一部分，并且已經集成到一些庫中，如deal.II、MFEM和HyTeG。

3.5 SuperLU/STRUMPACK

稀疏分解和對應的求解算法（如與LU分解相關的三角求解）通常是多物理學和多尺度模擬的線性系統最穩健的算法選擇。它們可以作為直接求解器被使用，例如作為多重網格算法中的粗網格求解器，或者作為迭代求解器的預處理子。SuperLU/STRUMPACK這兩個求解器都是純代數的，適用于各種各樣的應用領域。

SuperLU（https://portal.nersc.gov/project/sparse/superlu/）是一個通用的庫，用于直接求解大型稀疏非對稱的線性方程組，已有超過20年的發展歷史，它適應的體系架構從多核、分布式內存發展到最近的多核GPU 集群。該庫是用C 語言編寫的，可以從C 語言或Fortran程序中調用。它支持實數和復數數據類型，支持單精度和雙精度計算。該庫包含部分主元LU 分解以及通過向前和向后回代進行三角系統求解。LU 分解程序可以處理非方陣，但三角系統求解只針對方陣。SuperLU還包含預估條件數、計算相對向后誤差等。SuperLU 有超過20 年的發展歷史，從多核、分布式內存發展到最近的多核GPU集群，被能源部、工業界和學術界的眾多應用所使用。盡管它們具有廣泛的實用性，但傳統因式分解算法的一個基本擴展障礙是，所需的flops 和內存與問題大小不是線性關系。

STRUMPACK（https://portal.nersc.gov/project/sparse/strumpack/）是一個提供線性代數程序和線性系統求解器的軟件庫，它于2015年首次發布，目前擁有許多能源部的用戶。STRUMPACK用于稀疏和稠密秩結構的線性系統。STRUMPACK 使用multifrontal 算法實現稀疏LU 分解。他們在稀疏因子中的稠密塊（frontal矩陣）中引入了數據稀疏性，使用了幾種分層矩陣低秩壓縮格式，包括分層半分離（hierarchically semi-separable，HSS）、分層非對角線低秩（hierarchically off-diagonal low rank，HODLR）和塊低秩（block low rank，BLR）。除了秩結構壓縮，STRUMPACK 稀疏求解器還支持使用ZFP庫對因子進行壓縮，這是一種針對浮點數據的通用壓縮算法。這可以在指定的精度下使用，也可以使用無損壓縮。

STRUMPACK的稀疏求解器也可以作為精確的直接求解器使用，在這種情況下，它的功能與SuperLU或Superlu_Dist 等類似。STRUMPACK 還提供預條件GMRES和BiCGStab迭代求解器。

4 數值線性代數核心技術進展

新一代數值線性代數算法，除了最大限度地提高速度和準確性外，最大限度地降低能耗也是一個重要考慮標準。這些優化技術包括：提高緩存局部性和向量化程度以最大化利用體系結構特征；深入理解科學應用的動態行為，并在必要時適當降低計算性能和數值精度以滿足降低能效的目標；設計新的通信避免算法以減少節點間數據傳輸。為了利用好現代及未來的超級計算機，Dongarra等人[44]討論了一些高效可行的優化方法。

4.1 性能可移植性

之前的超級計算機均采用NVIDIA GPU，但是最近幾年AMD GPU和Intel GPU也陸續被一些超級計算機所采用，例如：Frontier 和EL Capitan 使用的是AMD GPU，Aurora 使用Intel GPU。表2 列出了主流數值線性代數解法器目前對各GPU 架構的支持情況。針對性能可移植性，目前采取的策略主要有：

表2 對主流GPU架構支持的情況Table 2 Support for mainstream GPU architectures

（1）主要采用隔離影響性能的代碼，編寫CUDA、HIP、SYCL（https://www.khronos.org/sycl）后端；

（2）使用Kokkos和RAJA（https://github.com/LLNL/RAJA）作為主要可移植層；

（3）結合前面兩種策略。

MAGMA 的最新版本是用于CUDA 和HIP 的MAGMA 2.7.1。MAGMA采用多種針對異構平臺的優化方法，例如混合同步、動態任務調度及通信避免。MAGMA 使用混合方法，算法被拆分為不同粒度的任務，并且被分配到相應可用的硬件組件上進行計算。調度可以是靜態的，也可以是動態的，通常在關鍵路徑上安排的小型不可并行化任務都調度到CPU上進行計算，而更大的可并行化的任務（通常是3級BLAS）調度在GPU上進行計算。

SLATE 是建立 在MPI 和OpenMP 等標準以及NVIDIA CUDA和AMD HIP等行業標準解決方案之上的。與ECP Exa MPI 和OMPI-X 項目的合作是為了提高消息傳遞能力，而與ECP SOLLVE 項目的合作旨在改進多線程能力。SLATE 項目還包括為BLAS、LAPACK 以及batch BLAS 設計和實現C++API。在這些API 下，將調用高度優化的供應商庫以實現最大性能，如Intel MKL、IBM ESSL、NVIDIA cuBLAS、AMD rocBLAS。

PETSc性能可移植性主要基于Kokkos和Kokkos Kernels，支持向量、矩陣和矩陣矩陣操作。PETSc/TAO中支持加速器的總體策略基于靈活性和關注點分離（separation of concerns），將用戶編程語言和編程模型的數據指針包裝在PETSc矢量和矩陣矢量對象中。這種方法使PETSc項目組能夠將精力集中在向量、矩陣向量、矩陣矩陣和其他計算的kernels 上，而使用PETSc 的人員可以專注于他們的應用程序。PETSc提供多個后端，并支持AMD、英特爾和NVIDIA加速器。

Trilinos 依靠基于庫的方法來實現可移植性，其中模版化分布式線性代數容器Tpetra 的底層基本都是基于Kokkos 實現。Kokkos 可以輕松地將Tpetra移植到新的計算機架構，并擴展其對并行計算內核和線程可擴展數據結構的使用。Tpetra有以下優點：（1）通過Kokkos 支持多核和眾核架構；（2）模板化的標量類型允許輕松切換精度，包括混合精度；（3）64位整數和指針支持非常大的計算規模。

Hypre 庫之前使用MPI/OpenMP 且不支持GPU。由于當今計算機體系架構以異構為主，Hypre團隊近幾年開始研究使Hypre 能夠利用GPU 的各種方案，如使用CUDA、OpenMP 4.5、RAJA 和Kokkos。他們將后三種方案限制在結構化接口和求解器上，使用了BoxLoops 循環宏。對非結構化接口采用了模塊化的方法，這種方法依賴于重組求解器組件以使用較小的kernels，這些kernels 的CUDA 版本正在實現中。由于HIP與CUDA相似，移植到AMD GPU 上是相當直接的。他們正在研究使用供應商提供的從CUDA到SYCL的轉換工具，將結構化和非結構化的求解器移植到帶有Intel GPU的E級計算機上。

Ginkgo 在線性代數求解器、預條件求解器和矩陣格式方面具有可擴展性。它主要是基于C++（C++11標準）實現，利用各種C++功能（如數據抽象、泛型編程和自動內存管理）來提高庫的性能，同時保持易用性和可維護性。為了最大限度地提高兼容性和可擴展性，Ginkgo 將包含算法的core 與特定架構的后端分開。使用特定架構的“執行器”，可以根據硬件架構和并行化策略的變化而添加、刪除或修改后端。Ginkgo 支持CUDA（用于NVIDIA GPU）、HIP（用于AMD GPU）和OpenMP（用于通用多核處理器，如英特爾、AMD或ARM的處理器）。Ginkgo與硬件供應商的庫（如hipSPARSE 和cuSPARSE）相比具有競爭力。

STRUMPACK使用混合MPI+OpenMP的并行方式。STRUMPACK的求解器和預處理器采用多層次的并行性：稀疏矩陣消去樹的并發性、HSS 層次的并發性和BLAS/（Sca）LAPACK程序的并發性。STRUMPACK 中的多層次并行性使得代碼移植到GPU 架構上非常具有挑戰性。STRUMPACK目前僅支持稀疏直接求解器中的GPU 加速，預條件和秩結構化求解器尚未支持GPU加速。

4.2 混合精度

多年來，科學計算一直在支持IEEE 1985算術標準的硬件上進行單精度（fp32）或雙精度（fp64）運算。2008 年，修訂后的IEEE 標準引入了半精度（fp16）和四倍精度（fp128）格式。半精度僅用于存儲，但現在一些硬件支持將其用于計算。除了IBM z13系統，四倍精度僅在軟件中可用。

期望浮點運算的成本與操作數的位數大致成正比，因此從雙精度到單精度或從單精度到半精度應該會有2 倍加速。較低精度的數據還可以降低存儲要求和數據移動成本。NVIDIA 的Volta 和Turing架構包含tensor cores，可以執行塊FMA（fused multiplyadd）運算。Tensor Cores 支持混合精度計算，動態調整計算以加快吞吐量，同時保持精度。

2019 年，LINPACK 基準和HPL 的提出者Dongarra 對LINPACK 進行擴展，即HPL-AI（https://icl.utk.edu/hpl-ai），以評估超級計算機在混合精度計算上的性能。HPL-AI在矩陣求解過程中放棄雙精度計算的要求，可以利用已有的和即將推出的AI 設備來加速整個求解過程。同時使用IR 算法，恢復在低精度運算中損失的精度。

Abdelfattah等人[45]總結了ECP xSDK-multiprecision開發和部署的混合精度技術進展。

PLASMA為求解一般線性方程組和對稱正定矩陣方程組實現了混合精度例程[9]。PLASMA 采用和LAPACK 一樣的混合精度算法。這些算法基于在迭代更新（IR）過程中以單精度對矩陣進行因式分解，并在迭代更新過程中采用雙精度。

MAGMA 庫支持多個精度，支持基于LU、Cholesky 和QR 分解的混合精度IR 算法，其中采用fp16或fp32作為低精度運算。

SLATE 中實現了兩個混合精度求解器[46]：基于LU 分解的gesvMixed()例程和基于Cholesky 分解的posvMixed()例程。這些例程將低精度和高精度作為模板參數。因式分解在低精度下進行，迭代更新在高精度下進行。

Ginkgo 采用三種機制實現混合精度算法的性能：（1）Ginkgo 具有封裝動態壓縮的內存訪問器，用于將內存精度與算術精度解耦。這會加快內存約束算法的性能，這些算法可以容忍內存操作中的一些信息丟失。（2）Ginkgo 允許不同精度格式存儲的向量不需要進行顯式格式轉換即可進行線性操作。（3）Ginkgo 包含許多有效的混合和多精度算法，如混合精度稀疏近似逆預處理或CB-GMRES方法。

Trilinos、Kokkos Core 和Kokkos Kernels 最近加入了一些混合精度和多精度功能。主要包括：（1）在Kokkos Core 中增加了半精度功能。這使得ECP 應用程序和軟件技術項目可以使用半精度的可移植算術特征實現，并使用半精度類型的Kokkos 編程模型范式。出于性能的考慮，Kokkos也支持對某些操作使用更高的精度。（2）支持半精度稀疏/稠密線性代數。（3）Trilinos在其線性求解器中具有單精度MueLu（多重網格）和IfPack2（基于因式分解的）預處理能力。

PETSc 中的混合精度方法側重于為供應商代數庫（如cuBLAS/cuSPARSE、rocBLAS/roc-SPARSE）和其他支持混合精度的庫（如Ginkgo）開發抽象層。PETSc團隊擬計劃在PETSc 4中增加多精度功能。

Hypre 庫目前可以構建成三種精度（單精度、雙精度和四倍精度）。精度類型的選擇由用戶在構建時通過設置選項--enable-來指定。Hypre團隊正在擴展構建系統，允許統一構建Hypre 支持的所有三個精度，以促進混合精度求解器的開發。

4.3 通信避免算法

通過提高數據可重用性并在數據位于cache 等上層存儲中時對數據執行盡可能多的計算，可以最大限度地減少低級別存儲的數據傳輸并提高性能。

Barron 等人[47]在1960 年首次嘗試減少稠密線性代數中的通信，他們使用EDSAC 2 計算機通過高斯消去法求解線性方程組。他們采用塊LU分解，這種塊分解是在大多數稠密線性代數庫中（如LAPACK、ScaLAPACK、PLASMA）實現的算法的基礎。近年來，為了減少甚至最小化數值線性代數中的通信，人們引入了一種不同的方法，即依靠不同的方式來計算相同的代數運算。例如，矩陣的LU分解是通過使用不同的選主元策略而不是使用列/行選主元來計算的，同時保持數值的穩定性。這些算法被稱為通信避免（communication-avoiding，CA）的算法。ECP PEEKS和Clover項目最近支持的算法重點也是通信避免方法。

SLATE[48]中使用Grigori 等設計的tournament 選主元LU方法（通信避免的LU分解，CALU）[49]。該方法已經被證實是有效的[50]。CALU 的主要思想是通過執行冗余算術運算來減少panel分解期間交換的消息數量。CALU 和經典LU（即列/行選主元LU）之間的主要區別在于panel 分解。在LU 中，處理器需要同步panel的每一列，而在CALU中，處理器只需要同步panel的每個列塊。

雖然Krylov 方法具有良好的收斂性，但它們通常依賴于全局內積（或正交化），這對于大規模的計算體系可能很昂貴。Trilinos 中Belos 軟件包現在支持single-reduce、pipelined 和s-step（CA）求解器[31]，這是減少通信和同步的三種不同方式。single-reduce方法通過將兩個或更多的all-reduce 操作合并為一個，將每次迭代中調用全局（all-）reduce 的數量減少到只有一個。pipelined方法將全局通信、局部通信和計算重疊在一起。s-step方法每s步通信一次，其中s是一個小整數。在CA-GMRES 中，正交化步驟被更有效的塊正交化所取代。代價是需要更多的內存來存儲“ghost”數據。未來的一個潛在改進是采用“矩陣冪kernel”來計算Krylov基。

為了降低算術復雜度，SuperLU團隊開發了一種新型的通信避免的因式分解和三角矩陣求解方法[51]，其中MPI 進程被布置為三維進程網格，沿進程網格的第三個維度選擇性復制數據，用少量增加的內存換取大大減少的每個進程上的通信。由于算術強度較低和任務依賴性較高，求解階段比因式分解更具挑戰性。他們引入了兩種主要技術來減少同步的數量：第一種技術使用基于異步二叉樹的通信方案，通過非阻塞的MPI函數實現，在4 000多個核心上實現了4倍的改進；第二種技術利用one-sided MPI通信功能來實現無同步任務隊列，允許更多的通信和計算重疊，從而在4 000多個核心上獲得額外的2倍改進[52]。

5 總結與展望

為了充分發揮現代超級計算機的計算性能，提升應用計算規模，降低能耗成本，主流稠密及稀疏線性代數解法器面向高性能計算體系架構進行了許多功能及性能方面的優化。本文對主流線性代數解法器進行了總結。針對新一代超級計算架構，解法器的發展主要有兩類表現：一類是新開發的數值線性代數軟件包，如SLATE 和Ginkgo 明確針對新型節點架構從頭設計實現新算法；另一類是已經存在的數值線性代數庫，如Hypre、PETSc和Trilinos等，正在模塊化集成、調整適配式發展。

面向現代計算架構，數值線性代數解法器也需要進行優化，這些優化技術包括：隔離異構計算模塊和設計新的統一編程框架，以實現軟件的性能可移植性；利用混合精度方法提升數值計算和數據存儲的性能水平；發展并實現通信避免算法，避免或減少效率低下的大規模數據通信。目前各解法器在這些方面都取得了階段性進展，但還需要根據E級計算的特點繼續完善。比如，對于硬件驅動的算法設計，了解到該行業將繼續生產針對其他市場（如ML應用程序）的加速器。在不久的將來，GPU加速器（NVIDIA、AMD/ATI、Intel）和具有寬矢量擴展的多核處理器（如ARM SVE和Intel的AVX512）將繼續占據主導地位。因此，還需要加大力度設計及實現適配這些加速器的數值線性代數算法。

面對我國E級超級計算機的特定體系架構，國家也推出一些高性能數值線性代數解法器，從而為頂層應用性能發揮提供了基礎支撐。目前我國高性能數值線性代數解法器要更好地發展，還需從幾個方面努力：首先，加強與應用領域專家合作，實現面向應用需求的解法器性能優化和功能擴充，特別是大規模計算應用驅動的軟件算法發展；其次，加強與國內計算數學領域專家合作，將最新的優秀算法集成到解法器里；最后，加強與具體超級計算機體系架構適配，充分挖掘硬件性能，實現面向特定體系結構的解法器性能優化和功能擴充。