Actual source code: narnoldi.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 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>         /*I "slepcnep.h" I*/
 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: {
 36:   PetscBool      istrivial;

 39:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
 40:   if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
 41:   if (!nep->max_it) nep->max_it = nep->nev*nep->ncv;
 42:   if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");
 43:   if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NARNOLDI only available for split operator");

 45:   RGIsTrivial(nep->rg,&istrivial);
 46:   if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");

 48:   NEPAllocateSolution(nep,0);
 49:   NEPSetWorkVecs(nep,3);
 50:   return(0);
 51: }

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

 70:   /* get initial space and shift */
 71:   NEPGetDefaultShift(nep,&sigma);
 72:   if (!nep->nini) {
 73:     BVSetRandomColumn(nep->V,0);
 74:     BVNormColumn(nep->V,0,NORM_2,&nrm);
 75:     BVScaleColumn(nep->V,0,1.0/nrm);
 76:     n = 1;
 77:   } else n = nep->nini;

 79:   if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
 80:   NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
 81:   NEPDeflationCreateBV(extop,nep->ncv,&Vext);

 83:   /* prepare linear solver */
 84:   NEPDeflationSolveSetUp(extop,sigma);

 86:   BVGetColumn(Vext,0,&f);
 87:   VecDuplicate(f,&r);
 88:   VecDuplicate(f,&u);
 89:   BVGetColumn(nep->V,0,&uu);
 90:   NEPDeflationCopyToExtendedVec(extop,uu,NULL,f,PETSC_FALSE);
 91:   BVRestoreColumn(nep->V,0,&uu);
 92:   VecCopy(f,r);
 93:   NEPDeflationFunctionSolve(extop,r,f);
 94:   VecNorm(f,NORM_2,&nrm);
 95:   VecScale(f,1.0/nrm);
 96:   BVRestoreColumn(Vext,0,&f);

 98:   DSCreate(PetscObjectComm((PetscObject)nep),&ds);
 99:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ds);
100:   DSSetType(ds,DSNEP);
101:   DSNEPSetFN(ds,nep->nt,nep->f);
102:   DSAllocate(ds,nep->ncv);
103:   DSGetSlepcSC(ds,&sc);
104:   sc->comparison    = nep->sc->comparison;
105:   sc->comparisonctx = nep->sc->comparisonctx;
106:   DSSetFromOptions(ds);

108:   /* build projected matrices for initial space */
109:   DSSetDimensions(ds,n,0,0,0);
110:   NEPDeflationProjectOperator(extop,Vext,ds,0,n);

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

