Actual source code: rii.c

slepc-3.20.1 2023-11-27
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: "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>
 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:   PetscBool herm;             /* whether the Hermitian version of the scalar equation must be used */
 33:   PetscReal deftol;           /* tolerance for the deflation (threshold) */
 34:   KSP       ksp;              /* linear solver object */
 35: } NEP_RII;

 37: static PetscErrorCode NEPSetUp_RII(NEP nep)
 38: {
 39:   PetscFunctionBegin;
 40:   if (nep->ncv!=PETSC_DEFAULT) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
 41:   nep->ncv = nep->nev;
 42:   if (nep->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
 43:   nep->mpd = nep->nev;
 44:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
 45:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
 46:   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
 47:   NEPCheckUnsupported(nep,NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
 48:   PetscCall(NEPAllocateSolution(nep,0));
 49:   PetscCall(NEPSetWorkVecs(nep,2));
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: static PetscErrorCode NEPSolve_RII(NEP nep)
 54: {
 55:   NEP_RII            *ctx = (NEP_RII*)nep->data;
 56:   Mat                T,Tp,H,A;
 57:   Vec                uu,u,r,delta,t;
 58:   PetscScalar        lambda,lambda2,sigma,a1,a2,corr;
 59:   PetscReal          nrm,resnorm=1.0,ktol=0.1,perr,rtol;
 60:   PetscBool          skip=PETSC_FALSE,lock=PETSC_FALSE;
 61:   PetscInt           inner_its,its=0;
 62:   NEP_EXT_OP         extop=NULL;
 63:   KSPConvergedReason kspreason;

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

 83:   /* prepare linear solver */
 84:   PetscCall(NEPDeflationSolveSetUp(extop,sigma));
 85:   PetscCall(KSPGetTolerances(ctx->ksp,&rtol,NULL,NULL,NULL));

 87:   PetscCall(VecCopy(u,r));
 88:   PetscCall(NEPDeflationFunctionSolve(extop,r,u));
 89:   PetscCall(VecNorm(u,NORM_2,&nrm));
 90:   PetscCall(VecScale(u,1.0/nrm));

 92:   /* Restart loop */
 93:   while (nep->reason == NEP_CONVERGED_ITERATING) {
 94:     its++;

 96:     /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
 97:        estimate as starting value. */
 98:     inner_its=0;
 99:     lambda2 = lambda;
100:     do {
101:       PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
102:       PetscCall(MatMult(T,u,r));
103:       if (!ctx->herm) {
104:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
105:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
106:         if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
107:         t = delta;
108:       } else t = r;
109:       PetscCall(VecDot(t,u,&a1));
110:       PetscCall(NEPDeflationComputeJacobian(extop,lambda,&Tp));
111:       PetscCall(MatMult(Tp,u,r));
112:       if (!ctx->herm) {
113:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
114:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
115:         if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
116:         t = delta;
117:       } else t = r;
118:       PetscCall(VecDot(t,u,&a2));
119:       corr = a1/a2;
120:       lambda = lambda - corr;
121:       inner_its++;
122:     } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);

124:     /* form residual,  r = T(lambda)*u */
125:     PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
126:     PetscCall(MatMult(T,u,r));

128:     /* convergence test */
129:     perr = nep->errest[nep->nconv];
130:     PetscCall(VecNorm(r,NORM_2,&resnorm));
131:     PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
132:     nep->eigr[nep->nconv] = lambda;
133:     if (its>1 && (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol)) {
134:       if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
135:         PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
136:         PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
137:         PetscCall(MatMult(T,u,r));
138:         PetscCall(VecNorm(r,NORM_2,&resnorm));
139:         PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
140:         if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
141:       } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
142:     }
143:     if (lock) {
144:       PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
145:       nep->nconv = nep->nconv + 1;
146:       PetscCall(NEPDeflationLocking(extop,u,lambda));
147:       nep->its += its;
148:       skip = PETSC_TRUE;
149:       lock = PETSC_FALSE;
150:     }
151:     PetscCall((*nep->stopping)(nep,nep->its+its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
152:     if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its+its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));

154:     if (nep->reason == NEP_CONVERGED_ITERATING) {
155:       if (!skip) {
156:         /* update preconditioner and set adaptive tolerance */
157:         if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) PetscCall(NEPDeflationSolveSetUp(extop,lambda2));
158:         if (!ctx->cctol) {
159:           ktol = PetscMax(ktol/2.0,rtol);
160:           PetscCall(KSPSetTolerances(ctx->ksp,ktol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
161:         }

163:         /* eigenvector correction: delta = T(sigma)\r */
164:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
165:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
166:         if (kspreason<0) {
167:           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its));
168:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
169:           break;
170:         }

172:         /* update eigenvector: u = u - delta */
173:         PetscCall(VecAXPY(u,-1.0,delta));

175:         /* normalize eigenvector */
176:         PetscCall(VecNormalize(u,NULL));
177:       } else {
178:         its = -1;
179:         PetscCall(NEPDeflationSetRandomVec(extop,u));
180:         PetscCall(NEPDeflationSolveSetUp(extop,sigma));
181:         PetscCall(VecCopy(u,r));
182:         PetscCall(NEPDeflationFunctionSolve(extop,r,u));
183:         PetscCall(VecNorm(u,NORM_2,&nrm));
184:         PetscCall(VecScale(u,1.0/nrm));
185:         lambda = sigma;
186:         skip = PETSC_FALSE;
187:       }
188:     }
189:   }
190:   PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
191:   PetscCall(DSSetType(nep->ds,DSNHEP));
192:   PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
193:   PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
194:   PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
195:   PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
196:   PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
197:   PetscCall(MatDestroy(&H));
198:   PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
199:   PetscCall(NEPDeflationReset(extop));
200:   PetscCall(VecDestroy(&u));
201:   PetscCall(VecDestroy(&r));
202:   PetscCall(VecDestroy(&delta));
203:   PetscFunctionReturn(PETSC_SUCCESS);
204: }

