分布式严格对角化

• 引言
• 实例与源码速览
  • 问题提出
  • 代码
  • benchmark
  • 提醒
• PETSc 与 SLEPc 简介
  • 程序基本结构
  • Vec
  • Mat
  • EPS
  • 程序编译运行
  • 天河2号安装运行记录
• 本征值问题黑魔法掠影
  • Krylov Methods
  • preconditioning
  • 本征值问题解决框架
  • 稀疏矩阵求解大量本征值的讨论
  • Harmonic extraction
  • 谱变换与本征值求解
• 线性代数各库的关系
• Reference

引言

本文从一篇关于 MBL 严格对角化 (ED) 计算的文章¹及其对应的代码²开始，对分布式内存的矩阵计算与对角化的技术，生态和代码进行一点梳理和记录。本文的重点在于对 PETSc 和 SLEPc 这两个分布式矩阵运算库的使用做一个简要的介绍。对于只是想了解分布式内存矩阵运算技术和相应库基本的用法的读者，可以跳过和物理背景相关的第一部分。

实例与源码速览

问题提出

问题的来源是处理多体局域化（MBL）问题，由于 MBL 没有什么好的解析手段，只能用严格对角化来处理，求得能谱。而这么做内存就极大的限制了能处理的系统大小。而 MBL 的标准模型就是简单的 XXZ model。对于一个 N sites 的 spin \(1/2\) 系统，在 \(S_z=0\) 的 sector，要想求得全部能谱，需要的内存是 \(C_0(C_N^{N/2})^2*8\) ，其中 \(C_0\) 是一个大于 1 的常数。因此对于通常的 32G 内存的工作站，想要计算 18 个 sites 的系统就十分勉强，甚至做不到了。本来，在这种情况下，大家通常求助于稀疏矩阵来减小内存使用，同时利用 Lanczos 等迭代手段，只求我们感兴趣的少数本征值。这种方法虽然适合求基态，但不适合 MBL 问题。因为对于 MBL 我们关注的是所谓的 infinite temperature，也就是能谱中间的情形。而由于这里 DOS 很大，所以很难适用于一般的迭代手段直接求解。

这篇文章则利用了 shift-invert 变换，改为计算 \((H-\lambda_\star)^{-1}\) 的能谱，将 Hamiltonian 指定位置 \(\lambda_\star\) 处的谱变得间距更大，从而实现了稀疏矩阵迭代求能谱。该方法将系统 size 扩大到了 24。看起来这似乎只是 6 个 sites 的提升，但换句话说，这在 Hilbert space 是 \(4^6\sim 4000\) 倍的提升！不过作者这一结果也是利用了2000个节点同时并行来计算的(每个节点24个核)，这一算力也很恐怖了。

代码

作者提供的 demo 代码²，由 Hamiltonian, Basis, Operator 三个子类组成，分别用来定义矩阵元，生成满足对称性的基，和对于波函数进行处理和操作。而 main 部分则利用了 SLEPc 库来进行分布式对角化操作。其中这里的自旋体系的 basis 用 std::bitset 表示非常方便。三个类里与最后利用的 SLEPc 有关的方法包括 Basis::get_mpi_ownership_range() 负责向不同进程分配对应矩阵存储的起止行数。XXZHamiltonian::calculate_nnz() 负责计算每个节点存储部分中每行的非零的对角块和非对角块的元素个数。XXZHamiltonian::calculate_sparse_rows() 负责矩阵分条生成，数据可以供 PETSc 的 MatSetValues 方便使用。 Operator::Siz() 则对波函数（分布式向量）作用算符 \(S_i^z\)，得到新的分布式向量。

至于 sinvert.cc 文件的 main 函数，由于和 SLEPc 的使用高度相关，故下文主要讨论 SLEPc 的使用，理解了 SLEPc 的工作方式，main 的代码就都显然了。

benchmark

