SLEPc

chimez

2021-07-03 13:08

源文件

SLEPc: PETSc 的本征值求解器扩展

SLEPc 是用来求 PETSc 矩阵的扩展库。

基本结构

SLEPc 有以下基本部分构成

EPS (Eigenvalue Problem Solver)
SVD (Singular Value Decomposition)
PEP (Polynomial Eigenvalue Problem)
NEP (Nonlinear Eigenvalue Problem)
MFN (Matrix Function)

还有一些辅助模块

ST (Spectral Transformation)
BV (Basis Vectors)
DS (Dense System)
RG (Region)
FN (Mathematical Function)

安装

下载时注意，SLEPc 的版本必须和 PETSc 版本一样

编译选项

   ./configure --download-arpack --download-blopex --download-hpddm --download-primme  --with-slepc4py=1

最后别忘了把环境变量 SLEPC_DIR 加上

基本使用

程序开始和结束时分别调用 SlepcInitialize() 和 SlepcFinalize(), 实际上就是 PetscInitialize() 和 PetscFinalize()

矩阵的使用跟 PETSc 完全一致

EPS: 矩阵本征值求解器

基本使用

EPSCreate(MPI_Comm comm,EPS *eps);
1. 创建求解器
EPSSetOperators(EPS eps,Mat A,Mat B);
1. 设置要解的算符
EPSSetFromOptions(EPS eps);
1. 添加运行时设置
EPSSolve(EPS eps);
1. 运行求解器
EPSDestroy(EPS *eps);
1. 销毁求解器
EPSSetUp(EPS eps);
1. 某些更具体的设置, 如果有需要的话，要在运行求解器之前设置
EPSGetST(EPS eps,ST *st);
1. 获得谱变换器(ST), EPS 会默认自动创建一个 ST 对象，如果需要设置谱变换，就要这样获得它
EPSView(EPS eps,PetscViewer viewer);
1. 打印求解器配置

    #include <slepceps.h>


    Mat A; /* problem matrix */
    /* do something to fill A */

    PetscErrorCode ierr;
    /* Create eigensolver context */
    EPS eps; /* eigenproblem solver context */
    ierr = EPSCreate(PETSC_COMM_WORLD, &eps);
    CHKERRQ(ierr);

    /* Set operators. In this case, it is a standard eigenvalue problem */
    ierr = EPSSetOperators(eps, A, NULL);
    CHKERRQ(ierr);
    ierr = EPSSetProblemType(eps, EPS_HEP);
    CHKERRQ(ierr);

    /* show EPS setup */
    ierr = EPSView(eps, PETSC_VIEWER_STDOUT_WORLD) ;
    CHKERRQ(ierr);

    /* Set solver parameters at runtime */
    ierr = EPSSetFromOptions(eps);
    CHKERRQ(ierr);

    /* Solve the eigensystem */
    ierr = EPSSolve(eps);
    CHKERRQ(ierr);

    /* Optional: Get some information from the solver and display it */
    PetscInt its;
    ierr = EPSGetIterationNumber(eps, &its);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Number of iterations of the method: %D\n", its);
    CHKERRQ(ierr);

    EPSType type;
    ierr = EPSGetType(eps, &type);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Solution method: %s\n\n", type);
    CHKERRQ(ierr);

    PetscInt nev;
    ierr = EPSGetDimensions(eps, &nev, NULL, NULL);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Number of requested eigenvalues: %D\n", nev);
    CHKERRQ(ierr);

    PetscReal tol;
    PetscInt maxit;
    ierr = EPSGetTolerances(eps, &tol, &maxit);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Stopping condition: tol=%.4g, maxit=%D\n", (double)tol, maxit);
    CHKERRQ(ierr);

    /* Display solution and clean up */
    /* Get number of converged approximate eigenpairs */
    PetscInt nconv;
    ierr = EPSGetConverged(eps, &nconv);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Number of converged eigenpairs: %D\n\n", nconv);
    CHKERRQ(ierr);

    if (nconv > 0) {
    /* Display eigenvalues and relative errors */
    ierr = PetscPrintf(PETSC_COMM_WORLD,
                        "           k          ||Ax-kx||/||kx||\n"
                        "   ----------------- ------------------\n");
    CHKERRQ(ierr);

    for (i = 0; i < nconv; i++) {
        /* Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and ki (imaginary part) */
        PetscScalar kr, ki;
        Vec xr, xi;
        ierr = MatCreateVecs(A, NULL, &xr);
        CHKERRQ(ierr);
        ierr = MatCreateVecs(A, NULL, &xi);
        CHKERRQ(ierr);
        ierr = EPSGetEigenpair(eps, i, &kr, &ki, xr, xi);
        CHKERRQ(ierr);

        /* Compute the relative error associated to each eigenpair */
        PetscReal error;
        ierr = EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error);
        CHKERRQ(ierr);

        PetscReal re = PetscRealPart(kr);
        PetscReal im = PetscImaginaryPart(kr);
        if (im != 0.0) {
            ierr = PetscPrintf(PETSC_COMM_WORLD, " %9f%+9fi %12g\n", (double)re, (double)im, (double)error);
            CHKERRQ(ierr);
        } else {
            ierr = PetscPrintf(PETSC_COMM_WORLD, "   %12f       %12g\n", (double)re, (double)error);
            CHKERRQ(ierr);
        }
    }
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");
    CHKERRQ(ierr);
    }


    /* Free work space */

    ierr = EPSDestroy(&eps);
    CHKERRQ(ierr);
    ierr = MatDestroy(&A);
    CHKERRQ(ierr);
    ierr = VecDestroy(&xr);
    CHKERRQ(ierr);
    ierr = VecDestroy(&xi);
    CHKERRQ(ierr);

