Actual source code: rii.c
slepc-3.21.0 2024-03-30
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: }