Actual source code: arnoldi.c

slepc-3.15.0 2021-03-31
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain

  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"

 13:    Method: Explicitly Restarted Arnoldi

 15:    Algorithm:

 17:        Arnoldi method with explicit restart and deflation.

 19:    References:

 21:        [1] "Arnoldi Methods in SLEPc", SLEPc Technical Report STR-4,
 22:            available at https://slepc.upv.es.
 23: */

 25: #include <slepc/private/epsimpl.h>

 27: typedef struct {
 28:   PetscBool delayed;
 29: } EPS_ARNOLDI;

 31: PetscErrorCode EPSSetUp_Arnoldi(EPS eps)
 32: {

 36:   EPSCheckDefinite(eps);
 37:   EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
 38:   if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
 39:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 40:   if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
 41:   if (eps->which==EPS_ALL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
 42:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_TWOSIDED);

 44:   EPSAllocateSolution(eps,1);
 45:   EPS_SetInnerProduct(eps);
 46:   DSSetType(eps->ds,DSNHEP);
 47:   if (eps->extraction==EPS_REFINED || eps->extraction==EPS_REFINED_HARMONIC) {
 48:     DSSetRefined(eps->ds,PETSC_TRUE);
 49:   }
 50:   DSSetExtraRow(eps->ds,PETSC_TRUE);
 51:   DSAllocate(eps->ds,eps->ncv+1);
 52:   return(0);
 53: }

 55: PetscErrorCode EPSSolve_Arnoldi(EPS eps)
 56: {
 57:   PetscErrorCode     ierr;
 58:   PetscInt           k,nv,ld;
 59:   Mat                U,Op;
 60:   PetscScalar        *H;
 61:   PetscReal          beta,gamma=1.0;
 62:   PetscBool          breakdown,harmonic,refined;
 63:   BVOrthogRefineType orthog_ref;
 64:   EPS_ARNOLDI        *arnoldi = (EPS_ARNOLDI*)eps->data;

 67:   DSGetLeadingDimension(eps->ds,&ld);
 68:   DSGetRefined(eps->ds,&refined);
 69:   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
 70:   BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL);

 72:   /* Get the starting Arnoldi vector */
 73:   EPSGetStartVector(eps,0,NULL);

 75:   /* Restart loop */
 76:   while (eps->reason == EPS_CONVERGED_ITERATING) {
 77:     eps->its++;

 79:     /* Compute an nv-step Arnoldi factorization */
 80:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 81:     DSSetDimensions(eps->ds,nv,0,eps->nconv,0);
 82:     DSGetArray(eps->ds,DS_MAT_A,&H);
 83:     if (!arnoldi->delayed) {
 84:       STGetOperator(eps->st,&Op);
 85:       BVMatArnoldi(eps->V,Op,H,ld,eps->nconv,&nv,&beta,&breakdown);
 86:       STRestoreOperator(eps->st,&Op);
 87:     } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
 88:       EPSDelayedArnoldi1(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
 89:     } else {
 90:       EPSDelayedArnoldi(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
 91:     }
 92:     DSRestoreArray(eps->ds,DS_MAT_A,&H);
 93:     DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
 94:     BVSetActiveColumns(eps->V,eps->nconv,nv);

 96:     /* Compute translation of Krylov decomposition if harmonic extraction used */
 97:     if (harmonic) {
 98:       DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,NULL,&gamma);
 99:     }

101:     /* Solve projected problem */
102:     DSSolve(eps->ds,eps->eigr,eps->eigi);
103:     DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
104:     DSUpdateExtraRow(eps->ds);
105:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

107:     /* Check convergence */
108:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
109:     if (refined) {
110:       DSGetMat(eps->ds,DS_MAT_X,&U);
111:       BVMultInPlace(eps->V,U,eps->nconv,k+1);
112:       MatDestroy(&U);
113:       BVOrthonormalizeColumn(eps->V,k,PETSC_FALSE,NULL,NULL);
114:     } else {
115:       DSGetMat(eps->ds,DS_MAT_Q,&U);
116:       BVMultInPlace(eps->V,U,eps->nconv,PetscMin(k+1,nv));
117:       MatDestroy(&U);
118:     }
119:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
120:     if (eps->reason == EPS_CONVERGED_ITERATING && breakdown) {
121:       PetscInfo2(eps,"Breakdown in Arnoldi method (it=%D norm=%g)\n",eps->its,(double)beta);
122:       EPSGetStartVector(eps,k,&breakdown);
123:       if (breakdown) {
124:         eps->reason = EPS_DIVERGED_BREAKDOWN;
125:         PetscInfo(eps,"Unable to generate more start vectors\n");
126:       }
127:     }
128:     eps->nconv = k;
129:     EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv);
130:   }
131:   DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);
132:   return(0);
133: }

