Actual source code: narnoldi.c

slepc-3.17.2 2022-08-09
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:   PetscFunctionReturn(0);
 44: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

190:   NEPDeflationGetInvariantPair(extop,NULL,&H);
191:   MatGetSize(H,NULL,&ldh);
192:   DSSetType(nep->ds,DSNHEP);
193:   DSAllocate(nep->ds,PetscMax(nep->nconv,1));
194:   DSGetLeadingDimension(nep->ds,&ldds);
195:   MatDenseGetArrayRead(H,&Hp);
196:   DSGetArray(nep->ds,DS_MAT_A,&Ap);
197:   for (j=0;j<nep->nconv;j++)
198:     for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
199:   DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
200:   MatDenseRestoreArrayRead(H,&Hp);
201:   MatDestroy(&H);
202:   DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv);
203:   DSSolve(nep->ds,nep->eigr,nep->eigi);
204:   NEPDeflationReset(extop);
205:   VecDestroy(&u);
206:   VecDestroy(&r);
207:   BVDestroy(&Vext);
208:   DSDestroy(&ds);
209:   PetscFree(eigr);
210:   PetscFunctionReturn(0);
211: }

213: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
214: {
215:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

218:   ctx->lag = lag;
219:   PetscFunctionReturn(0);
220: }

222: /*@
223:    NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
224:    nonlinear solve.

226:    Logically Collective on nep

228:    Input Parameters:
229: +  nep - nonlinear eigenvalue solver
230: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
231:           computed within the nonlinear iteration, 2 means every second time
232:           the Jacobian is built, etc.

234:    Options Database Keys:
235: .  -nep_narnoldi_lag_preconditioner <lag> - the lag value

237:    Notes:
238:    The default is 1.
239:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

241:    Level: intermediate

243: .seealso: NEPNArnoldiGetLagPreconditioner()
244: @*/
245: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
246: {
249:   PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
250:   PetscFunctionReturn(0);
251: }

253: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
254: {
255:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

257:   *lag = ctx->lag;
258:   PetscFunctionReturn(0);
259: }

261: /*@
262:    NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

264:    Not Collective

266:    Input Parameter:
267: .  nep - nonlinear eigenvalue solver

269:    Output Parameter:
270: .  lag - the lag parameter

272:    Level: intermediate

274: .seealso: NEPNArnoldiSetLagPreconditioner()
275: @*/
276: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
277: {
280:   PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
281:   PetscFunctionReturn(0);
282: }

284: PetscErrorCode NEPSetFromOptions_NArnoldi(PetscOptionItems *PetscOptionsObject,NEP nep)
285: {
286:   PetscInt       i;
287:   PetscBool      flg;
288:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

290:   PetscOptionsHead(PetscOptionsObject,"NEP N-Arnoldi Options");
291:     i = 0;
292:     PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg);
293:     if (flg) NEPNArnoldiSetLagPreconditioner(nep,i);

295:   PetscOptionsTail();

297:   if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
298:   KSPSetFromOptions(ctx->ksp);
299:   PetscFunctionReturn(0);
300: }

302: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
303: {
304:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

306:   PetscObjectReference((PetscObject)ksp);
307:   KSPDestroy(&ctx->ksp);
308:   ctx->ksp = ksp;
309:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
310:   nep->state = NEP_STATE_INITIAL;
311:   PetscFunctionReturn(0);
312: }

314: /*@
315:    NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
316:    eigenvalue solver.

318:    Collective on nep

320:    Input Parameters:
321: +  nep - eigenvalue solver
322: -  ksp - the linear solver object

324:    Level: advanced

326: .seealso: NEPNArnoldiGetKSP()
327: @*/
328: PetscErrorCode NEPNArnoldiSetKSP(NEP nep,KSP ksp)
329: {
333:   PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
334:   PetscFunctionReturn(0);
335: }

337: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
338: {
339:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

341:   if (!ctx->ksp) {
342:     KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
343:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
344:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
345:     KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_");
346:     PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
347:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
348:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
349:     KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
350:   }
351:   *ksp = ctx->ksp;
352:   PetscFunctionReturn(0);
353: }

355: /*@
356:    NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
357:    the nonlinear eigenvalue solver.

359:    Not Collective

361:    Input Parameter:
362: .  nep - nonlinear eigenvalue solver

364:    Output Parameter:
365: .  ksp - the linear solver object

367:    Level: advanced

369: .seealso: NEPNArnoldiSetKSP()
370: @*/
371: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
372: {
375:   PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
376:   PetscFunctionReturn(0);
377: }

379: PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
380: {
381:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;
382:   PetscBool      isascii;

384:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
385:   if (isascii) {
386:     if (ctx->lag) PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag);
387:     if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
388:     PetscViewerASCIIPushTab(viewer);
389:     KSPView(ctx->ksp,viewer);
390:     PetscViewerASCIIPopTab(viewer);
391:   }
392:   PetscFunctionReturn(0);
393: }

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

399:   KSPReset(ctx->ksp);
400:   PetscFunctionReturn(0);
401: }

403: PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
404: {
405:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

407:   KSPDestroy(&ctx->ksp);
408:   PetscFree(nep->data);
409:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL);
410:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL);
411:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL);
412:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL);
413:   PetscFunctionReturn(0);
414: }

416: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
417: {
418:   NEP_NARNOLDI   *ctx;

420:   PetscNewLog(nep,&ctx);
421:   nep->data = (void*)ctx;
422:   ctx->lag  = 1;

424:   nep->useds = PETSC_TRUE;

426:   nep->ops->solve          = NEPSolve_NArnoldi;
427:   nep->ops->setup          = NEPSetUp_NArnoldi;
428:   nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
429:   nep->ops->reset          = NEPReset_NArnoldi;
430:   nep->ops->destroy        = NEPDestroy_NArnoldi;
431:   nep->ops->view           = NEPView_NArnoldi;
432:   nep->ops->computevectors = NEPComputeVectors_Schur;

434:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi);
435:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi);
436:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi);
437:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi);
438:   PetscFunctionReturn(0);
439: }