Arpack++2 大型本征方程库

                     

贡献者： addis

　　 Arpack（ARnoldi PACKage）最初是一个 Fortran 语言编写的，用于求解大型本征方程的程序库，其主要算法是 Arnoldi 循环。注意如果要求密矩阵的本征问题，那么没必要使用 Arpack，直接使用 Lapack 即可。下面是 Arpack 的一些特性：

• 可以求解广义本征问题（General）$ \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{x}} = \lambda \boldsymbol{\mathbf{B}} \boldsymbol{\mathbf{x}} $ 或者标准本征问题（Standard）$ \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{x}} = \lambda \boldsymbol{\mathbf{x}} $。
• 矩阵的数据类型支持密矩阵，系数的 CSC 矩阵。也可以不需要提供矩阵本身而是提供矩阵—矢量乘法。也就是你只需要提供一个函数 v = F(u) 给 Arpack 即可，其中 v 就是矩阵乘以 u 的结果。
• Regular 模式：它可以求出本征值最小（大）的指定个数的本征矢，或者本征值的绝对值最小（大）的指定个数的本征矢。
• Shift 模式：可以求出指定值附近的若干本征值和本征矢。但这可能需要提供矩阵的逆（见下文）。

　　 Arpack++ 是 Arpack 的 C++ 模板库，但实际上只是一个 C++ 接口（只提供头文件），仍然调用 Fortran 版的动态链接库。官网提供的版本已经多年没有维护（最后更新于 2000 年，版本 1.2），用目前的 g++ 编译器无法正常编译。麻省理工在 GitHub 上维护了一个新版本：Arpack++2。本文使用的是当前最新的 release 2.3。另见 github 上的 Arpack-ng，这个应该是原 Fortran 版本的更新，Arpack++2 依赖于它，ubuntu 中也可以用 sudo apt install libarpack2-dev 安装。Arpack++2 只包含头文件，需要依赖于 fortran 的 libarpack2-dev。

　　 另一个由 Arpack 改写的 C++ 头文件库是 Spectra，基于 Eigen。

编译 arpack-ng

• 详见 SLISC 中的 *.dockerfile
• 安装命令：wget https://github.com/opencollab/arpack-ng/archive/refs/tags/3.8.0.tar.gz，tar -xzf 3.8.0.tar.gz，cd arpack-ng-3.8.0，sh bootstrap，./configure，make -j12，make check，make install。
• 安装的文件：/usr/local/include/arpack/ 中的 arpackdef.h，debug.h，stat.h。库文件：/usr/local/lib/ 中的 libarpack.la，libarpack.so，libarpack.so.2，libarpack.so.2.1.0

1. Ubuntu/Debian 安装

　　 Ubuntu/Debian 中可以直接用 apt 安装，如果没有装 gfortran 要先 apt install gfortran。然后安装 apt install libarpack++2-dev（当前最新是 2.3 版），其实我们只需要它安装的二进制文件和 dependency。注意 dependency 中会自动安装 OpenBlas，如果已经装了其他版本的 Blas 如 CBLAS 可以将其手动卸载¹：sudo apt purge libopenblas-dev。

　　 安装好后，可以用 dpkg -L libarpack++2-dev 查看安装的文件，其中头文件和文档可以删除，使用 github 上的头文件和文档。doc 中的 pdf 文档是重要的参考。例如要编译 examples/areig/nonsym/simple.cc 中的例子，在该目录运行 make simple 即可编译。用 ./simple 运行。

2. 使用

　　 第三章讲解了 Classes that require user-defined matrix-vector products，即不直接提供矩阵而是提供矩阵—矢量乘法。Real symmetric standard problems 讲解了具体怎么用 ARSymStdEig 类。如果要求解（数值或绝对值）最大或者最小的本征值，那么只需要矩阵—矢量乘法即可，但如果要算某个值 $\sigma$ 附近的本征值，那么需要提供 $(A - \sigma)^{-1}$ 的乘法。这有时候会大大增加难度，甚至不可能做到。

　　 Arpack 解决的本征问题分为两种，一种是 General 本征问题，就是 $Ax = \lambda Bx$，例如例程 symg...cc 代表对称矩阵的 General 本征问题，没有 g 代表是 Standard 本征问题。

　　 example/product/sym 文件夹中是通过提供矩阵—矢量乘法来解对称矩阵的方程。其中没有 g 的例子有两个，分别是 symreg.cc 和 symshft.cc，使用的都是 ARSymStdEig 类。其中 Sym 代表对称，Std 代表标准本征值问题。使用方法：

#include "arssym.h"
// 求 4 个绝对值最小的本征值的本征矢。
ARSymStdEig<T, SymMatrixA<T> >
    dprob(A.ncols(), 4, &A, &SymMatrixA<T>::MultMv, "SM");
// 求 4 个最接近 1 的本征值的本征矢。
ARSymStdEig<T, SymMatrixB<T> >
    dprob(A.ncols(), 4, &A, &SymMatrixB<T>::MultOPv, 1.);

其中 SymMatrixB<ART>::MultMv(T* u, T* v) 提供 $ \boldsymbol{\mathbf{v}} = \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{u}} $。MultOPv(T* u, T*v) 提供 $ \boldsymbol{\mathbf{v}} = ( \boldsymbol{\mathbf{A}} - \sigma)^{-1} \boldsymbol{\mathbf{u}} $，例子中 $\sigma = 1$。

3. 量子力学应用：波函数本征态

　　 根据笔者的经验，除非是一维问题，否则企图直接用 Arpack 在网格基底下解出一系列波函数的本征态的成功率很小，大部分情况下很难收敛。一种方法是使用更物理一些的基底，例如氢原子束缚态作为基底等，然后算出哈密顿矩阵 $ \left\langle i \middle| H \middle| j \right\rangle $，再用 Lapack 求解。

---

1. ^ OpenBlas 和 CBLAS 的接口不完全相同，前者的函数使用 double * 输入复数，而后者使用 void *，兼容性更强。