这里复述一下作者 benchmark 的一些数据，以供大家参考。经过测试，作者指出对于该问题高效的 solver 是 MUMPS 和 STRUMPACK。该算法在 “reasonable” 的资源下应该可以计算到 22 sites。不过作者计算的本征值数量比较少，只取中间10个。对于 \(L\leq 18\) 的系统，最佳策略是单进程计算。这里的并行策略有两种，因为同时具有 MPI based 的分布式多进程和 OpenMP based 多线程。对于 MUMPS solver， full MPI 策略更好一些，而 STRUMPACK 则反之。对于更大的系统，STRUMPACK 在内存和时间上更优一些。

另一方面，这一对角化所需的时间基本和系统在相图的位置无关，使得该方法更适合大范围的研究 MBL 系统。从计算时间考虑，一般我们不计算超过 100 个本征值。另一方面，作者也测试了 float 精度数据的对角化，这可能造成很小一部分本征向量出现混合。考虑到其只能 reduce 一半内存，可能没什么必要使用。

提醒

原则上，按照原文章 step by step 就可以成功编译运行对应代码。不过还是有一些细节，文章有点问题或没有提。

1. 编译程序时，cmake 之前，需要先添加环境变量，. ./src/conf/linux.cfg，注意要先把该文件中的对应变量改成符合自己机器的设置。
2. 此外需要在 CMakeList.txt 中添加一行 set (CMAKE_CXX_STANDARD 11) 确保 c++11 的相关语法可以编译通过。此外 cmake 要求的版本较高，linux 发行版包管理提供的默认版本可能不足，有可能需要手动编译安装最新的 cmake。
3. 文章给出的 slepc.options 示例之中，最后一项需要改为 st_pc_factor_mat_solver_type。(注：这里有一个 API 变更，老版本v3.8的选项仍为 -st_pc_factor_mat_solver_package)

PETSc 与 SLEPc 简介

PETSc 是 Portable，Extensible Toolkit for Scientific Computation 的缩写。其重点关注的是分布式解决 PDE 问题，不过由于其提供了分布式向量和矩阵的存储与操作的基础设施，因此是上述代码解决分布式对角化问题所必需的。所以为了进行分布式对角化，我们需要关注的主要模块是 Vec 和 Mat 这种基础设施的声明，构造和组装。当然由于求本征值和线性方程组问题的相似之处，以及谱变换方法解决本征值问题时，利用到了解线性方程组，KSP 对应的解线性方程组模块和相关联的 PC 对应的 precondition 模块的原理和选项也最好有些了解。

SLEPc 是 Scalable Library for Eigenvalue Problem Computations 的缩写，这一库构造在 PETSc 提供的基础设施之上，专注于解决分布式内存稀疏矩阵的本征值问题及相关问题。对应的模块为 EPS，其中还隐含使用了谱变换模块 ST。该库提供了和 PETSc 风格高度一致的 API。

下面我们结合两者的用户手册³⁴， 对我们可能使用的模块和 API 进行一下简要梳理，其中涉及算法原理的一些复杂的 API 和选项可能略过或放在下节，结合算法具体讨论。这里只会给出一些 API 的关键词，具体的用法 manual 里都写的比较清楚。

程序基本结构

由于我们讨论的是解决本征值问题，因此程序框架都是调用的 SLEPc 中的 API，其中很多 API 都是 PETSc 中 API 的 wrapper，如果只使用 PETSc 工作的话，也可以将对应的 API 替换为 PETSc 版本。一个程序最基本的骨架如下所示。如果只利用 PETSc，需要包含的头文件为 <petscksp.h>

# include <slepceps.h> // other heads in slepc and petsc are auto included

static char help[] = "this is help content";

int main(int argc, char **argv)
{
    SlepcInitialize(&argc, &argv, "filename", help); // initilize both slepc, petsc, and mpi. PETSC_COMM_WORLD is just alias for MPI_COMM_WORLD later
    PestcOptionsGetInt(NULL, NULL, "-n", &n, NULL); // read options from file and command lines, (using this default option database due to the first NULL parameter)
	// main program here
	SlepcFinalize();
}