114:   /* Restart loop */
115:   while (nep->reason == NEP_CONVERGED_ITERATING) {
116:     nep->its++;

118:     /* solve projected problem */
119:     DSSetDimensions(ds,n,0,0,0);
120:     DSSetState(ds,DS_STATE_RAW);
121:     DSSolve(ds,eigr,NULL);
122:     DSSynchronize(ds,eigr,NULL);
123:     if (nep->its>1) lambda2 = lambda;
124:     lambda = eigr[0];
125:     nep->eigr[nep->nconv] = lambda;

127:     /* compute Ritz vector, x = V*s */
128:     DSGetArray(ds,DS_MAT_X,&X);
129:     BVSetActiveColumns(Vext,0,n);
130:     BVMultVec(Vext,1.0,0.0,u,X);
131:     DSRestoreArray(ds,DS_MAT_X,&X);

133:     /* compute the residual, r = T(lambda)*x */
134:     NEPDeflationComputeFunction(extop,lambda,&T);
135:     MatMult(T,u,r);

137:     /* convergence test */
138:     VecNorm(r,NORM_2,&resnorm);
139:     if (nep->its>1) perr=nep->errest[nep->nconv];
140:     (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
141:     if (nep->errest[nep->nconv]<=nep->tol) {
142:       nep->nconv = nep->nconv + 1;
143:       NEPDeflationLocking(extop,u,lambda);
144:       skip = PETSC_TRUE;
145:     }
146:     (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
147:     if (!skip || nep->reason>0) {
148:       NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
149:     }

151:     if (nep->reason == NEP_CONVERGED_ITERATING) {
152:       if (!skip) {
153:         if (n>=nep->ncv) {
154:           nep->reason = NEP_DIVERGED_SUBSPACE_EXHAUSTED;
155:           break;
156:         }
157:         if (ctx->lag && !(nep->its%ctx->lag) && nep->its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) {
158:           NEPDeflationSolveSetUp(extop,lambda2);
159:         }

161:         /* continuation vector: f = T(sigma)\r */
162:         BVGetColumn(Vext,n,&f);
163:         NEPDeflationFunctionSolve(extop,r,f);
164:         BVRestoreColumn(Vext,n,&f);
165:         KSPGetConvergedReason(ctx->ksp,&kspreason);
166:         if (kspreason<0) {
167:           PetscInfo1(nep,"iter=%D, linear solve failed, stopping solve\n",nep->its);
168:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
169:           break;
170:         }

172:         /* orthonormalize */
173:         BVOrthonormalizeColumn(Vext,n,PETSC_FALSE,&beta,&breakdown);
174:         if (breakdown || beta==0.0) {
175:           PetscInfo1(nep,"iter=%D, orthogonalization failed, stopping solve\n",nep->its);
176:           nep->reason = NEP_DIVERGED_BREAKDOWN;
177:           break;
178:         }

180:         /* update projected matrices */
181:         DSSetDimensions(ds,n+1,0,0,0);
182:         NEPDeflationProjectOperator(extop,Vext,ds,n,n+1);
183:         n++;
184:       } else {
185:         nep->its--;  /* do not count this as a full iteration */
186:         BVGetColumn(Vext,0,&f);
187:         NEPDeflationSetRandomVec(extop,f);
188:         NEPDeflationSolveSetUp(extop,sigma);
189:         VecCopy(f,r);
190:         NEPDeflationFunctionSolve(extop,r,f);
191:         VecNorm(f,NORM_2,&nrm);
192:         VecScale(f,1.0/nrm);
193:         BVRestoreColumn(Vext,0,&f);
194:         n = 1;
195:         DSSetDimensions(ds,n,0,0,0);
196:         NEPDeflationProjectOperator(extop,Vext,ds,n-1,n);
197:         skip = PETSC_FALSE;
198:       }
199:     }
200:   }

202:   NEPDeflationGetInvariantPair(extop,NULL,&H);
203:   MatGetSize(H,NULL,&ldh);
204:   DSSetType(nep->ds,DSNHEP);
205:   DSAllocate(nep->ds,PetscMax(nep->nconv,1));
206:   DSGetLeadingDimension(nep->ds,&ldds);
207:   MatDenseGetArray(H,&Hp);
208:   DSGetArray(nep->ds,DS_MAT_A,&Ap);
209:   for (j=0;j<nep->nconv;j++)
210:     for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
211:   DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
212:   MatDenseRestoreArray(H,&Hp);
213:   MatDestroy(&H);
214:   DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
215:   DSSolve(nep->ds,nep->eigr,nep->eigi);
216:   NEPDeflationReset(extop);
217:   VecDestroy(&u);
218:   VecDestroy(&r);
219:   BVDestroy(&Vext);
220:   DSDestroy(&ds);
221:   PetscFree(eigr);
222:   return(0);
223: }

225: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
226: {
227:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

230:   if (lag<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
231:   ctx->lag = lag;
232:   return(0);
233: }

235: /*@
236:    NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
237:    nonlinear solve.

239:    Logically Collective on NEP

241:    Input Parameters:
242: +  nep - nonlinear eigenvalue solver
243: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
244:           computed within the nonlinear iteration, 2 means every second time
245:           the Jacobian is built, etc.

247:    Options Database Keys:
248: .  -nep_narnoldi_lag_preconditioner <lag>

250:    Notes:
251:    The default is 1.
252:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

254:    Level: intermediate

256: .seealso: NEPNArnoldiGetLagPreconditioner()
257: @*/
258: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
259: {

265:   PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
266:   return(0);
267: }

269: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
270: {
271:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

274:   *lag = ctx->lag;
275:   return(0);
276: }

278: /*@
279:    NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

281:    Not Collective

283:    Input Parameter:
284: .  nep - nonlinear eigenvalue solver

286:    Output Parameter:
287: .  lag - the lag parameter

289:    Level: intermediate

291: .seealso: NEPNArnoldiSetLagPreconditioner()
292: @*/
293: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
294: {

300:   PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
301:   return(0);
302: }

304: PetscErrorCode NEPSetFromOptions_NArnoldi(PetscOptionItems *PetscOptionsObject,NEP nep)
305: {
307:   PetscInt       i;
308:   PetscBool      flg;
309:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

312:   PetscOptionsHead(PetscOptionsObject,"NEP N-Arnoldi Options");
313:     i = 0;
314:     PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg);
315:     if (flg) { NEPNArnoldiSetLagPreconditioner(nep,i); }

317:   PetscOptionsTail();

319:   if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
320:   KSPSetFromOptions(ctx->ksp);
321:   return(0);
322: }

324: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
325: {
327:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

330:   PetscObjectReference((PetscObject)ksp);
331:   KSPDestroy(&ctx->ksp);
332:   ctx->ksp = ksp;
333:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
334:   nep->state = NEP_STATE_INITIAL;
335:   return(0);
336: }

338: /*@
339:    NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
340:    eigenvalue solver.

342:    Collective on NEP

344:    Input Parameters:
345: +  nep - eigenvalue solver
346: -  ksp - the linear solver object

348:    Level: advanced

350: .seealso: NEPNArnoldiGetKSP()
351: @*/
352: PetscErrorCode NEPNArnoldiSetKSP(NEP nep,KSP ksp)
353: {

360:   PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
361:   return(0);
362: }

364: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
365: {
367:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

370:   if (!ctx->ksp) {
371:     KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
372:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
373:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
374:     KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_");
375:     PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
376:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
377:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
378:     KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
379:   }
380:   *ksp = ctx->ksp;
381:   return(0);
382: }

384: /*@
385:    NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
386:    the nonlinear eigenvalue solver.

388:    Not Collective

390:    Input Parameter:
391: .  nep - nonlinear eigenvalue solver

393:    Output Parameter:
394: .  ksp - the linear solver object

396:    Level: advanced

398: .seealso: NEPNArnoldiSetKSP()
399: @*/
400: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
401: {

407:   PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
408:   return(0);
409: }

411: PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
412: {
414:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;
415:   PetscBool      isascii;

418:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
419:   if (isascii) {
420:     if (ctx->lag) {
421:       PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %D iterations\n",ctx->lag);
422:     }
423:     if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
424:     PetscViewerASCIIPushTab(viewer);
425:     KSPView(ctx->ksp,viewer);
426:     PetscViewerASCIIPopTab(viewer);
427:   }
428:   return(0);
429: }

431: PetscErrorCode NEPReset_NArnoldi(NEP nep)
432: {
434:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

437:   KSPReset(ctx->ksp);
438:   return(0);
439: }

441: PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
442: {
444:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

447:   KSPDestroy(&ctx->ksp);
448:   PetscFree(nep->data);
449:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL);
450:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL);
451:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL);
452:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL);
453:   return(0);
454: }

456: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
457: {
459:   NEP_NARNOLDI   *ctx;

462:   PetscNewLog(nep,&ctx);
463:   nep->data = (void*)ctx;
464:   ctx->lag  = 1;

466:   nep->useds = PETSC_TRUE;

468:   nep->ops->solve          = NEPSolve_NArnoldi;
469:   nep->ops->setup          = NEPSetUp_NArnoldi;
470:   nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
471:   nep->ops->reset          = NEPReset_NArnoldi;
472:   nep->ops->destroy        = NEPDestroy_NArnoldi;
473:   nep->ops->view           = NEPView_NArnoldi;
474:   nep->ops->computevectors = NEPComputeVectors_Schur;

476:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi);
477:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi);
478:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi);
479:   PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi);
480:   return(0);
481: }