Actual source code: rii.c
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_DETERMINE) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
41: nep->ncv = nep->nev;
42: if (nep->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
43: nep->mpd = nep->nev;
44: if (nep->max_it==PETSC_DETERMINE) 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_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
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 || its==PETSC_DECIDE) 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 - the nonlinear eigensolver context
261: - its - maximum inner iterations
263: Options Database Key:
264: . -nep_rii_max_it its - set the maximum number of inner iterations
266: Level: advanced
268: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetMaximumIterations()`
269: @*/
270: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
271: {
272: PetscFunctionBegin;
275: PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
280: {
281: NEP_RII *ctx = (NEP_RII*)nep->data;
283: PetscFunctionBegin;
284: *its = ctx->max_inner_it;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.
291: Not Collective
293: Input Parameter:
294: . nep - the nonlinear eigensolver context
296: Output Parameter:
297: . its - maximum inner iterations
299: Level: advanced
301: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetMaximumIterations()`
302: @*/
303: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
304: {
305: PetscFunctionBegin;
307: PetscAssertPointer(its,2);
308: PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
313: {
314: NEP_RII *ctx = (NEP_RII*)nep->data;
316: PetscFunctionBegin;
317: PetscCheck(lag>=0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
318: ctx->lag = lag;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
324: nonlinear solve.
326: Logically Collective
328: Input Parameters:
329: + nep - the nonlinear eigensolver context
330: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
331: computed within the nonlinear iteration, 2 means every second time
332: the Jacobian is built, etc.
334: Options Database Key:
335: . -nep_rii_lag_preconditioner lag - the lag value
337: Notes:
338: The default is 1.
339: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
341: Level: intermediate
343: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetLagPreconditioner()`
344: @*/
345: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
346: {
347: PetscFunctionBegin;
350: PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
355: {
356: NEP_RII *ctx = (NEP_RII*)nep->data;
358: PetscFunctionBegin;
359: *lag = ctx->lag;
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /*@
364: NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
366: Not Collective
368: Input Parameter:
369: . nep - the nonlinear eigensolver context
371: Output Parameter:
372: . lag - the lag parameter
374: Level: intermediate
376: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetLagPreconditioner()`
377: @*/
378: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
379: {
380: PetscFunctionBegin;
382: PetscAssertPointer(lag,2);
383: PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
388: {
389: NEP_RII *ctx = (NEP_RII*)nep->data;
391: PetscFunctionBegin;
392: ctx->cctol = cct;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
398: in the linear solver constant.
400: Logically Collective
402: Input Parameters:
403: + nep - the nonlinear eigensolver context
404: - cct - if `PETSC_FALSE`, the `KSP` relative tolerance is set to $2^{-i}$
406: Options Database Key:
407: . -nep_rii_const_correction_tol (true|false) - set a constant or dynamic stopping criterion
409: Notes:
410: By default, an exponentially decreasing tolerance is set in the `KSP` used
411: within the nonlinear iteration, so that each Newton iteration requests
412: better accuracy than the previous one. The constant correction tolerance
413: flag stops this behavior.
415: Level: intermediate
417: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetConstCorrectionTol()`, `NEPRIIGetKSP()`
418: @*/
419: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
420: {
421: PetscFunctionBegin;
424: PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
429: {
430: NEP_RII *ctx = (NEP_RII*)nep->data;
432: PetscFunctionBegin;
433: *cct = ctx->cctol;
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: /*@
438: NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.
440: Not Collective
442: Input Parameter:
443: . nep - the nonlinear eigensolver context
445: Output Parameter:
446: . cct - the value of the constant tolerance flag
448: Level: intermediate
450: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetConstCorrectionTol()`
451: @*/
452: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
453: {
454: PetscFunctionBegin;
456: PetscAssertPointer(cct,2);
457: PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode NEPRIISetHermitian_RII(NEP nep,PetscBool herm)
462: {
463: NEP_RII *ctx = (NEP_RII*)nep->data;
465: PetscFunctionBegin;
466: ctx->herm = herm;
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: NEPRIISetHermitian - Sets a flag to indicate if the Hermitian version of the
472: scalar nonlinear equation must be used by the solver.
474: Logically Collective
476: Input Parameters:
477: + nep - the nonlinear eigensolver context
478: - herm - `PETSC_TRUE` if the Hermitian version is preferred
480: Options Database Key:
481: . -nep_rii_hermitian (true|false) - toggle the Hermitian version
483: Notes:
484: By default, the scalar nonlinear equation $x^*T(\sigma)^{-1}T(z)x=0$ is solved
485: at each step of the nonlinear iteration. When this flag is set the simpler
486: form $x^*T(z)x=0$ is used, which is supposed to be valid only for Hermitian
487: problems.
489: Level: intermediate
491: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetHermitian()`
492: @*/
493: PetscErrorCode NEPRIISetHermitian(NEP nep,PetscBool herm)
494: {
495: PetscFunctionBegin;
498: PetscTryMethod(nep,"NEPRIISetHermitian_C",(NEP,PetscBool),(nep,herm));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static PetscErrorCode NEPRIIGetHermitian_RII(NEP nep,PetscBool *herm)
503: {
504: NEP_RII *ctx = (NEP_RII*)nep->data;
506: PetscFunctionBegin;
507: *herm = ctx->herm;
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /*@
512: NEPRIIGetHermitian - Returns the flag about using the Hermitian version of
513: the scalar nonlinear equation.
515: Not Collective
517: Input Parameter:
518: . nep - the nonlinear eigensolver context
520: Output Parameter:
521: . herm - the value of the Hermitian flag
523: Level: intermediate
525: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetHermitian()`
526: @*/
527: PetscErrorCode NEPRIIGetHermitian(NEP nep,PetscBool *herm)
528: {
529: PetscFunctionBegin;
531: PetscAssertPointer(herm,2);
532: PetscUseMethod(nep,"NEPRIIGetHermitian_C",(NEP,PetscBool*),(nep,herm));
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: static PetscErrorCode NEPRIISetDeflationThreshold_RII(NEP nep,PetscReal deftol)
537: {
538: NEP_RII *ctx = (NEP_RII*)nep->data;
540: PetscFunctionBegin;
541: ctx->deftol = deftol;
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: /*@
546: NEPRIISetDeflationThreshold - Sets the threshold value used to switch between
547: deflated and non-deflated iteration.
549: Logically Collective
551: Input Parameters:
552: + nep - the nonlinear eigensolver context
553: - deftol - the threshold value
555: Options Database Key:
556: . -nep_rii_deflation_threshold deftol - set the threshold
558: Note:
559: Normally, the solver iterates on the extended problem in order to deflate
560: previously converged eigenpairs. If this threshold is set to a nonzero value,
561: then once the residual error is below this threshold the solver will
562: continue the iteration without deflation. The intention is to be able to
563: improve the current eigenpair further, despite having previous eigenpairs
564: with somewhat bad precision.
566: Level: advanced
568: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetDeflationThreshold()`
569: @*/
570: PetscErrorCode NEPRIISetDeflationThreshold(NEP nep,PetscReal deftol)
571: {
572: PetscFunctionBegin;
575: PetscTryMethod(nep,"NEPRIISetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode NEPRIIGetDeflationThreshold_RII(NEP nep,PetscReal *deftol)
580: {
581: NEP_RII *ctx = (NEP_RII*)nep->data;
583: PetscFunctionBegin;
584: *deftol = ctx->deftol;
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: /*@
589: NEPRIIGetDeflationThreshold - Returns the threshold value that controls deflation.
591: Not Collective
593: Input Parameter:
594: . nep - the nonlinear eigensolver context
596: Output Parameter:
597: . deftol - the threshold
599: Level: advanced
601: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetDeflationThreshold()`
602: @*/
603: PetscErrorCode NEPRIIGetDeflationThreshold(NEP nep,PetscReal *deftol)
604: {
605: PetscFunctionBegin;
607: PetscAssertPointer(deftol,2);
608: PetscUseMethod(nep,"NEPRIIGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
613: {
614: NEP_RII *ctx = (NEP_RII*)nep->data;
616: PetscFunctionBegin;
617: PetscCall(PetscObjectReference((PetscObject)ksp));
618: PetscCall(KSPDestroy(&ctx->ksp));
619: ctx->ksp = ksp;
620: nep->state = NEP_STATE_INITIAL;
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: /*@
625: NEPRIISetKSP - Associate a linear solver object (`KSP`) to the nonlinear
626: eigenvalue solver.
628: Collective
630: Input Parameters:
631: + nep - the nonlinear eigensolver context
632: - ksp - the linear solver object
634: Level: advanced
636: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetKSP()`
637: @*/
638: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
639: {
640: PetscFunctionBegin;
643: PetscCheckSameComm(nep,1,ksp,2);
644: PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
649: {
650: NEP_RII *ctx = (NEP_RII*)nep->data;
652: PetscFunctionBegin;
653: if (!ctx->ksp) {
654: PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
655: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
656: PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
657: PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_"));
658: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
659: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
660: PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
661: }
662: *ksp = ctx->ksp;
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@
667: NEPRIIGetKSP - Retrieve the linear solver object (`KSP`) associated with
668: the nonlinear eigenvalue solver.
670: Collective
672: Input Parameter:
673: . nep - the nonlinear eigensolver context
675: Output Parameter:
676: . ksp - the linear solver object
678: Level: advanced
680: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetKSP()`
681: @*/
682: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
683: {
684: PetscFunctionBegin;
686: PetscAssertPointer(ksp,2);
687: PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
692: {
693: NEP_RII *ctx = (NEP_RII*)nep->data;
694: PetscBool isascii;
696: PetscFunctionBegin;
697: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
698: if (isascii) {
699: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of inner iterations: %" PetscInt_FMT "\n",ctx->max_inner_it));
700: if (ctx->cctol) PetscCall(PetscViewerASCIIPrintf(viewer," using a constant tolerance for the linear solver\n"));
701: if (ctx->herm) PetscCall(PetscViewerASCIIPrintf(viewer," using the Hermitian version of the scalar nonlinear equation\n"));
702: if (ctx->lag) PetscCall(PetscViewerASCIIPrintf(viewer," updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag));
703: if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer," deflation threshold: %g\n",(double)ctx->deftol));
704: if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
705: PetscCall(PetscViewerASCIIPushTab(viewer));
706: PetscCall(KSPView(ctx->ksp,viewer));
707: PetscCall(PetscViewerASCIIPopTab(viewer));
708: }
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: static PetscErrorCode NEPReset_RII(NEP nep)
713: {
714: NEP_RII *ctx = (NEP_RII*)nep->data;
716: PetscFunctionBegin;
717: PetscCall(KSPReset(ctx->ksp));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode NEPDestroy_RII(NEP nep)
722: {
723: NEP_RII *ctx = (NEP_RII*)nep->data;
725: PetscFunctionBegin;
726: PetscCall(KSPDestroy(&ctx->ksp));
727: PetscCall(PetscFree(nep->data));
728: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL));
729: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL));
730: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL));
731: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL));
732: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL));
733: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL));
734: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NULL));
735: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NULL));
736: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NULL));
737: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NULL));
738: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL));
739: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*MC
744: NEPRII - NEPRII = "rii" - Simple residual inverse iteration with
745: varying shift.
747: Notes:
748: This is the default solver, although it is very basic. Users are
749: advised to try other solvers.
751: This solver is based on the modification proposed by {cite:t}`Neu85`
752: of the classical residual inverse iteration method. At each step,
753: this method has to solve a system of linear equations. Call
754: `NEPRIIGetKSP()` to configure the `KSP` object used for this.
755: When using iterative linear solvers, the preconditioner will be updated
756: in each step or not, depending on `NEPRIISetLagPreconditioner()`, and
757: the tolerance will be constant or not depending on
758: `NEPRIISetConstCorrectionTol()`.
760: The solver incorporates deflation, so that several eigenpairs con be
761: computed. Details of the implementation in SLEPc can be found in
762: {cite:p}`Cam21`.
764: Level: beginner
766: .seealso: [](ch:nep), `NEP`, `NEPType`, `NEPSetType()`, `NEPRIIGetKSP()`, `NEPRIISetLagPreconditioner()`, `NEPRIISetConstCorrectionTol()`
767: M*/
768: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
769: {
770: NEP_RII *ctx;
772: PetscFunctionBegin;
773: PetscCall(PetscNew(&ctx));
774: nep->data = (void*)ctx;
775: ctx->max_inner_it = 10;
776: ctx->lag = 1;
777: ctx->cctol = PETSC_FALSE;
778: ctx->herm = PETSC_FALSE;
779: ctx->deftol = 0.0;
781: nep->useds = PETSC_TRUE;
783: nep->ops->solve = NEPSolve_RII;
784: nep->ops->setup = NEPSetUp_RII;
785: nep->ops->setfromoptions = NEPSetFromOptions_RII;
786: nep->ops->reset = NEPReset_RII;
787: nep->ops->destroy = NEPDestroy_RII;
788: nep->ops->view = NEPView_RII;
789: nep->ops->computevectors = NEPComputeVectors_Schur;
791: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII));
792: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII));
793: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII));
794: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII));
795: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII));
796: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII));
797: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NEPRIISetHermitian_RII));
798: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NEPRIIGetHermitian_RII));
799: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NEPRIISetDeflationThreshold_RII));
800: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NEPRIIGetDeflationThreshold_RII));
801: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII));
802: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII));
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }