Arpack++2 大型本征方程库

                     

贡献者: addis

预备知识 C++ 基础,密矩阵,稀疏矩阵

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

   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

1. Ubuntu/Debian 安装

   Ubuntu/Debian 中可以直接用 apt 安装,如果没有装 gfortran 要先 apt install gfortran。然后安装 apt install libarpack++2-dev(当前最新是 2.3 版),其实我们只需要它安装的二进制文件和 dependency。注意 dependency 中会自动安装 OpenBlas,如果已经装了其他版本的 Blas 如 CBLAS 可以将其手动卸载1sudo 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.ccsymshft.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 *,兼容性更强。

                     

© 小时科技 保留一切权利