求解器配置

问题类型

问题类型使用 EPSSetProblemType(EPS eps,EPSProblemType type); 进行设置

可选的有

问题类型	EPSProblemTyle	命令行选项
Hermitian	EPS_HEP	-eps_hermitian
Non-Hermitian	EPS_NHEP	-eps_non_hermitian
Generalized Hermitian	EPS_GHEP	-eps_gen_hermitian
Generalized Hermitian indefinite	EPS_GHIEP	-eps_gen_indefinite
Generalized Non-Hermitian	EPS_GNHEP	-eps_gen_non_hermitian
GNHEP with positive (semi-)definite	BEPS_PGNHEP	-eps_pos_gen_non_hermitian

还可以用一系列函数来判断设置的问题类型

EPSIsGeneralized(EPS eps,PetscBool *gen);
EPSIsHermitian(EPS eps,PetscBool *her);
EPSIsPositive(EPS eps,PetscBool *pos);

求解本征值个数
使用 EPSSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd); 设置要求的本征值个数
1. PetscInt nev 本征值个数
2. PetscInt ncv 最大工作空间的维数，指的是用多少个中间向量, ncv 至少和 nev 一样多，最好是它的两倍以上
3. PetscInt mpd 最大投影空间维数，用来计算很多本征值个数的时候，设置它可以减少工作空间的需求

求解本征值位置

EPSSetWhichEigenpairs(EPS eps,EPSWhich which); 设置要求本征值的位置
EPSSetTarget(EPS eps,PetscScalar target); 设置计算距离某个值最近的本征值
EPSSetInterval(EPS eps,PetscScalar a,PetscScalar b); 计算 \(\lambda \in [a, b]\) 的所有本征值
EPSSetTwoSided(EPS eps,PetscBool twosided); 计算非厄米问题的左本征矢

可选的本征值位置有

EPSWhich	命令行	说明
EPS_LARGEST_MAGNITUDE	-eps_largest_magnitude	最大 \(\vert \lambda\vert\)
EPS_SMALLEST_MAGNITUDE	-eps_smallest_magnitude	最小 \(\vert\lambda\vert\)
EPS_LARGEST_REAL	-eps_largest_real	最大 \(\mathrm{Re}(\lambda)\)
EPS_SMALLEST_REAL	-eps_smallest_real	最小 \(\mathrm{Re}(\lambda)\)
EPS_LARGEST_IMAGINARY	-eps_largest_imaginary	最大 \(\mathrm{Im}(\lambda)\)
EPS_SMALLEST_IMAGINARY	-eps_smallest_imaginary	最小 \(\mathrm{Im}(\lambda)\)
EPS_TARGET_MAGNITUDE	-eps_target_magnitude	最小 \(\vert\lambda - \tau\vert\)
EPS_TARGET_REAL	-eps_target_real	最小 \(\vert\mathrm{Re}(\lambda-\tau)\vert\)
EPS_TARGET_IMAGINARY	-eps_target_imaginary	最小 \(\vert\mathrm{Im}(\lambda-\tau)\vert\)
EPS_ALL	-eps_all	所有 \(\lambda \in [a,b]\)
EPS_WHICH_USER		用户定义