稍微提一句这两个库 API 的风格。虽然有很强的面向对象的特征，但由于需要和 C 兼容， API 的调用还是完全命令式的。返回值也基本采取函数变量指针的形式来获取。函数的命名以对象开始，比如 VecCreate() 按照 c++ 风格对应的就是 Vec.Create()。因此整个 API 风格是藏在过程式背后的面向对象。

Vec

PETSc 提供 sequential 和 parallel 两种 Vec，其中后者是基于 MPI 的，也即分布式计算。其 API 分别是 VecCreateSeq(PESTC_COMM_SELF, PestcInt m, Vec *x) 和 VecCreateMPI(PETSC_COMM_WORLD，PestcInt m, PestInt M， Vec *x)，其中的 m 指的是在本进程中存储的元素数，M 是向量的总长。一般来说，指定一个即可，另一个直接用 PETSC_DECIDE 宏来代替，其可直接选用合适的默认值。另一方面，也可以把 create 过程拆开，拆成 VecCreate, VecSetSizes, 和 VecSetFromOptions 三个 API，而向量的类型，则有命令行选项 -vec_type mpi 给定。

对于向量的赋值，则有 VecSet 和 VecSetValues 两个 API，前者将向量所有分量赋为同一个值，后者则可按指标赋值。当赋值完成时，需要调用 VecAssemblyBegin 和 VecAssmeblyEnd 才可以使得向量完成。

对于向量的取值，可以通过 VecGetValues 来取得本地进程存储的向量分量的值。向量本地存储的范围可以通过 VecGetOwnershipRange 来获得，如果只想获得本地向量的大小，可以利用 VecGetLocalSize。如果想拿到向量分量对象本身，需要使用 VecGetArray 和 VecRestoreArray 来完成生成好的向量的后续修改。注意这里返回的指针也是只能访问本地存储的部分向量的。如果想把分布式的向量收集到某个单独的进程上，则需要一系列和 VecScatter 相关的 API。

VecDuplicate 可以完成向量的复制。VecDestroy 可以完成向量的回收。

向量运算的主要 API 包括： VecNorm， VecDot, VecScale, VecSum, VecPointwiseMult, VecMax 等等。其功能看名字也知道大半，具体细节和其他向量计算函数可以查阅 manual。

Mat

矩阵创建的 API 和赋值的API 都和向量类似，只不过把 Vec 改成 Mat 即可。稀疏矩阵默认的存储格式是 AIJ。矩阵通常分条（一些矩阵行）存储在不同进程的内存里。对于稀疏矩阵，一切没有赋值的元素在 Assembly 过程中都从内存中去除，后续再添加非零元素的 cost 很大。因此可以提前用 MatSetValues 插入 0 占位。MatSetOption API 可以接受一个已经 assemble 好的矩阵，从而快速复用其存储和通讯结构以及非0 pattern。另一方面，直接使用 MatCreateAIJ 可以更好的控制矩阵底层的内存分配（直接指定非零元素位置）。矩阵的具体存储结构如下图所示。关于矩阵运算的 API 请参阅文档。

EPS

EPS 是用来解决本征值问题的 context （C style 的面向对象），由 SLEPc 提供。其基本使用流程如下：

// declaration of variables
EPS eps; /* eigensolver context */
Mat A; /* matrix of Ax=kx */
Vec xr, xi; /* eigenvector real part and img part, x */
PetscScalar kr, ki; /* eigenvalue real part and img par, k */
/* difference between PestcScalar and PestcReal: scalar type can be adjust by options --with-scalar-type, --with-precision */
PetscInt j, nconv; /*nconv: number of converged eigenvalues*/
PetscReal error;