206: static PetscErrorCode NEPSetFromOptions_RII(NEP nep,PetscOptionItems *PetscOptionsObject)
207: {
208:   NEP_RII        *ctx = (NEP_RII*)nep->data;
209:   PetscBool      flg;
210:   PetscInt       i;
211:   PetscReal      r;

213:   PetscFunctionBegin;
214:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP RII Options");

216:     i = 0;
217:     PetscCall(PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&i,&flg));
218:     if (flg) PetscCall(NEPRIISetMaximumIterations(nep,i));

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

222:     PetscCall(PetscOptionsBool("-nep_rii_hermitian","Use Hermitian version of the scalar nonlinear equation","NEPRIISetHermitian",ctx->herm,&ctx->herm,NULL));

224:     i = 0;
225:     PetscCall(PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg));
226:     if (flg) PetscCall(NEPRIISetLagPreconditioner(nep,i));

228:     r = 0.0;
229:     PetscCall(PetscOptionsReal("-nep_rii_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPRIISetDeflationThreshold",ctx->deftol,&r,&flg));
230:     if (flg) PetscCall(NEPRIISetDeflationThreshold(nep,r));

232:   PetscOptionsHeadEnd();

234:   if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
235:   PetscCall(KSPSetFromOptions(ctx->ksp));
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
240: {
241:   NEP_RII *ctx = (NEP_RII*)nep->data;

243:   PetscFunctionBegin;
244:   if (its==PETSC_DEFAULT) ctx->max_inner_it = 10;
245:   else {
246:     PetscCheck(its>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations must be >0");
247:     ctx->max_inner_it = its;
248:   }
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: /*@
253:    NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
254:    used in the RII solver. These are the Newton iterations related to the computation
255:    of the nonlinear Rayleigh functional.

257:    Logically Collective

259:    Input Parameters:
260: +  nep - nonlinear eigenvalue solver
261: -  its - maximum inner iterations

263:    Level: advanced

265: .seealso: NEPRIIGetMaximumIterations()
266: @*/
267: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
268: {
269:   PetscFunctionBegin;
272:   PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
277: {
278:   NEP_RII *ctx = (NEP_RII*)nep->data;

280:   PetscFunctionBegin;
281:   *its = ctx->max_inner_it;
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: /*@
286:    NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.

288:    Not Collective

290:    Input Parameter:
291: .  nep - nonlinear eigenvalue solver

293:    Output Parameter:
294: .  its - maximum inner iterations

296:    Level: advanced

298: .seealso: NEPRIISetMaximumIterations()
299: @*/
300: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
301: {
302:   PetscFunctionBegin;
304:   PetscAssertPointer(its,2);
305:   PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }

309: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
310: {
311:   NEP_RII *ctx = (NEP_RII*)nep->data;

313:   PetscFunctionBegin;
314:   PetscCheck(lag>=0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
315:   ctx->lag = lag;
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

319: /*@
320:    NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
321:    nonlinear solve.

323:    Logically Collective

325:    Input Parameters:
326: +  nep - nonlinear eigenvalue solver
327: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
328:           computed within the nonlinear iteration, 2 means every second time
329:           the Jacobian is built, etc.

331:    Options Database Keys:
332: .  -nep_rii_lag_preconditioner <lag> - the lag value

334:    Notes:
335:    The default is 1.
336:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

338:    Level: intermediate

340: .seealso: NEPRIIGetLagPreconditioner()
341: @*/
342: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
343: {
344:   PetscFunctionBegin;
347:   PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
352: {
353:   NEP_RII *ctx = (NEP_RII*)nep->data;

355:   PetscFunctionBegin;
356:   *lag = ctx->lag;
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: /*@
361:    NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

363:    Not Collective

365:    Input Parameter:
366: .  nep - nonlinear eigenvalue solver

368:    Output Parameter:
369: .  lag - the lag parameter

371:    Level: intermediate

373: .seealso: NEPRIISetLagPreconditioner()
374: @*/
375: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
376: {
377:   PetscFunctionBegin;
379:   PetscAssertPointer(lag,2);
380:   PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
385: {
386:   NEP_RII *ctx = (NEP_RII*)nep->data;

388:   PetscFunctionBegin;
389:   ctx->cctol = cct;
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }

393: /*@
394:    NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
395:    in the linear solver constant.

397:    Logically Collective

399:    Input Parameters:
400: +  nep - nonlinear eigenvalue solver
401: -  cct - a boolean value

403:    Options Database Keys:
404: .  -nep_rii_const_correction_tol <bool> - set the boolean flag

406:    Notes:
407:    By default, an exponentially decreasing tolerance is set in the KSP used
408:    within the nonlinear iteration, so that each Newton iteration requests
409:    better accuracy than the previous one. The constant correction tolerance
410:    flag stops this behaviour.

412:    Level: intermediate

414: .seealso: NEPRIIGetConstCorrectionTol()
415: @*/
416: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
417: {
418:   PetscFunctionBegin;
421:   PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
422:   PetscFunctionReturn(PETSC_SUCCESS);
423: }

425: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
426: {
427:   NEP_RII *ctx = (NEP_RII*)nep->data;

429:   PetscFunctionBegin;
430:   *cct = ctx->cctol;
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: /*@
435:    NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.

437:    Not Collective

439:    Input Parameter:
440: .  nep - nonlinear eigenvalue solver

442:    Output Parameter:
443: .  cct - the value of the constant tolerance flag

445:    Level: intermediate

447: .seealso: NEPRIISetConstCorrectionTol()
448: @*/
449: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
450: {
451:   PetscFunctionBegin;
453:   PetscAssertPointer(cct,2);
454:   PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: static PetscErrorCode NEPRIISetHermitian_RII(NEP nep,PetscBool herm)
459: {
460:   NEP_RII *ctx = (NEP_RII*)nep->data;

462:   PetscFunctionBegin;
463:   ctx->herm = herm;
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: /*@
468:    NEPRIISetHermitian - Sets a flag to indicate if the Hermitian version of the
469:    scalar nonlinear equation must be used by the solver.

471:    Logically Collective

473:    Input Parameters:
474: +  nep  - nonlinear eigenvalue solver
475: -  herm - a boolean value

477:    Options Database Keys:
478: .  -nep_rii_hermitian <bool> - set the boolean flag

480:    Notes:
481:    By default, the scalar nonlinear equation x'*inv(T(sigma))*T(z)*x=0 is solved
482:    at each step of the nonlinear iteration. When this flag is set the simpler
483:    form x'*T(z)*x=0 is used, which is supposed to be valid only for Hermitian
484:    problems.

486:    Level: intermediate

488: .seealso: NEPRIIGetHermitian()
489: @*/
490: PetscErrorCode NEPRIISetHermitian(NEP nep,PetscBool herm)
491: {
492:   PetscFunctionBegin;
495:   PetscTryMethod(nep,"NEPRIISetHermitian_C",(NEP,PetscBool),(nep,herm));
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

499: static PetscErrorCode NEPRIIGetHermitian_RII(NEP nep,PetscBool *herm)
500: {
501:   NEP_RII *ctx = (NEP_RII*)nep->data;

503:   PetscFunctionBegin;
504:   *herm = ctx->herm;
505:   PetscFunctionReturn(PETSC_SUCCESS);
506: }

508: /*@
509:    NEPRIIGetHermitian - Returns the flag about using the Hermitian version of
510:    the scalar nonlinear equation.

512:    Not Collective

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

517:    Output Parameter:
518: .  herm - the value of the hermitian flag

520:    Level: intermediate

522: .seealso: NEPRIISetHermitian()
523: @*/
524: PetscErrorCode NEPRIIGetHermitian(NEP nep,PetscBool *herm)
525: {
526:   PetscFunctionBegin;
528:   PetscAssertPointer(herm,2);
529:   PetscUseMethod(nep,"NEPRIIGetHermitian_C",(NEP,PetscBool*),(nep,herm));
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: static PetscErrorCode NEPRIISetDeflationThreshold_RII(NEP nep,PetscReal deftol)
534: {
535:   NEP_RII *ctx = (NEP_RII*)nep->data;

537:   PetscFunctionBegin;
538:   ctx->deftol = deftol;
539:   PetscFunctionReturn(PETSC_SUCCESS);
540: }

542: /*@
543:    NEPRIISetDeflationThreshold - Sets the threshold value used to switch between
544:    deflated and non-deflated iteration.

546:    Logically Collective

548:    Input Parameters:
549: +  nep    - nonlinear eigenvalue solver
550: -  deftol - the threshold value

552:    Options Database Keys:
553: .  -nep_rii_deflation_threshold <deftol> - set the threshold

555:    Notes:
556:    Normally, the solver iterates on the extended problem in order to deflate
557:    previously converged eigenpairs. If this threshold is set to a nonzero value,
558:    then once the residual error is below this threshold the solver will
559:    continue the iteration without deflation. The intention is to be able to
560:    improve the current eigenpair further, despite having previous eigenpairs
561:    with somewhat bad precision.

563:    Level: advanced

565: .seealso: NEPRIIGetDeflationThreshold()
566: @*/
567: PetscErrorCode NEPRIISetDeflationThreshold(NEP nep,PetscReal deftol)
568: {
569:   PetscFunctionBegin;
572:   PetscTryMethod(nep,"NEPRIISetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: static PetscErrorCode NEPRIIGetDeflationThreshold_RII(NEP nep,PetscReal *deftol)
577: {
578:   NEP_RII *ctx = (NEP_RII*)nep->data;

580:   PetscFunctionBegin;
581:   *deftol = ctx->deftol;
582:   PetscFunctionReturn(PETSC_SUCCESS);
583: }

585: /*@
586:    NEPRIIGetDeflationThreshold - Returns the threshold value that controls deflation.

588:    Not Collective

590:    Input Parameter:
591: .  nep - nonlinear eigenvalue solver

593:    Output Parameter:
594: .  deftol - the threshold

596:    Level: advanced

598: .seealso: NEPRIISetDeflationThreshold()
599: @*/
600: PetscErrorCode NEPRIIGetDeflationThreshold(NEP nep,PetscReal *deftol)
601: {
602:   PetscFunctionBegin;
604:   PetscAssertPointer(deftol,2);
605:   PetscUseMethod(nep,"NEPRIIGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
606:   PetscFunctionReturn(PETSC_SUCCESS);
607: }

609: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
610: {
611:   NEP_RII        *ctx = (NEP_RII*)nep->data;

613:   PetscFunctionBegin;
614:   PetscCall(PetscObjectReference((PetscObject)ksp));
615:   PetscCall(KSPDestroy(&ctx->ksp));
616:   ctx->ksp = ksp;
617:   nep->state = NEP_STATE_INITIAL;
618:   PetscFunctionReturn(PETSC_SUCCESS);
619: }

621: /*@
622:    NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
623:    eigenvalue solver.

625:    Collective

627:    Input Parameters:
628: +  nep - eigenvalue solver
629: -  ksp - the linear solver object

631:    Level: advanced

633: .seealso: NEPRIIGetKSP()
634: @*/
635: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
636: {
637:   PetscFunctionBegin;
640:   PetscCheckSameComm(nep,1,ksp,2);
641:   PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
646: {
647:   NEP_RII        *ctx = (NEP_RII*)nep->data;

649:   PetscFunctionBegin;
650:   if (!ctx->ksp) {
651:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
652:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
653:     PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
654:     PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_"));
655:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
656:     PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
657:     PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
658:   }
659:   *ksp = ctx->ksp;
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: /*@
664:    NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
665:    the nonlinear eigenvalue solver.

667:    Collective

669:    Input Parameter:
670: .  nep - nonlinear eigenvalue solver

672:    Output Parameter:
673: .  ksp - the linear solver object

675:    Level: advanced

677: .seealso: NEPRIISetKSP()
678: @*/
679: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
680: {
681:   PetscFunctionBegin;
683:   PetscAssertPointer(ksp,2);
684:   PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
685:   PetscFunctionReturn(PETSC_SUCCESS);
686: }

688: static PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
689: {
690:   NEP_RII        *ctx = (NEP_RII*)nep->data;
691:   PetscBool      isascii;

693:   PetscFunctionBegin;
694:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
695:   if (isascii) {
696:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum number of inner iterations: %" PetscInt_FMT "\n",ctx->max_inner_it));
697:     if (ctx->cctol) PetscCall(PetscViewerASCIIPrintf(viewer,"  using a constant tolerance for the linear solver\n"));
698:     if (ctx->herm) PetscCall(PetscViewerASCIIPrintf(viewer,"  using the Hermitian version of the scalar nonlinear equation\n"));
699:     if (ctx->lag) PetscCall(PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag));
700:     if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer,"  deflation threshold: %g\n",(double)ctx->deftol));
701:     if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
702:     PetscCall(PetscViewerASCIIPushTab(viewer));
703:     PetscCall(KSPView(ctx->ksp,viewer));
704:     PetscCall(PetscViewerASCIIPopTab(viewer));
705:   }
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

709: static PetscErrorCode NEPReset_RII(NEP nep)
710: {
711:   NEP_RII        *ctx = (NEP_RII*)nep->data;

713:   PetscFunctionBegin;
714:   PetscCall(KSPReset(ctx->ksp));
715:   PetscFunctionReturn(PETSC_SUCCESS);
716: }

718: static PetscErrorCode NEPDestroy_RII(NEP nep)
719: {
720:   NEP_RII        *ctx = (NEP_RII*)nep->data;

722:   PetscFunctionBegin;
723:   PetscCall(KSPDestroy(&ctx->ksp));
724:   PetscCall(PetscFree(nep->data));
725:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL));
726:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL));
727:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL));
728:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL));
729:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL));
730:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL));
731:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NULL));
732:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NULL));
733:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NULL));
734:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NULL));
735:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL));
736:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL));
737:   PetscFunctionReturn(PETSC_SUCCESS);
738: }

