贡献者: addis
Arpack(ARnoldi PACKage)最初是一个 Fortran 语言编写的,用于求解大型本征方程的程序库,其主要算法是 Arnoldi 循环。注意如果要求密矩阵的本征问题,那么没必要使用 Arpack,直接使用 Lapack 即可。下面是 Arpack 的一些特性:
v = F(u)
给 Arpack 即可,其中 v
就是矩阵乘以 u
的结果。
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。
*.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
Ubuntu/Debian 中可以直接用 apt 安装,如果没有装 gfortran 要先 apt install gfortran
。然后安装 apt install libarpack++2-dev
(当前最新是 2.3 版),其实我们只需要它安装的二进制文件和 dependency。注意 dependency 中会自动安装 OpenBlas,如果已经装了其他版本的 Blas 如 CBLAS 可以将其手动卸载1:sudo apt purge libopenblas-dev
。
安装好后,可以用 dpkg -L libarpack++2-dev
查看安装的文件,其中头文件和文档可以删除,使用 github 上的头文件和文档。doc 中的 pdf 文档是重要的参考。例如要编译 examples/areig/nonsym/simple.cc
中的例子,在该目录运行 make simple
即可编译。用 ./simple
运行。
第三章讲解了 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$。
根据笔者的经验,除非是一维问题,否则企图直接用 Arpack 在网格基底下解出一系列波函数的本征态的成功率很小,大部分情况下很难收敛。一种方法是使用更物理一些的基底,例如氢原子束缚态作为基底等,然后算出哈密顿矩阵 $ \left\langle i \middle| H \middle| j \right\rangle $,再用 Lapack 求解。
1. ^ OpenBlas 和 CBLAS 的接口不完全相同,前者的函数使用 double *
输入复数,而后者使用 void *
,兼容性更强。