Actual source code: narnoldi.c

slepc-3.18.0 2022-10-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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 nonlinear eigensolver: "narnoldi"

 13:    Method: Nonlinear Arnoldi

 15:    Algorithm:

 17:        Arnoldi for nonlinear eigenproblems.

 19:    References:

 21:        [1] H. Voss, "An Arnoldi method for nonlinear eigenvalue problems",
 22:            BIT 44:387-401, 2004.
 23: */

 25: #include <slepc/private/nepimpl.h>
 26: #include <../src/nep/impls/nepdefl.h>

 28: typedef struct {
 29:   PetscInt lag;             /* interval to rebuild preconditioner */
 30:   KSP      ksp;             /* linear solver object */
 31: } NEP_NARNOLDI;

 33: PetscErrorCode NEPSetUp_NArnoldi(NEP nep)
 34: {
 35:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
 37:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = nep->nev*nep->ncv;
 38:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
 40:   NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
 41:   NEPAllocateSolution(nep,0);
 42:   NEPSetWorkVecs(nep,3);
 43:   return 0;
 44: }

 46: PetscErrorCode NEPSolve_NArnoldi(NEP nep)
 47: {
 48:   NEP_NARNOLDI       *ctx = (NEP_NARNOLDI*)nep->data;
 49:   Mat                T,H,A;
 50:   Vec                f,r,u,uu;
 51:   PetscScalar        *X,lambda=0.0,lambda2=0.0,*eigr,sigma;
 52:   PetscReal          beta,resnorm=0.0,nrm,perr=0.0;
 53:   PetscInt           n;
 54:   PetscBool          breakdown,skip=PETSC_FALSE;
 55:   BV                 Vext;
 56:   DS                 ds;
 57:   NEP_EXT_OP         extop=NULL;
 58:   SlepcSC            sc;
 59:   KSPConvergedReason kspreason;

 61:   /* get initial space and shift */
 62:   NEPGetDefaultShift(nep,&sigma);
 63:   if (!nep->nini) {
 64:     BVSetRandomColumn(nep->V,0);
 65:     BVNormColumn(nep->V,0,NORM_2,&nrm);
 66:     BVScaleColumn(nep->V,0,1.0/nrm);
 67:     n = 1;
 68:   } else n = nep->nini;

 70:   if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
 71:   NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
 72:   NEPDeflationCreateBV(extop,nep->ncv,&Vext);

 74:   /* prepare linear solver */
 75:   NEPDeflationSolveSetUp(extop,sigma);

 77:   BVGetColumn(Vext,0,&f);
 78:   VecDuplicate(f,&r);
 79:   VecDuplicate(f,&u);
 80:   BVGetColumn(nep->V,0,&uu);
 81:   NEPDeflationCopyToExtendedVec(extop,uu,NULL,f,PETSC_FALSE);
 82:   BVRestoreColumn(nep->V,0,&uu);
 83:   VecCopy(f,r);
 84:   NEPDeflationFunctionSolve(extop,r,f);
 85:   VecNorm(f,NORM_2,&nrm);
 86:   VecScale(f,1.0/nrm);
 87:   BVRestoreColumn(Vext,0,&f);

 89:   DSCreate(PetscObjectComm((PetscObject)nep),&ds);
 90:   DSSetType(ds,DSNEP);
 91:   DSNEPSetFN(ds,nep->nt,nep->f);
 92:   DSAllocate(ds,nep->ncv);
 93:   DSGetSlepcSC(ds,&sc);
 94:   sc->comparison    = nep->sc->comparison;
 95:   sc->comparisonctx = nep->sc->comparisonctx;
 96:   DSSetFromOptions(ds);

 98:   /* build projected matrices for initial space */
 99:   DSSetDimensions(ds,n,0,0);
100:   NEPDeflationProjectOperator(extop,Vext,ds,0,n);

102:   PetscMalloc1(nep->ncv,&eigr);

104:   /* Restart loop */
105:   while (nep->reason == NEP_CONVERGED_ITERATING) {
106:     nep->its++;

108:     /* solve projected problem */
109:     DSSetDimensions(ds,n,0,0);
110:     DSSetState(ds,DS_STATE_RAW);
111:     DSSolve(ds,eigr,NULL);
112:     DSSynchronize(ds,eigr,NULL);
113:     if (nep->its>1) lambda2 = lambda;
114:     lambda = eigr[0];
115:     nep->eigr[nep->nconv] = lambda;

117:     /* compute Ritz vector, x = V*s */
118:     DSGetArray(ds,DS_MAT_X,&X);
119:     BVSetActiveColumns(Vext,0,n);
120:     BVMultVec(Vext,1.0,0.0,u,X);
121:     DSRestoreArray(ds,DS_MAT_X,&X);

123:     /* compute the residual, r = T(lambda)*x */
124:     NEPDeflationComputeFunction(extop,lambda,&T);
125:     MatMult(T,u,r);

127:     /* convergence test */
128:     VecNorm(r,NORM_2,&resnorm);
129:     if (nep->its>1) perr=nep->errest[nep->nconv];
130:     (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
131:     if (nep->errest[nep->nconv]<=nep->tol) {
132:       nep->nconv = nep->nconv + 1;
133:       NEPDeflationLocking(extop,u,lambda);
134:       skip = PETSC_TRUE;
135:     }
136:     (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
137:     if (!skip || nep->reason>0) NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);

139:     if (nep->reason == NEP_CONVERGED_ITERATING) {
140:       if (!skip) {
141:         if (n>=nep->ncv) {
142:           nep->reason = NEP_DIVERGED_SUBSPACE_EXHAUSTED;
143:           break;
144:         }
145:         if (ctx->lag && !(nep->its%ctx->lag) && nep->its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) NEPDeflationSolveSetUp(extop,lambda2);

147:         /* continuation vector: f = T(sigma)\r */
148:         BVGetColumn(Vext,n,&f);
149:         NEPDeflationFunctionSolve(extop,r,f);
150:         BVRestoreColumn(Vext,n,&f);
151:         KSPGetConvergedReason(ctx->ksp,&kspreason);
152:         if (kspreason<0) {
153:           PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its);
154:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
155:           break;
156:         }

158:         /* orthonormalize */
159:         BVOrthonormalizeColumn(Vext,n,PETSC_FALSE,&beta,&breakdown);
160:         if (breakdown || beta==0.0) {
161:           PetscInfo(nep,"iter=%" PetscInt_FMT ", orthogonalization failed, stopping solve\n",nep->its);
162:           nep->reason = NEP_DIVERGED_BREAKDOWN;
163:           break;
164:         }

166:         /* update projected matrices */
167:         DSSetDimensions(ds,n+1,0,0);
168:         NEPDeflationProjectOperator(extop,Vext,ds,n,n+1);
169:         n++;
170:       } else {
171:         nep->its--;  /* do not count this as a full iteration */
172:         BVGetColumn(Vext,0,&f);
173:         NEPDeflationSetRandomVec(extop,f);
174:         NEPDeflationSolveSetUp(extop,sigma);
175:         VecCopy(f,r);
176:         NEPDeflationFunctionSolve(extop,r,f);
177:         VecNorm(f,NORM_2,&nrm);
178:         VecScale(f,1.0/nrm);
179:         BVRestoreColumn(Vext,0,&f);
180:         n = 1;
181:         DSSetDimensions(ds,n,0,0);
182:         NEPDeflationProjectOperator(extop,Vext,ds,n-1,n);
183:         skip = PETSC_FALSE;
184:       }
185:     }
186:   }

188:   NEPDeflationGetInvariantPair(extop,NULL,&H);
189:   DSSetType(nep->ds,DSNHEP);
190:   DSAllocate(nep->ds,PetscMax(nep->nconv,1));
191:   DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv);
192:   DSGetMat(nep->ds,DS_MAT_A,&A);
193:   MatCopy(H,A,SAME_NONZERO_PATTERN);
194:   DSRestoreMat(nep->ds,DS_MAT_A,&A);
195:   MatDestroy(&H);
196:   DSSolve(nep->ds,nep->eigr,nep->eigi);
197:   NEPDeflationReset(extop);
198:   VecDestroy(&u);
199:   VecDestroy(&r);
200:   BVDestroy(&Vext);
201:   DSDestroy(&ds);
202:   PetscFree(eigr);
203:   return 0;
204: }