//solve the problem
EPSCreate( PETSC_COMM_WORLD, &eps );
EPSSetOperators( eps, A, NULL ); // for problem Ax=kBx, B=I as indicated by NULL here
EPSSetProblemType( eps, EPS_NHEP ); // the second parameter indicates the matrix type, symmetries and Hermicity
EPSSetFromOptions( eps ); // configure options from CLI or file
EPSSolve( eps ); // main part to solve eigenvalue problem
EPSGetConverged( eps, &nconv ); // learn about how many eigenvalues are converged
for (j=0; j<nconv; j++) {
EPSGetEigenpair( eps, j, &kr, &ki, xr, xi ); // get eigenvalues and eigen functions
EPSComputeError( eps, j, EPS_ERROR_RELATIVE, &error );
}
EPSDestroy( &eps ); // clear the context

了解了以上基本流程后，我们就可以把重点放在 EPSSetFromOptions 可以给定的选项里，这些选项构成了解决本征值问题利用方法的丰富性。这些选项还包括了 EPS 隐藏的 ST 对象的设置，这一对象的选项又包含了对 PETSc 中 KSP 和 PC 对象的选项设置。需要注意，这些参数我们都可以不去管，也会有默认值可以正常计算，所以设置这些选项都属于进阶需求。

-eps_nev ：这一选项给出我们希望获得的本征值的数目，-eps_ncv： 这一选项给出我们用来计算的子空间的最大维度，一般要取为 nev 的两倍以上。通过 EPSSetWhichEigenairs ， EPSSetTarget 和 EPSSetInterval API，我们可以对想要求的本征值的范围进行指定，比如最小的本征值，离某个具体值最近的本征值等。这些 API 都有对应的命令行选项，因此可以在 runtime 设置，具体参数参考 manual。值得一提的是如果想要求能谱中间的某个本征值，往往普通的迭代方法由于收敛很慢，都不好用，因此需要辅以特定的算法，比如 harmonic extraction 或者 shift-invert spectral transformation。这些算法在 SLEPc 中都有内置的支持。

本征值主算法的选项，由 -eps_type 指定，默认的是 krylovschur。同时 SLEPc 也提供一些 LAPACK solver，但这些 solver 只能以 dense matrix 和单进程上运行，因此不推荐。不同的算法适用于求不同部分的本征值和不同类型的矩阵，具体对比参看 manual。此外 EPS 模块也提供了一组控制 solver 误差和判断迭代终止的 API。对于迭代运算工作的子空间，EPS 也可以指定初始的子空间，或者不允许的子空间，后者被称为 deflate subspace，API 分别为 EPSSetInitialSpace，EPSSetDeflationSpace。

程序编译运行

PETSc 和 SLEPc 两个库的安装和其他库区别不大，但要注意先安装 PETSc 后安装 SLEPc，同时要注意两者的版本匹配问题。我就曾经因为安装了 master 分支 head 的 PETSc 库，从而发行版的 SLEPc 库无法安装的情况，报错稳定版与测试版不能匹配。此外需要注意指定 PETSC_DIR 和 SLEPC_DIR 的环境变量，用于指定安装位置。此外还可以指定 PETSC_ARCH 这一环境变量，其主要是用于区分不同版本，配置和编译器生成的 petsc 库 binary。即使不指定，安装也会默认生成一个该变量名，对应的就是 PETSC_DIR 安装版的文件夹名，如 arch-linux2-c-debug。对应的 SLEPc 的 PETSC_ARCH 变量必须一致。也就是 SLEPc 库和 PETSc 库是一一对应的关系。具体的安装调试过程，可以参考¹文章的附录。需要注意的时，最好提前配置好 cmake 和 mpi 等工具，否则问题也不大。PETSc 的 configure 工具非常智能，随时检查配置的兼容性，并且报错可读性强，按照提示来，很容易安装成功。

