Actual source code: rii.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: "rii"

 13:    Method: Residual inverse iteration

 15:    Algorithm:

 17:        Simple residual inverse iteration with varying shift.

 19:    References:

 21:        [1] A. Neumaier, "Residual inverse iteration for the nonlinear
 22:            eigenvalue problem", SIAM J. Numer. Anal. 22(5):914-923, 1985.
 23: */

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

 28: typedef struct {
 29:   PetscInt  max_inner_it;     /* maximum number of Newton iterations */
 30:   PetscInt  lag;              /* interval to rebuild preconditioner */
 31:   PetscBool cctol;            /* constant correction tolerance */
 32:   KSP       ksp;              /* linear solver object */
 33: } NEP_RII;

 35: PetscErrorCode NEPSetUp_RII(NEP nep)
 36: {
 38:   PetscBool      istrivial;

 41:   if (nep->ncv) { PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"); }
 42:   nep->ncv = nep->nev;
 43:   if (nep->mpd) { PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"); }
 44:   nep->mpd = nep->nev;
 45:   if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
 46:   if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");

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

 51:   NEPAllocateSolution(nep,0);
 52:   NEPSetWorkVecs(nep,2);
 53:   return(0);
 54: }

 56: PetscErrorCode NEPSolve_RII(NEP nep)
 57: {
 58:   PetscErrorCode     ierr;
 59:   NEP_RII            *ctx = (NEP_RII*)nep->data;
 60:   Mat                T,Tp,H;
 61:   Vec                uu,u,r,delta;
 62:   PetscScalar        lambda,lambda2,sigma,a1,a2,corr,*Hp,*Ap;
 63:   PetscReal          nrm,resnorm=1.0,ktol=0.1,perr;
 64:   PetscBool          skip=PETSC_FALSE;
 65:   PetscInt           inner_its,its=0,ldh,ldds,i,j;
 66:   NEP_EXT_OP         extop=NULL;
 67:   KSPConvergedReason kspreason;

 70:   /* get initial approximation of eigenvalue and eigenvector */
 71:   NEPGetDefaultShift(nep,&sigma);
 72:   lambda = sigma;
 73:   if (!nep->nini) {
 74:     BVSetRandomColumn(nep->V,0);
 75:     BVNormColumn(nep->V,0,NORM_2,&nrm);
 76:     BVScaleColumn(nep->V,0,1.0/nrm);
 77:   }
 78:   if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
 79:   NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
 80:   NEPDeflationCreateVec(extop,&u);
 81:   VecDuplicate(u,&r);
 82:   VecDuplicate(u,&delta);
 83:   BVGetColumn(nep->V,0,&uu);
 84:   NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE);
 85:   BVRestoreColumn(nep->V,0,&uu);

 87:   /* prepare linear solver */
 88:   NEPDeflationSolveSetUp(extop,sigma);

 90:   VecCopy(u,r);
 91:   NEPDeflationFunctionSolve(extop,r,u);
 92:   VecNorm(u,NORM_2,&nrm);
 93:   VecScale(u,1.0/nrm);

 95:   /* Restart loop */
 96:   while (nep->reason == NEP_CONVERGED_ITERATING) {
 97:     its++;

 99:     /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
100:        estimate as starting value. */
101:     inner_its=0;
102:     lambda2 = lambda;
103:     do {
104:       NEPDeflationComputeFunction(extop,lambda,&T);
105:       MatMult(T,u,r);
106:       VecDot(r,u,&a1);
107:       if (inner_its && PetscAbsScalar(a1)/PetscAbsScalar(lambda)<=PETSC_SQRT_MACHINE_EPSILON) break;
108:       NEPDeflationComputeJacobian(extop,lambda,&Tp);
109:       MatMult(Tp,u,r);
110:       VecDot(r,u,&a2);
111:       corr = a1/a2;
112:       lambda = lambda - corr;
113:       inner_its++;
114:     } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);

116:     /* form residual,  r = T(lambda)*u */
117:     NEPDeflationComputeFunction(extop,lambda,&T);
118:     MatMult(T,u,r);

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

136:     if (nep->reason == NEP_CONVERGED_ITERATING) {
137:       if (!skip) {
138:         /* update preconditioner and set adaptive tolerance */
139:         if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) {
140:           NEPDeflationSolveSetUp(extop,lambda2);
141:         }
142:         if (!ctx->cctol) {
143:           ktol = PetscMax(ktol/2.0,PETSC_MACHINE_EPSILON*10.0);
144:           KSPSetTolerances(ctx->ksp,ktol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
145:         }

147:         /* eigenvector correction: delta = T(sigma)\r */
148:         NEPDeflationFunctionSolve(extop,r,delta);
149:         KSPGetConvergedReason(ctx->ksp,&kspreason);
150:         if (kspreason<0) {
151:           PetscInfo1(nep,"iter=%D, linear solve failed, stopping solve\n",nep->its);
152:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
153:           break;
154:         }

156:         /* update eigenvector: u = u - delta */
157:         VecAXPY(u,-1.0,delta);

159:         /* normalize eigenvector */
160:         VecNormalize(u,NULL);
161:       } else {
162:         its = -1;
163:         NEPDeflationSetRandomVec(extop,u);
164:         NEPDeflationSolveSetUp(extop,sigma);
165:         VecCopy(u,r);
166:         NEPDeflationFunctionSolve(extop,r,u);
167:         VecNorm(u,NORM_2,&nrm);
168:         VecScale(u,1.0/nrm);
169:         lambda = sigma;
170:         skip = PETSC_FALSE;
171:       }
172:     }
173:   }
174:   NEPDeflationGetInvariantPair(extop,NULL,&H);
175:   MatGetSize(H,NULL,&ldh);
176:   DSSetType(nep->ds,DSNHEP);
177:   DSAllocate(nep->ds,PetscMax(nep->nconv,1));
178:   DSGetLeadingDimension(nep->ds,&ldds);
179:   MatDenseGetArray(H,&Hp);
180:   DSGetArray(nep->ds,DS_MAT_A,&Ap);
181:   for (j=0;j<nep->nconv;j++)
182:     for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
183:   DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
184:   MatDenseRestoreArray(H,&Hp);
185:   MatDestroy(&H);
186:   DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
187:   DSSolve(nep->ds,nep->eigr,nep->eigi);
188:   NEPDeflationReset(extop);
189:   VecDestroy(&u);
190:   VecDestroy(&r);
191:   VecDestroy(&delta);
192:   return(0);
193: }

