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.

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

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 *,兼容性更强.例如本书的 SLISC0 库用 OpenBlas 就会编译出错(待更新).


致读者: 小时百科一直以来坚持所有内容免费,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利