运行用户程序时，自然首先源代码要包含相应的头文件，此外需安装时的那三个环境变量 PETSC_DIR, PETSC_ARCH, SLEPC_DIR 还在即可正常的编译（利用 mpicc，但推荐使用下文介绍 PETSc 自己提供的 makefile 的方式，否则可能动态库链接的让人抓狂）和运行（利用 mpirun）。运行时加上 -log_view 的命令行选项，可以得到非常丰富的性能调试信息。此外 PestcViewer 用法也值得关注，其为程序提供了一个统一便捷的 IO 接口。

天河2号安装运行记录

首先天河通过 module load 的形式可以提供几个版本的 PETSc 和 SLEPc 的预装。但直接使用有一些问题，第一是很难保证 PETSc 和 SLEPc 的版本匹配，此外两者的版本都过于古老，已经有大量 API 发生了改变。最重要的是，这些预装的 PETSc 库都没有 configure 外部 solver，而这在 KSP 高效解线性方程组时是必须的，因此手动安装两个库是跑不掉的。由于天河的节点配置和其不能链接外网的特点（以前能用的 proxy 已经挂掉了），再加上登录节点上慢到令人发指的 configure 速度，使得顺利安装运行这两个库还是得费一番力气，下面是一个记录。

首先是将匹配版本的两个库的发行版压缩包 scp 到超算，解压 tar -zxvf *.tar.gz 后，先设置环境变量 PETSC_DIR，SLEPC_DIR 为两个库的安装路径，再设置 PETSC_ARCH 为任意名字作为安装版的区分。这些变量之后编译用户程序时还需要调用，建议统一存为 env.cfg 如下：

export PETSC_DIR=<petsc path>
export SLEPC_DIR=<slepc path>
export PETSC_ARCH=<branch name>

在正式 ./configure 配置 PETSc 之前，先激活环境变量 . env.cfg, 最好再 module load cmake/3.8.1 以利用新版本的 cmake 编译配置。这之后，还要将想要安装的底层 solver 库的压缩包先行下载，并 scp 到天河（因为直接下载不联网）。要注意这些外部库的下载版本问题，推荐的方式是去 PETSc 源码中的 config/BuildSystem/config/packages 文件夹下，寻找对应的外部库的配置 py 文件，并找到与对应版本 PETSc 兼容的外部 solver 的下载地址，也即对应文件的 Configure 对象的 self.download 属性。例如 PETSc v3.8 对应的 MUMPS 库的下载位置为：https://bitbucket.org/petsc/pkg-mumps/get/v5.1.1-p3.tar.gz。注意国内下载 bitbucket 内容极慢，建议科学上网。

将这些外部库存放在超算之后，配置 PETSc 的命令如下：

./configure  --with-pic=1 --with-shared-libraries=1 CFLAGS=-fPIC CXXFLAGS=-fPIC FFLAGS=-fPIC FCFLAGS=-fPIC F90FLAGS=-fPIC F77FLAGS=-fPIC --with-blas-lib=/WORK/app/BLAS/3.5.0-icc/lib/libblas.a --with-lapack-lib=/WORK/app/LAPACK/3.5.0-icc/lib/liblapack.a  --download-mumps=<mumps path>/mumps.tar.gz --download-metis=<metis path>/metis.tar.gz  --download-scalapack=<scalapack path>/scalapack.tgz

这段配置凝结了很多血泪教训，修改需谨慎，你可能需要配置 20 分钟才会收到报错。像什么 --with-64-bit-indices=1 和 MUMPS 安装冲突，BLAS 和 LAPACK 位置需要手动指定等坑真是踩都踩不完。当然有些没什么冲突的配置项也可以添加，比如关掉 debug 来提高效率 --with-debugging=0，和设置更多加速编译的 CFLAGS 等。配置完成后，直接 make all 就完成了 PETSc 库的用户级别的安装。

之后 SLEPc 库的安装就很容易，因为 configure 时会直接读取 PETSc 的配置，因此不需要额外配置。只需简单的 ./configure, make 即可完成安装。