195: PetscErrorCode NEPSetFromOptions_RII(PetscOptionItems *PetscOptionsObject,NEP nep)
196: {
198:   NEP_RII        *ctx = (NEP_RII*)nep->data;
199:   PetscBool      flg;
200:   PetscInt       i;

203:   PetscOptionsHead(PetscOptionsObject,"NEP RII Options");

205:     PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&ctx->max_inner_it,NULL);

207:     PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL);

209:     i = 0;
210:     PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg);
211:     if (flg) { NEPRIISetLagPreconditioner(nep,i); }

213:   PetscOptionsTail();

215:   if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
216:   KSPSetFromOptions(ctx->ksp);
217:   return(0);
218: }

220: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
221: {
222:   NEP_RII *ctx = (NEP_RII*)nep->data;

225:   ctx->max_inner_it = its;
226:   return(0);
227: }

229: /*@
230:    NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
231:    used in the RII solver. These are the Newton iterations related to the computation
232:    of the nonlinear Rayleigh functional.

234:    Logically Collective on NEP

236:    Input Parameters:
237: +  nep - nonlinear eigenvalue solver
238: -  its - maximum inner iterations

240:    Level: advanced

242: .seealso: NEPRIIGetMaximumIterations()
243: @*/
244: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
245: {

251:   PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
252:   return(0);
253: }

255: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
256: {
257:   NEP_RII *ctx = (NEP_RII*)nep->data;

260:   *its = ctx->max_inner_it;
261:   return(0);
262: }

264: /*@
265:    NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.

267:    Not Collective

269:    Input Parameter:
270: .  nep - nonlinear eigenvalue solver

272:    Output Parameter:
273: .  its - maximum inner iterations

275:    Level: advanced

277: .seealso: NEPRIISetMaximumIterations()
278: @*/
279: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
280: {

286:   PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
287:   return(0);
288: }