选择算法

通过 EPSSetType(EPS eps,EPSType method); 设置求解器使用的算法

SLEPc 支持的算法有

基本算法
- Power Itration, Rayleigh Quotient iteration(RQI)
- Subspace Iteration with Rayleigh-Ritz projection and locking
- Arnoldi method with explicit restart and deflation
- Lanczos with explicit restart
Krylov-Schur, thick-restart Lanczos method (默认)
Generalized Davidson
Jacobi-Davidson
RQCG
LOBPCG
CISS
Lyapunov inverse iteration

注意：

由于实现支持有限，只有 arnoldi, krylov-schur, gd, jd, arpack 支持所有类型的问题求解

方法	EPSType	选项名
Power / Inverse / RQI	EPSPOWER	power
Subspace Iteration	EPSSUBSPACE	subspace
Arnoldi	EPSARNOLDI	arnoldi
Lanczos	EPSLANCZOS	lanczos
Krylov-Schur	EPSKRYLOVSCHUR	krylovschur
Generalized Davidson	EPSGD	gd
Jacobi-Davidson	EPSJD	jd
Rayleigh quotient CG	EPSRQCG	rqcg
LOBPCG	EPSLOBPCG	lobpcg
Contour integral SS	EPSCISS	ciss
Lyapunov Inverse Iteration	EPSLYAPII	lyapii
LAPACK solver	EPSLAPACK	lapack
Wrapper to arpack	EPSARPACK	arpack
Wrapper to primme	EPSPRIMME	primme
Wrapper to evsl	EPSEVSL	evsl
Wrapper to trlan	EPSTRLAN	trlan
Wrapper to blopex	EPSBLOPEX	blopex
Wrapper to scalapack	EPSSCALAPACK	scalapack
Wrapper to elpa	EPSELPA	elpa
Wrapper to elemental	EPSELEMENTAL	elemental
Wrapper to feast	EPSFEAST	feast

获得结果
1. EPSGetConverged(EPS eps,PetscInt *nconv); 获得收敛的解个数
2. EPSGetEigenpair(EPS eps,PetscInt j,PetscScalar *kr,PetscScalar *ki,Vec xr,Vec xi); 获得本征值和本征矢
3. EPSGetLeftEigenvector(EPS eps,PetscInt j,Vec yr,Vec yi); 获得左本征矢，如果设置求解器计算左本征矢的话
4. EPSComputeError(EPS eps,PetscInt j,EPSErrorType type,PetscReal *error); 获得结果的误差
  1. 可选的误差类型有 \(||r||\) EPS_ERROR_ABSOLUTE, \(||r||/|\lambda|\) EPS_ERROR_RELATIVE, \(||r||/(||A||+|\lambda| ||B||)\) EPS_ERROR_BACKWARD
5. EPSGetIterationNumber(EPS eps,PetscInt *its); 获得迭代次数
6. EPSSetTolerances(EPS eps,PetscReal tol,PetscInt max_it); 设置误差和最大迭代次数
  1. 默认误差 \(10^{-8}\)
7. EPSSetTrueResidual(EPS eps,PetscBool trueres); 设置使用真实残差计算收敛条件，而不是默认的简化形式
8. EPSSetConvergenceTest(EPS eps,EPSConv conv); 设置收敛条件
  1. 可选的误差类型有 \(||r||\) EPS_CONV_ABS, \(||r||/|\lambda|\) EPS_CONV_REL, \(||r||/(||A||+|\lambda| ||B||)\) EPS_CONV_NORM
9. EPSGetErrorEstimate(EPS eps,PetscInt j,PetscReal *errest); 获得误差估计
注意：
1. 如果编译时使用实数模式，那么 kr, ki 等分别是实部和虚部
2. 如果编译时使用复数模式，那么 kr, xr 中就保存复数结果， ki,xi 不使用（设为全零）
3. 如果设置命令行选项 -eps_monitor 那么会在每次迭代过程中打印计算过程
4. 其它命令行选项甚至可以画图（ -eps_monitor draw:draw_lg -draw_pause .2 ）
其它设置
- 初始猜解
  
  EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[]);
- 处理简并
- 选择正交化方法
- 选择滤波算法
- 计算大量本征值的方法

SLEPc: PETSc 的本征值求解器扩展

基本结构

安装

基本使用

EPS: 矩阵本征值求解器

基本使用

求解器配置

评论