之后是关于用户代码的编译和提交。首先需要注意这阶段遇到的问题，不一定是程序编译的锅。源码与安装库对应 API 不同，和天河残废的 c++ 11 支持，都可能导致报错，需要具体问题具体分析。而编译流程，则按照官方 manual 推荐的 makefile 方法，由于调用了库本身的 makefile，因此写起来非常轻松，makefile 如下：

CFLAGS = 
helloworld: helloworld.o chkopts
	-${CLINKER} -o hellworld helloworld.o ${SLEPC_SYS_LIB}
	
include ${SLEPC_DIR}/lib/slepc/conf/slepc_common

makefile 通过 include 语句调用库内置的 makefile 逻辑，从而自动处理编译过程，我们需要手动维护的只有安装时那三个环境变量即可，也即 make 之前，可以执行一遍 . env.cfg。生成可执行文件后，通过 slurm 提交任务和运行的方式就和普通程序一样了。

本征值问题黑魔法掠影

对于本征值问题 \(Ax=\lambda x\)，最简单的想法当然是直接求 \(\det(A-\lambda I)\) 的根，不过这种做法不仅速度无法接受，在数值上也不稳定。由于无论在科学上还是工程上，本征值问题都可能是最重要的一个数值问题（没有之一的那种），因此自然发展出了大量的本征值求解的数值算法。其方法多样，适用条件有别，在历史发展的长河中就形成了很多外行一脸蒙比的黑魔法和 dirty trick，这也是为什么这类问题的底层往往是大量的祖传 Fortran 老代码。本文无意事无巨细的讲解这些本征问题算法的细节，但鉴于求解本征值问题实践中，这些黑魔法暴露出了很多选项和参数到上层用户，因此对本征值求解算法有一个大概的了解还是必须的，不然面对那么多的选项和参数，可能就束手无策了。当然想高效的解决某种系统的本征值问题，需要大量的参数调整与尝试，可以说是一门和神经网络调超参一个级别的手艺了。

我们先从解决线性方程组的问题 \(Ax=b\) 看起，这涉及 PETSc 的 KSP 和 PC 模块，其中 PC 被封装在 KSP 中，往往不直接把对象暴露给用户。KSP 的 API 几乎和 SLEPc 中的 EPS API 完全一致。一般来讲，解决线性方程组问题，需要选取一个 -ksp_type 作为主算法和一个 -pc_type 参数作为 precondition 的算法。其中的 precondition 步骤在调用 KSPSetOperators() 时被计算，而主算法则在 KSPSolve() 中调用。

Krylov Methods

一类最主要的解线性方程组的算法就是 Krylov methods。其主要是在一个子空间中迭代，来寻求线性系统的近似解向量。该子空间被称为 Krylov subspace，该族算法的具体介绍可以参考⁵。其包含了一大类的具体算法，每种算法还有更多的参数可调。这主要通过 KSPSetType() 和 KSP[alg]Set[opt]() 的 API 来控制。

另一方面由于 PETSc 的高度可扩展性，我们还可以选取外部库提供的 linear solver，不过这需要在安装 PETSc 的时候就做好 configure， 如 ./cofigure --download-mumps 等。同时使用这些外部 solver 完全不需要更改代码，这需要提供 runtime 命令行选项即可，如 -pc_factor_mat_solver_type <packagename>。这里的 package 包含大量选项，比如 matlab, superlu, elemental, mkl_pardiso, mumps 等，其强大的胶水功能也体现了 PETSc 库的目标之一 Extensible。

preconditioning

由于 Krylov 方法的收敛速度依赖于矩阵 A 的谱的情况，因此有必要对矩阵进行一些预处理，使得 Krylov 方法收敛更快，这就是 preconditioning。具体的如下式

