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