LCOV - code coverage report
Current view: top level - eps/impls/krylov/arnoldi - arnoldi.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 123 135 91.1 %
Date: 2024-03-29 00:34:52 Functions: 9 10 90.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    SLEPc eigensolver: "arnoldi"
      12             : 
      13             :    Method: Explicitly Restarted Arnoldi
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Arnoldi method with explicit restart and deflation.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] "Arnoldi Methods in SLEPc", SLEPc Technical Report STR-4,
      22             :            available at https://slepc.upv.es.
      23             : */
      24             : 
      25             : #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/
      26             : 
      27             : typedef struct {
      28             :   PetscBool delayed;
      29             : } EPS_ARNOLDI;
      30             : 
      31          29 : static PetscErrorCode EPSSetUp_Arnoldi(EPS eps)
      32             : {
      33          29 :   PetscFunctionBegin;
      34          29 :   EPSCheckDefinite(eps);
      35          29 :   PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
      36          29 :   PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
      37          29 :   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
      38          29 :   if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
      39          29 :   PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
      40          29 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_TWOSIDED);
      41             : 
      42          29 :   PetscCall(EPSAllocateSolution(eps,1));
      43          29 :   PetscCall(EPS_SetInnerProduct(eps));
      44          29 :   PetscCall(DSSetType(eps->ds,DSNHEP));
      45          29 :   if (eps->extraction==EPS_REFINED || eps->extraction==EPS_REFINED_HARMONIC) PetscCall(DSSetRefined(eps->ds,PETSC_TRUE));
      46          29 :   PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
      47          29 :   PetscCall(DSAllocate(eps->ds,eps->ncv+1));
      48          29 :   PetscFunctionReturn(PETSC_SUCCESS);
      49             : }
      50             : 
      51          29 : static PetscErrorCode EPSSolve_Arnoldi(EPS eps)
      52             : {
      53          29 :   PetscInt           k,nv,ld;
      54          29 :   Mat                U,Op,H;
      55          29 :   PetscScalar        *Harray;
      56          29 :   PetscReal          beta,gamma=1.0;
      57          29 :   PetscBool          breakdown,harmonic,refined;
      58          29 :   BVOrthogRefineType orthog_ref;
      59          29 :   EPS_ARNOLDI        *arnoldi = (EPS_ARNOLDI*)eps->data;
      60             : 
      61          29 :   PetscFunctionBegin;
      62          29 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
      63          29 :   PetscCall(DSGetRefined(eps->ds,&refined));
      64          29 :   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
      65          29 :   PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
      66             : 
      67             :   /* Get the starting Arnoldi vector */
      68          29 :   PetscCall(EPSGetStartVector(eps,0,NULL));
      69             : 
      70             :   /* Restart loop */
      71         853 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
      72         824 :     eps->its++;
      73             : 
      74             :     /* Compute an nv-step Arnoldi factorization */
      75         824 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
      76         824 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
      77         824 :     if (!arnoldi->delayed) {
      78         788 :       PetscCall(STGetOperator(eps->st,&Op));
      79         788 :       PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
      80         788 :       PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv,&nv,&beta,&breakdown));
      81         788 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
      82         788 :       PetscCall(STRestoreOperator(eps->st,&Op));
      83          36 :     } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
      84          18 :       PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
      85          18 :       PetscCall(EPSDelayedArnoldi1(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
      86          18 :       PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
      87             :     } else {
      88          18 :       PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
      89          18 :       PetscCall(EPSDelayedArnoldi(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
      90          18 :       PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
      91             :     }
      92         824 :     PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
      93         824 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
      94             : 
      95             :     /* Compute translation of Krylov decomposition if harmonic extraction used */
      96         824 :     if (harmonic) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,NULL,&gamma));
      97             : 
      98             :     /* Solve projected problem */
      99         824 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     100         824 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     101         824 :     PetscCall(DSUpdateExtraRow(eps->ds));
     102         824 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     103             : 
     104             :     /* Check convergence */
     105         824 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
     106         824 :     if (refined) {
     107          63 :       PetscCall(DSGetMat(eps->ds,DS_MAT_X,&U));
     108          63 :       PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+1));
     109          63 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&U));
     110          63 :       PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_FALSE,NULL,NULL));
     111             :     } else {
     112         761 :       PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
     113         761 :       PetscCall(BVMultInPlace(eps->V,U,eps->nconv,PetscMin(k+1,nv)));
     114         761 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
     115             :     }
     116         824 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     117         824 :     if (eps->reason == EPS_CONVERGED_ITERATING && breakdown) {
     118           0 :       PetscCall(PetscInfo(eps,"Breakdown in Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
     119           0 :       PetscCall(EPSGetStartVector(eps,k,&breakdown));
     120           0 :       if (breakdown) {
     121           0 :         eps->reason = EPS_DIVERGED_BREAKDOWN;
     122           0 :         PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
     123             :       }
     124             :     }
     125         824 :     eps->nconv = k;
     126         853 :     PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
     127             :   }
     128          29 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     129          29 :   PetscFunctionReturn(PETSC_SUCCESS);
     130             : }
     131             : 
     132          21 : static PetscErrorCode EPSSetFromOptions_Arnoldi(EPS eps,PetscOptionItems *PetscOptionsObject)
     133             : {
     134          21 :   PetscBool      set,val;
     135          21 :   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;
     136             : 
     137          21 :   PetscFunctionBegin;
     138          21 :   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Arnoldi Options");
     139             : 
     140          21 :     PetscCall(PetscOptionsBool("-eps_arnoldi_delayed","Use delayed reorthogonalization","EPSArnoldiSetDelayed",arnoldi->delayed,&val,&set));
     141          21 :     if (set) PetscCall(EPSArnoldiSetDelayed(eps,val));
     142             : 
     143          21 :   PetscOptionsHeadEnd();
     144          21 :   PetscFunctionReturn(PETSC_SUCCESS);
     145             : }
     146             : 
     147           4 : static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
     148             : {
     149           4 :   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
     150             : 
     151           4 :   PetscFunctionBegin;
     152           4 :   arnoldi->delayed = delayed;
     153           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     154             : }
     155             : 
     156             : /*@
     157             :    EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
     158             :    in the Arnoldi iteration.
     159             : 
     160             :    Logically Collective
     161             : 
     162             :    Input Parameters:
     163             : +  eps - the eigenproblem solver context
     164             : -  delayed - boolean flag
     165             : 
     166             :    Options Database Key:
     167             : .  -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi
     168             : 
     169             :    Note:
     170             :    Delayed reorthogonalization is an aggressive optimization for the Arnoldi
     171             :    eigensolver than may provide better scalability, but sometimes makes the
     172             :    solver converge less than the default algorithm.
     173             : 
     174             :    Level: advanced
     175             : 
     176             : .seealso: EPSArnoldiGetDelayed()
     177             : @*/
     178           4 : PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
     179             : {
     180           4 :   PetscFunctionBegin;
     181           4 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     182          16 :   PetscValidLogicalCollectiveBool(eps,delayed,2);
     183           4 :   PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
     184           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     185             : }
     186             : 
     187           5 : static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
     188             : {
     189           5 :   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
     190             : 
     191           5 :   PetscFunctionBegin;
     192           5 :   *delayed = arnoldi->delayed;
     193           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     194             : }
     195             : 
     196             : /*@
     197             :    EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
     198             :    iteration.
     199             : 
     200             :    Not Collective
     201             : 
     202             :    Input Parameter:
     203             : .  eps - the eigenproblem solver context
     204             : 
     205             :    Output Parameter:
     206             : .  delayed - boolean flag indicating if delayed reorthogonalization has been enabled
     207             : 
     208             :    Level: advanced
     209             : 
     210             : .seealso: EPSArnoldiSetDelayed()
     211             : @*/
     212           5 : PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
     213             : {
     214           5 :   PetscFunctionBegin;
     215           5 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     216           5 :   PetscAssertPointer(delayed,2);
     217           5 :   PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
     218           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     219             : }
     220             : 
     221          22 : static PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
     222             : {
     223          22 :   PetscFunctionBegin;
     224          22 :   PetscCall(PetscFree(eps->data));
     225          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL));
     226          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL));
     227          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     228             : }
     229             : 
     230           0 : static PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
     231             : {
     232           0 :   PetscBool      isascii;
     233           0 :   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;
     234             : 
     235           0 :   PetscFunctionBegin;
     236           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     237           0 :   if (isascii && arnoldi->delayed) PetscCall(PetscViewerASCIIPrintf(viewer,"  using delayed reorthogonalization\n"));
     238           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     239             : }
     240             : 
     241          22 : SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
     242             : {
     243          22 :   EPS_ARNOLDI    *ctx;
     244             : 
     245          22 :   PetscFunctionBegin;
     246          22 :   PetscCall(PetscNew(&ctx));
     247          22 :   eps->data = (void*)ctx;
     248             : 
     249          22 :   eps->useds = PETSC_TRUE;
     250             : 
     251          22 :   eps->ops->solve          = EPSSolve_Arnoldi;
     252          22 :   eps->ops->setup          = EPSSetUp_Arnoldi;
     253          22 :   eps->ops->setupsort      = EPSSetUpSort_Default;
     254          22 :   eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
     255          22 :   eps->ops->destroy        = EPSDestroy_Arnoldi;
     256          22 :   eps->ops->view           = EPSView_Arnoldi;
     257          22 :   eps->ops->backtransform  = EPSBackTransform_Default;
     258          22 :   eps->ops->computevectors = EPSComputeVectors_Schur;
     259             : 
     260          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi));
     261          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi));
     262          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     263             : }

Generated by: LCOV version 1.14