precondition 步骤就是通过某些算法构造矩阵 \(M_L,M_R\)，使得 \(M_L^{-1}AM_R^{-1}\) 更适合应用 Krylov 方法求解。其中如果 \(M_L(M_R)=I\)，则被称为 right(left) preconditioning. 这些算法都是被 PETSc 中的 PC 模块控制的。大部分 KSP 实现使用的是 left preconditioning。这可以通过 KSPSetPCSide() API 或 -ksp_pc_side left 来设置。注意和 KSP 情况相同，不同的 PC type 对应的算法也还有更多的参数可调。

本征值问题解决框架

大多数本征值的数值方法都利用了 Rayleigh-Ritz 投影，通过把问题限制在某个合适的子空间来求解。所以不同的 solver 的区别主要在于如何构建子空间，和如何加快收敛和减小内存使用。其中利用一个向量及其生成向量 \(A^n v\) 张成的子空间就回到了上文所述的 Krylov 方法，包括著名的 Arnoldi 和 Lanczos 算法。这些方法的收敛性和本征值谱的分布有关，越密的地方越难得到好的结果，因此用来计算极大和极小本征值效果更好。为了能够计算其他位置的本征值，通常需要一些谱变换和如上文所述的 preconditioning。precondition 可以看成廉价版的谱变换（主要是 invert），这些被称为 preconditioned eigensolvers，也就是说本征值求解中的 precondition 是和 solver 绑在一起的，指定对应的 solver 就自带 precondition， 比如 JD 和 GD。有时为了求得内部的本征值，需要先求的固定外侧本征值，并优化子空间重新计算，这通常被称为 restart。关于本征值问题和利用 SLEPc 解决该问题的介绍，可以参考⁶。

稀疏矩阵求解大量本征值的讨论

大量本征值是相对的概念，很多时候问题只需要我们求的一个或者几个本征值。一般来讲需要求得几千个本征值的问题，就算是大量本征值了。首先，这种情形下没必要再把 ncv 设置到 nev 的两倍以上，而是略大即可。关于这种问题 cost 的总结：

1. 空间上，至少要存储 ncv 个矩阵维度 n 的向量
2. 时间上，向量的 Schmit 正交化需要时间 \(O(ncv^2 n)\)
3. 空间上，至少一个 ncv 乘 ncv 的 dense matrix 需要被存储
4. 时间上，子空间本征值问题需要时间 \(O(ncv^3)\)，在迭代过程中这一本征值问题需要求解多次

其中 1，2两种 cost 在并行分布式环境中被平摊，问题不大。而在 SLEPc 中，3，4两步的子空间本征值问题是无法并行的，因此大量本征值求解中，子空间本征值求解会成为瓶颈。这一问题有一个 mpd 参数可以缓解，这里不再详述。由以上分析可以看出，稀疏矩阵对于求解全部本征值是没有意义的，不会有任何内存和时间的节约。

Harmonic extraction

为了克服求中间本征值收敛慢的问题，应用 Harmonic extraction 算法理论上先收敛到离 target 最近的本征值。其对应的选项是 -eps_harmonic。此时需要提高 ncv 的数值。

谱变换与本征值求解

ST 模块控制了谱变换和 precondition 等选项，作为计算本征值之前对矩阵的预处理。一般计算本征值时，只需在选项输入即可，不需额外改变代码。其通过 -st_type 指定，最常用的即为 sinvert 对应 shift-invert 谱变换。这些谱变换的参数和 EPSSetTarget 的参数密切相关。对于大多数谱变换，需要处理逆矩阵对向量的乘法，这通常通过内部的 KSP 模块以解线性方程的方式得到，而不会直接去求逆矩阵。为了修改这里对应的 KSP 的参数，只需在 KSP 和 PC 原参数前加 -st_ 即可，这就表示这些线性解模块是为了谱变换处理服务的。默认的选项是 -st_ksp_type preonly -st_pc_type lu。ST 模块包含的 KSP 对象，同样用来服务于具有 precondition 的 eigensolver，因为其 precondition 也涉及线性方程的解。在整个本征值求解问题中，谱变换是完全透明的，用户只需要指定对角化原始矩阵即可，得到的结果也自动会转换回原始矩阵。相关的模块在本征值问题下的依赖关系和逻辑如下图。