206: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
207: {
208:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

211:   ctx->lag = lag;
212:   return 0;
213: }

215: /*@
216:    NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
217:    nonlinear solve.

219:    Logically Collective on nep

221:    Input Parameters:
222: +  nep - nonlinear eigenvalue solver
223: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
224:           computed within the nonlinear iteration, 2 means every second time
225:           the Jacobian is built, etc.

227:    Options Database Keys:
228: .  -nep_narnoldi_lag_preconditioner <lag> - the lag value

230:    Notes:
231:    The default is 1.
232:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

234:    Level: intermediate

236: .seealso: NEPNArnoldiGetLagPreconditioner()
237: @*/
238: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
239: {
242:   PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
243:   return 0;
244: }

246: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
247: {
248:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

250:   *lag = ctx->lag;
251:   return 0;
252: }

254: /*@
255:    NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

257:    Not Collective

259:    Input Parameter:
260: .  nep - nonlinear eigenvalue solver

262:    Output Parameter:
263: .  lag - the lag parameter

265:    Level: intermediate

267: .seealso: NEPNArnoldiSetLagPreconditioner()
268: @*/
269: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
270: {
273:   PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
274:   return 0;
275: }

277: PetscErrorCode NEPSetFromOptions_NArnoldi(NEP nep,PetscOptionItems *PetscOptionsObject)
278: {
279:   PetscInt       i;
280:   PetscBool      flg;
281:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

283:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP N-Arnoldi Options");
284:     i = 0;
285:     PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg);
286:     if (flg) NEPNArnoldiSetLagPreconditioner(nep,i);