740: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
741: {
742:   NEP_RII        *ctx;

744:   PetscFunctionBegin;
745:   PetscCall(PetscNew(&ctx));
746:   nep->data = (void*)ctx;
747:   ctx->max_inner_it = 10;
748:   ctx->lag          = 1;
749:   ctx->cctol        = PETSC_FALSE;
750:   ctx->herm         = PETSC_FALSE;
751:   ctx->deftol       = 0.0;

753:   nep->useds = PETSC_TRUE;

755:   nep->ops->solve          = NEPSolve_RII;
756:   nep->ops->setup          = NEPSetUp_RII;
757:   nep->ops->setfromoptions = NEPSetFromOptions_RII;
758:   nep->ops->reset          = NEPReset_RII;
759:   nep->ops->destroy        = NEPDestroy_RII;
760:   nep->ops->view           = NEPView_RII;
761:   nep->ops->computevectors = NEPComputeVectors_Schur;

763:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII));
764:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII));
765:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII));
766:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII));
767:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII));
768:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII));
769:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NEPRIISetHermitian_RII));
770:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NEPRIIGetHermitian_RII));
771:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NEPRIISetDeflationThreshold_RII));
772:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NEPRIIGetDeflationThreshold_RII));
773:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII));
774:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII));
775:   PetscFunctionReturn(PETSC_SUCCESS);
776: }