Actual source code: arnoldi.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 http://slepc.upv.es.
 23: */

 25: #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/

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

 31: PetscErrorCode EPSSetUp_Arnoldi(EPS eps)
 32: {

 36:   EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
 37:   if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
 38:   if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 39:   if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
 40:   if (eps->which==EPS_ALL || (eps->ishermitian && eps->ispositive && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY))) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");

 42:   if (!eps->extraction) {
 43:     EPSSetExtraction(eps,EPS_RITZ);
 44:   }
 45:   if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");

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

 56:   if (eps->isgeneralized && eps->ishermitian && !eps->ispositive) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Requested method does not work for indefinite problems");
 57:   return(0);
 58: }

 60: PetscErrorCode EPSSolve_Arnoldi(EPS eps)
 61: {
 62:   PetscErrorCode     ierr;
 63:   PetscInt           k,nv,ld;
 64:   Mat                U;
 65:   PetscScalar        *H;
 66:   PetscReal          beta,gamma=1.0;
 67:   PetscBool          breakdown,harmonic,refined;
 68:   BVOrthogRefineType orthog_ref;
 69:   EPS_ARNOLDI        *arnoldi = (EPS_ARNOLDI*)eps->data;

 72:   DSGetLeadingDimension(eps->ds,&ld);
 73:   DSGetRefined(eps->ds,&refined);
 74:   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
 75:   BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL);

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

 80:   /* Restart loop */
 81:   while (eps->reason == EPS_CONVERGED_ITERATING) {
 82:     eps->its++;

 84:     /* Compute an nv-step Arnoldi factorization */
 85:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 86:     DSSetDimensions(eps->ds,nv,0,eps->nconv,0);
 87:     DSGetArray(eps->ds,DS_MAT_A,&H);
 88:     if (!arnoldi->delayed) {
 89:       EPSBasicArnoldi(eps,PETSC_FALSE,H,ld,eps->nconv,&nv,&beta,&breakdown);
 90:     } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
 91:       EPSDelayedArnoldi1(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
 92:     } else {
 93:       EPSDelayedArnoldi(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
 94:     }
 95:     DSRestoreArray(eps->ds,DS_MAT_A,&H);
 96:     DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
 97:     BVSetActiveColumns(eps->V,eps->nconv,nv);

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

104:     /* Solve projected problem */
105:     DSSolve(eps->ds,eps->eigr,eps->eigi);
106:     DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
107:     DSUpdateExtraRow(eps->ds);
108:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

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

135:   /* truncate Schur decomposition and change the state to raw so that
136:      DSVectors() computes eigenvectors from scratch */
137:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
138:   DSSetState(eps->ds,DS_STATE_RAW);
139:   return(0);
140: }

142: PetscErrorCode EPSSetFromOptions_Arnoldi(PetscOptionItems *PetscOptionsObject,EPS eps)
143: {
145:   PetscBool      set,val;
146:   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;

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

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

154:   PetscOptionsTail();
155:   return(0);
156: }

158: static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
159: {
160:   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;

163:   arnoldi->delayed = delayed;
164:   return(0);
165: }

167: /*@
168:    EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
169:    in the Arnoldi iteration.

171:    Logically Collective on EPS

173:    Input Parameters:
174: +  eps - the eigenproblem solver context
175: -  delayed - boolean flag

177:    Options Database Key:
178: .  -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi

180:    Note:
181:    Delayed reorthogonalization is an aggressive optimization for the Arnoldi
182:    eigensolver than may provide better scalability, but sometimes makes the
183:    solver converge less than the default algorithm.

185:    Level: advanced

187: .seealso: EPSArnoldiGetDelayed()
188: @*/
189: PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
190: {

196:   PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
197:   return(0);
198: }

200: static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
201: {
202:   EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;

205:   *delayed = arnoldi->delayed;
206:   return(0);
207: }

209: /*@
210:    EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
211:    iteration.

213:    Not Collective

215:    Input Parameter:
216: .  eps - the eigenproblem solver context

218:    Input Parameter:
219: .  delayed - boolean flag indicating if delayed reorthogonalization has been enabled

221:    Level: advanced

223: .seealso: EPSArnoldiSetDelayed()
224: @*/
225: PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
226: {

232:   PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
233:   return(0);
234: }

236: PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
237: {

241:   PetscFree(eps->data);
242:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL);
243:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL);
244:   return(0);
245: }

247: PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
248: {
250:   PetscBool      isascii;
251:   EPS_ARNOLDI    *arnoldi = (EPS_ARNOLDI*)eps->data;

254:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
255:   if (isascii && arnoldi->delayed) {
256:     PetscViewerASCIIPrintf(viewer,"  using delayed reorthogonalization\n");
257:   }
258:   return(0);
259: }

261: SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
262: {
263:   EPS_ARNOLDI    *ctx;

267:   PetscNewLog(eps,&ctx);
268:   eps->data = (void*)ctx;

270:   eps->useds = PETSC_TRUE;

272:   eps->ops->solve          = EPSSolve_Arnoldi;
273:   eps->ops->setup          = EPSSetUp_Arnoldi;
274:   eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
275:   eps->ops->destroy        = EPSDestroy_Arnoldi;
276:   eps->ops->view           = EPSView_Arnoldi;
277:   eps->ops->backtransform  = EPSBackTransform_Default;
278:   eps->ops->computevectors = EPSComputeVectors_Schur;

280:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi);
281:   PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi);
282:   return(0);
283: }