线性代数各库的关系

这部分简要说明一下一些常见的线性代数库的关系。因为新手往往被满天飞的和线性代数有关的库的名字吓到⁷。这里区分它们主要有以下几点：是底层库，还是中间层的 wrapper；能否处理稀疏矩阵；能否处理分布式内存问题等。

一般来讲线性代数计算库都是基于 BLAS (Basic Linear Algebra Subprograms) 和 LAPACK (Linear Algebra Package) 的，它们用来处理 dense matrix 的问题。BLAS 是一些矩阵乘法相关的 subroutine，其 code 可能非常偏向考虑各种可能的硬件优化和 cpu 加速。BLAS 实际上只是一些 API 规范，其实现非常之多。除了官方参考实现（非常简陋，性能不行）之外，还包括 OpenBLAS, GSL（GNU Scientific Library） BLAS，MKL BLAS（Intel 实现版本），ACML BLAS （AMD 版本），甚至 CUDA SDK 中也有对应的 BLAS API 实现。LAPACK 基于 BLAS API，解决一些更高层次的问题，包括矩阵的 LU，QR 等分解与求本征值，解线性方程等。同样的 LAPACK 也有诸如 MKL 等版本的实现。ATLAS (Automatically Tuned Linear Algebra Software) 则提供了 BLAS 和部分 LAPACK API 的实现，可以一定程度上作为一种 BLAS/LAPACK 的实现。以下几乎所有库，几乎全部严格依赖 BLAS 和 LAPACK 才能运行，反应了 BLAS/LAPACK 在科学计算中无可撼动的基础地位。

BLACS (Basic Linear Algebra Communication Subprograms) 是一组基于 MPI 的消息传播 API。该库作为底层，为其他分布式的线性代数计算奠定了基础。其主要用于 dense matrix 的通信。

ScaLAPACK (Scalable LAPACK) 则是建立在 BLACS 上的分布式的 LAPACK 库，其中还缺少了一些 LAPACK API 的支持。

以上这些库就属于 dense matrix 的基础设施，下面我们看一些 sparse matrix 的基础设施。

ARPACK (Arnoldi Package) 以 matrix free 的方式，通过 Arnoldi 算法来求解部分本征值。

SuperLU 提供了稀疏矩阵求解线性方程的方案。SuperLU_DIST 提供了分布式的版本。

此外还有 Elemental，可以支持分布式的 dense 和 sparse 矩阵计算。

这些线性代数的基础设施，往往 vendor 都会提供一些商业的 bundle，并且声称效率在相应硬件上更高，比如 Intel MKL。

除了这些还有其他很多提供比较具体算法和问题解决方案的库，它们可能提供某种 eigensolver 或线性系统的 solver，比如 MUMPS（Maultifrontal Massively Parallel sparse direct Solver）。

在这些基础设施之上，就是一些胶水和 wrapper 库，致力于向上提供更易用的 API，向下有统一的接口可以调用不同的 backend。这方面做的比较好的包括 armadillo （不支持分布式内存）和本文的 PETSc/SLEPc （支持分布式稀疏矩阵操作，并且后端极度灵活）。

当然更高层的，还有提供非 C/C++/Fortran 接口的库。最著名的莫过于 Python 中的 numpy/scipy，其中 numpy 提供了 dense matrix 的解决方案，而 scipy.sparse 提供了稀疏矩阵的解决方案。其底层可以以 LAPACK 为后端。

这里列举的和矩阵操作有关的库当然很不完整，但有了整个逻辑框架之后，再遇到其他库，可以按照这里的逻辑，将其放在合适的位置来认识和使用，从而不会陷入一堆令人混淆的名字的海洋里。

Reference