290: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
291: {
292:   NEP_RII *ctx = (NEP_RII*)nep->data;

295:   if (lag<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
296:   ctx->lag = lag;
297:   return(0);
298: }

300: /*@
301:    NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
302:    nonlinear solve.

304:    Logically Collective on NEP

306:    Input Parameters:
307: +  nep - nonlinear eigenvalue solver
308: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
309:           computed within the nonlinear iteration, 2 means every second time
310:           the Jacobian is built, etc.

312:    Options Database Keys:
313: .  -nep_rii_lag_preconditioner <lag>

315:    Notes:
316:    The default is 1.
317:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

319:    Level: intermediate

321: .seealso: NEPRIIGetLagPreconditioner()
322: @*/
323: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
324: {

330:   PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
331:   return(0);
332: }

334: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
335: {
336:   NEP_RII *ctx = (NEP_RII*)nep->data;

339:   *lag = ctx->lag;
340:   return(0);
341: }

343: /*@
344:    NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

346:    Not Collective

348:    Input Parameter:
349: .  nep - nonlinear eigenvalue solver

351:    Output Parameter:
352: .  lag - the lag parameter

354:    Level: intermediate

356: .seealso: NEPRIISetLagPreconditioner()
357: @*/
358: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
359: {

365:   PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
366:   return(0);
367: }

369: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
370: {
371:   NEP_RII *ctx = (NEP_RII*)nep->data;

374:   ctx->cctol = cct;
375:   return(0);
376: }

378: /*@
379:    NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
380:    in the linear solver constant.

382:    Logically Collective on NEP

384:    Input Parameters:
385: +  nep - nonlinear eigenvalue solver
386: -  cct - a boolean value

388:    Options Database Keys:
389: .  -nep_rii_const_correction_tol <bool> - set the boolean flag

391:    Notes:
392:    By default, an exponentially decreasing tolerance is set in the KSP used
393:    within the nonlinear iteration, so that each Newton iteration requests
394:    better accuracy than the previous one. The constant correction tolerance
395:    flag stops this behaviour.

397:    Level: intermediate

399: .seealso: NEPRIIGetConstCorrectionTol()
400: @*/
401: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
402: {

408:   PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
409:   return(0);
410: }

412: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
413: {
414:   NEP_RII *ctx = (NEP_RII*)nep->data;

417:   *cct = ctx->cctol;
418:   return(0);
419: }

421: /*@
422:    NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.

424:    Not Collective

426:    Input Parameter:
427: .  nep - nonlinear eigenvalue solver

429:    Output Parameter:
430: .  cct - the value of the constant tolerance flag

432:    Level: intermediate

434: .seealso: NEPRIISetConstCorrectionTol()
435: @*/
436: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
437: {

443:   PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
444:   return(0);
445: }

447: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
448: {
450:   NEP_RII        *ctx = (NEP_RII*)nep->data;

453:   PetscObjectReference((PetscObject)ksp);
454:   KSPDestroy(&ctx->ksp);
455:   ctx->ksp = ksp;
456:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
457:   nep->state = NEP_STATE_INITIAL;
458:   return(0);
459: }

461: /*@
462:    NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
463:    eigenvalue solver.

465:    Collective on NEP

467:    Input Parameters:
468: +  nep - eigenvalue solver
469: -  ksp - the linear solver object

471:    Level: advanced

473: .seealso: NEPRIIGetKSP()
474: @*/
475: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
476: {

483:   PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
484:   return(0);
485: }

487: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
488: {
490:   NEP_RII        *ctx = (NEP_RII*)nep->data;

493:   if (!ctx->ksp) {
494:     KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
495:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
496:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
497:     KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_");
498:     PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
499:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
500:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
501:     KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
502:   }
503:   *ksp = ctx->ksp;
504:   return(0);
505: }

507: /*@
508:    NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
509:    the nonlinear eigenvalue solver.

511:    Not Collective

513:    Input Parameter:
514: .  nep - nonlinear eigenvalue solver

516:    Output Parameter:
517: .  ksp - the linear solver object

519:    Level: advanced

521: .seealso: NEPRIISetKSP()
522: @*/
523: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
524: {

530:   PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
531:   return(0);
532: }

534: PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
535: {
537:   NEP_RII        *ctx = (NEP_RII*)nep->data;
538:   PetscBool      isascii;

541:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
542:   if (isascii) {
543:     if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
544:     PetscViewerASCIIPrintf(viewer,"  maximum number of inner iterations: %D\n",ctx->max_inner_it);
545:     if (ctx->cctol) {
546:       PetscViewerASCIIPrintf(viewer,"  using a constant tolerance for the linear solver\n");
547:     }
548:     if (ctx->lag) {
549:       PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %D iterations\n",ctx->lag);
550:     }
551:     PetscViewerASCIIPushTab(viewer);
552:     KSPView(ctx->ksp,viewer);
553:     PetscViewerASCIIPopTab(viewer);
554:   }
555:   return(0);
556: }

558: PetscErrorCode NEPReset_RII(NEP nep)
559: {
561:   NEP_RII        *ctx = (NEP_RII*)nep->data;

564:   KSPReset(ctx->ksp);
565:   return(0);
566: }

568: PetscErrorCode NEPDestroy_RII(NEP nep)
569: {
571:   NEP_RII        *ctx = (NEP_RII*)nep->data;

574:   KSPDestroy(&ctx->ksp);
575:   PetscFree(nep->data);
576:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL);
577:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL);
578:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL);
579:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL);
580:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL);
581:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL);
582:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL);
583:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL);
584:   return(0);
585: }

587: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
588: {
590:   NEP_RII        *ctx;

593:   PetscNewLog(nep,&ctx);
594:   nep->data = (void*)ctx;
595:   ctx->max_inner_it = 10;
596:   ctx->lag          = 1;
597:   ctx->cctol        = PETSC_FALSE;

599:   nep->useds = PETSC_TRUE;

601:   nep->ops->solve          = NEPSolve_RII;
602:   nep->ops->setup          = NEPSetUp_RII;
603:   nep->ops->setfromoptions = NEPSetFromOptions_RII;
604:   nep->ops->reset          = NEPReset_RII;
605:   nep->ops->destroy        = NEPDestroy_RII;
606:   nep->ops->view           = NEPView_RII;
607:   nep->ops->computevectors = NEPComputeVectors_Schur;

609:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII);
610:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII);
611:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII);
612:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII);
613:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII);
614:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII);
615:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII);
616:   PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII);
617:   return(0);
618: }