135: PetscErrorCode EPSSetFromOptions_Arnoldi(PetscOptionItems *PetscOptionsObject,EPS eps)
136: {
138:   PetscBool      set,val;
139:   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;

142:   PetscOptionsHead(PetscOptionsObject,"EPS Arnoldi Options");

144:     PetscOptionsBool("-eps_arnoldi_delayed","Use delayed reorthogonalization","EPSArnoldiSetDelayed",arnoldi->delayed,&val,&set);
145:     if (set) { EPSArnoldiSetDelayed(eps,val); }

147:   PetscOptionsTail();
148:   return(0);
149: }

151: static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
152: {
153:   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;

156:   arnoldi->delayed = delayed;
157:   return(0);
158: }

160: /*@
161:    EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
162:    in the Arnoldi iteration.

164:    Logically Collective on eps

166:    Input Parameters:
167: +  eps - the eigenproblem solver context
168: -  delayed - boolean flag

170:    Options Database Key:
171: .  -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi

173:    Note:
174:    Delayed reorthogonalization is an aggressive optimization for the Arnoldi
175:    eigensolver than may provide better scalability, but sometimes makes the
176:    solver converge less than the default algorithm.

178:    Level: advanced

180: .seealso: EPSArnoldiGetDelayed()
181: @*/
182: PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
183: {

189:   PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
190:   return(0);
191: }

193: static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
194: {
195:   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;

198:   *delayed = arnoldi->delayed;
199:   return(0);
200: }

202: /*@
203:    EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
204:    iteration.

206:    Not Collective

208:    Input Parameter:
209: .  eps - the eigenproblem solver context

211:    Input Parameter:
212: .  delayed - boolean flag indicating if delayed reorthogonalization has been enabled

214:    Level: advanced

216: .seealso: EPSArnoldiSetDelayed()
217: @*/
218: PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
219: {

225:   PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
226:   return(0);
227: }

229: PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
230: {

234:   PetscFree(eps->data);
235:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL);
236:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL);
237:   return(0);
238: }

240: PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
241: {
243:   PetscBool      isascii;
244:   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;

247:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
248:   if (isascii && arnoldi->delayed) {
249:     PetscViewerASCIIPrintf(viewer,"  using delayed reorthogonalization\n");
250:   }
251:   return(0);
252: }

254: SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
255: {
256:   EPS_ARNOLDI    *ctx;

260:   PetscNewLog(eps,&ctx);
261:   eps->data = (void*)ctx;

263:   eps->useds = PETSC_TRUE;

265:   eps->ops->solve          = EPSSolve_Arnoldi;
266:   eps->ops->setup          = EPSSetUp_Arnoldi;
267:   eps->ops->setupsort      = EPSSetUpSort_Default;
268:   eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
269:   eps->ops->destroy        = EPSDestroy_Arnoldi;
270:   eps->ops->view           = EPSView_Arnoldi;
271:   eps->ops->backtransform  = EPSBackTransform_Default;
272:   eps->ops->computevectors = EPSComputeVectors_Schur;

274:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi);
275:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi);
276:   return(0);
277: }