288:   PetscOptionsHeadEnd();

290:   if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
291:   KSPSetFromOptions(ctx->ksp);
292:   return 0;
293: }

295: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
296: {
297:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

299:   PetscObjectReference((PetscObject)ksp);
300:   KSPDestroy(&ctx->ksp);
301:   ctx->ksp = ksp;
302:   nep->state = NEP_STATE_INITIAL;
303:   return 0;
304: }

306: /*@
307:    NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
308:    eigenvalue solver.

310:    Collective on nep

312:    Input Parameters:
313: +  nep - eigenvalue solver
314: -  ksp - the linear solver object

316:    Level: advanced

318: .seealso: NEPNArnoldiGetKSP()
319: @*/
320: PetscErrorCode NEPNArnoldiSetKSP(NEP nep,KSP ksp)
321: {
325:   PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
326:   return 0;
327: }

329: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
330: {
331:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

333:   if (!ctx->ksp) {
334:     KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
335:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
336:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
337:     KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_");
338:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
339:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
340:     KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
341:   }
342:   *ksp = ctx->ksp;
343:   return 0;
344: }

346: /*@
347:    NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
348:    the nonlinear eigenvalue solver.

350:    Not Collective

352:    Input Parameter:
353: .  nep - nonlinear eigenvalue solver

355:    Output Parameter:
356: .  ksp - the linear solver object

358:    Level: advanced

360: .seealso: NEPNArnoldiSetKSP()
361: @*/
362: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
363: {
366:   PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
367:   return 0;
368: }

370: PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
371: {
372:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;
373:   PetscBool      isascii;

375:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
376:   if (isascii) {
377:     if (ctx->lag) PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag);
378:     if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
379:     PetscViewerASCIIPushTab(viewer);
380:     KSPView(ctx->ksp,viewer);
381:     PetscViewerASCIIPopTab(viewer);
382:   }
383:   return 0;
384: }

386: PetscErrorCode NEPReset_NArnoldi(NEP nep)
387: {
388:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

390:   KSPReset(ctx->ksp);
391:   return 0;
392: }

394: PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
395: {
396:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

398:   KSPDestroy(&ctx->ksp);
399:   PetscFree(nep->data);
400:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL);
401:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL);
402:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL);
403:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL);
404:   return 0;
405: }

407: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
408: {
409:   NEP_NARNOLDI   *ctx;

411:   PetscNew(&ctx);
412:   nep->data = (void*)ctx;
413:   ctx->lag  = 1;

415:   nep->useds = PETSC_TRUE;

417:   nep->ops->solve          = NEPSolve_NArnoldi;
418:   nep->ops->setup          = NEPSetUp_NArnoldi;
419:   nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
420:   nep->ops->reset          = NEPReset_NArnoldi;
421:   nep->ops->destroy        = NEPDestroy_NArnoldi;
422:   nep->ops->view           = NEPView_NArnoldi;
423:   nep->ops->computevectors = NEPComputeVectors_Schur;

425:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi);
426:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi);
427:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi);
428:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi);
429:   return 0;
430: }