Actual source code: rii.c
slepc-3.17.0 2022-03-31
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: PetscErrorCode NEPSetUp_RII(NEP nep)
38: {
39: if (nep->ncv!=PETSC_DEFAULT) PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n");
40: nep->ncv = nep->nev;
41: if (nep->mpd!=PETSC_DEFAULT) PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n");
42: nep->mpd = nep->nev;
43: if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
44: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
46: NEPCheckUnsupported(nep,NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
47: NEPAllocateSolution(nep,0);
48: NEPSetWorkVecs(nep,2);
49: PetscFunctionReturn(0);
50: }
52: PetscErrorCode NEPSolve_RII(NEP nep)
53: {
54: NEP_RII *ctx = (NEP_RII*)nep->data;
55: Mat T,Tp,H;
56: Vec uu,u,r,delta,t;
57: PetscScalar lambda,lambda2,sigma,a1,a2,corr,*Ap;
58: const PetscScalar *Hp;
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,ldh,ldds,i,j;
62: NEP_EXT_OP extop=NULL;
63: KSPConvergedReason kspreason;
65: /* get initial approximation of eigenvalue and eigenvector */
66: NEPGetDefaultShift(nep,&sigma);
67: lambda = sigma;
68: if (!nep->nini) {
69: BVSetRandomColumn(nep->V,0);
70: BVNormColumn(nep->V,0,NORM_2,&nrm);
71: BVScaleColumn(nep->V,0,1.0/nrm);
72: }
73: if (!ctx->ksp) NEPRIIGetKSP(nep,&ctx->ksp);
74: NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
75: NEPDeflationCreateVec(extop,&u);
76: VecDuplicate(u,&r);
77: VecDuplicate(u,&delta);
78: BVGetColumn(nep->V,0,&uu);
79: NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE);
80: BVRestoreColumn(nep->V,0,&uu);
82: /* prepare linear solver */
83: NEPDeflationSolveSetUp(extop,sigma);
84: KSPGetTolerances(ctx->ksp,&rtol,NULL,NULL,NULL);
86: VecCopy(u,r);
87: NEPDeflationFunctionSolve(extop,r,u);
88: VecNorm(u,NORM_2,&nrm);
89: VecScale(u,1.0/nrm);
91: /* Restart loop */
92: while (nep->reason == NEP_CONVERGED_ITERATING) {
93: its++;
95: /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
96: estimate as starting value. */
97: inner_its=0;
98: lambda2 = lambda;
99: do {
100: NEPDeflationComputeFunction(extop,lambda,&T);
101: MatMult(T,u,r);
102: if (!ctx->herm) {
103: NEPDeflationFunctionSolve(extop,r,delta);
104: KSPGetConvergedReason(ctx->ksp,&kspreason);
105: if (kspreason<0) PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its);
106: t = delta;
107: } else t = r;
108: VecDot(t,u,&a1);
109: NEPDeflationComputeJacobian(extop,lambda,&Tp);
110: MatMult(Tp,u,r);
111: if (!ctx->herm) {
112: NEPDeflationFunctionSolve(extop,r,delta);
113: KSPGetConvergedReason(ctx->ksp,&kspreason);
114: if (kspreason<0) PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its);
115: t = delta;
116: } else t = r;
117: VecDot(t,u,&a2);
118: corr = a1/a2;
119: lambda = lambda - corr;
120: inner_its++;
121: } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);
123: /* form residual, r = T(lambda)*u */
124: NEPDeflationComputeFunction(extop,lambda,&T);
125: MatMult(T,u,r);
127: /* convergence test */
128: perr = nep->errest[nep->nconv];
129: VecNorm(r,NORM_2,&resnorm);
130: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
131: nep->eigr[nep->nconv] = lambda;
132: if (its>1 && (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol)) {
133: if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
134: NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds);
135: NEPDeflationSetRefine(extop,PETSC_TRUE);
136: MatMult(T,u,r);
137: VecNorm(r,NORM_2,&resnorm);
138: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
139: if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
140: } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
141: }
142: if (lock) {
143: NEPDeflationSetRefine(extop,PETSC_FALSE);
144: nep->nconv = nep->nconv + 1;
145: NEPDeflationLocking(extop,u,lambda);
146: nep->its += its;
147: skip = PETSC_TRUE;
148: lock = PETSC_FALSE;
149: }
150: (*nep->stopping)(nep,nep->its+its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
151: if (!skip || nep->reason>0) NEPMonitor(nep,nep->its+its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
153: if (nep->reason == NEP_CONVERGED_ITERATING) {
154: if (!skip) {
155: /* update preconditioner and set adaptive tolerance */
156: if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) NEPDeflationSolveSetUp(extop,lambda2);
157: if (!ctx->cctol) {
158: ktol = PetscMax(ktol/2.0,rtol);
159: KSPSetTolerances(ctx->ksp,ktol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
160: }
162: /* eigenvector correction: delta = T(sigma)\r */
163: NEPDeflationFunctionSolve(extop,r,delta);
164: KSPGetConvergedReason(ctx->ksp,&kspreason);
165: if (kspreason<0) {
166: PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its);
167: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
168: break;
169: }
171: /* update eigenvector: u = u - delta */
172: VecAXPY(u,-1.0,delta);
174: /* normalize eigenvector */
175: VecNormalize(u,NULL);
176: } else {
177: its = -1;
178: NEPDeflationSetRandomVec(extop,u);
179: NEPDeflationSolveSetUp(extop,sigma);
180: VecCopy(u,r);
181: NEPDeflationFunctionSolve(extop,r,u);
182: VecNorm(u,NORM_2,&nrm);
183: VecScale(u,1.0/nrm);
184: lambda = sigma;
185: skip = PETSC_FALSE;
186: }
187: }
188: }
189: NEPDeflationGetInvariantPair(extop,NULL,&H);
190: MatGetSize(H,NULL,&ldh);
191: DSSetType(nep->ds,DSNHEP);
192: DSReset(nep->ds);
193: DSAllocate(nep->ds,PetscMax(nep->nconv,1));
194: DSGetLeadingDimension(nep->ds,&ldds);
195: MatDenseGetArrayRead(H,&Hp);
196: DSGetArray(nep->ds,DS_MAT_A,&Ap);
197: for (j=0;j<nep->nconv;j++)
198: for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
199: DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
200: MatDenseRestoreArrayRead(H,&Hp);
201: MatDestroy(&H);
202: DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv);
203: DSSolve(nep->ds,nep->eigr,nep->eigi);
204: NEPDeflationReset(extop);
205: VecDestroy(&u);
206: VecDestroy(&r);
207: VecDestroy(&delta);
208: PetscFunctionReturn(0);
209: }
211: PetscErrorCode NEPSetFromOptions_RII(PetscOptionItems *PetscOptionsObject,NEP nep)
212: {
213: NEP_RII *ctx = (NEP_RII*)nep->data;
214: PetscBool flg;
215: PetscInt i;
216: PetscReal r;
218: PetscOptionsHead(PetscOptionsObject,"NEP RII Options");
220: i = 0;
221: PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&i,&flg);
222: if (flg) NEPRIISetMaximumIterations(nep,i);
224: PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL);
226: PetscOptionsBool("-nep_rii_hermitian","Use Hermitian version of the scalar nonlinear equation","NEPRIISetHermitian",ctx->herm,&ctx->herm,NULL);
228: i = 0;
229: PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg);
230: if (flg) NEPRIISetLagPreconditioner(nep,i);
232: r = 0.0;
233: PetscOptionsReal("-nep_rii_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPRIISetDeflationThreshold",ctx->deftol,&r,&flg);
234: if (flg) NEPRIISetDeflationThreshold(nep,r);
236: PetscOptionsTail();
238: if (!ctx->ksp) NEPRIIGetKSP(nep,&ctx->ksp);
239: KSPSetFromOptions(ctx->ksp);
240: PetscFunctionReturn(0);
241: }
243: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
244: {
245: NEP_RII *ctx = (NEP_RII*)nep->data;
247: if (its==PETSC_DEFAULT) ctx->max_inner_it = 10;
248: else {
250: ctx->max_inner_it = its;
251: }
252: PetscFunctionReturn(0);
253: }
255: /*@
256: NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
257: used in the RII solver. These are the Newton iterations related to the computation
258: of the nonlinear Rayleigh functional.
260: Logically Collective on nep
262: Input Parameters:
263: + nep - nonlinear eigenvalue solver
264: - its - maximum inner iterations
266: Level: advanced
268: .seealso: NEPRIIGetMaximumIterations()
269: @*/
270: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
271: {
274: PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
275: PetscFunctionReturn(0);
276: }
278: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
279: {
280: NEP_RII *ctx = (NEP_RII*)nep->data;
282: *its = ctx->max_inner_it;
283: PetscFunctionReturn(0);
284: }
286: /*@
287: NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.
289: Not Collective
291: Input Parameter:
292: . nep - nonlinear eigenvalue solver
294: Output Parameter:
295: . its - maximum inner iterations
297: Level: advanced
299: .seealso: NEPRIISetMaximumIterations()
300: @*/
301: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
302: {
305: PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
306: PetscFunctionReturn(0);
307: }
309: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
310: {
311: NEP_RII *ctx = (NEP_RII*)nep->data;
314: ctx->lag = lag;
315: PetscFunctionReturn(0);
316: }
318: /*@
319: NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
320: nonlinear solve.
322: Logically Collective on nep
324: Input Parameters:
325: + nep - nonlinear eigenvalue solver
326: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
327: computed within the nonlinear iteration, 2 means every second time
328: the Jacobian is built, etc.
330: Options Database Keys:
331: . -nep_rii_lag_preconditioner <lag> - the lag value
333: Notes:
334: The default is 1.
335: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
337: Level: intermediate
339: .seealso: NEPRIIGetLagPreconditioner()
340: @*/
341: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
342: {
345: PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
346: PetscFunctionReturn(0);
347: }
349: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
350: {
351: NEP_RII *ctx = (NEP_RII*)nep->data;
353: *lag = ctx->lag;
354: PetscFunctionReturn(0);
355: }
357: /*@
358: NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
360: Not Collective
362: Input Parameter:
363: . nep - nonlinear eigenvalue solver
365: Output Parameter:
366: . lag - the lag parameter
368: Level: intermediate
370: .seealso: NEPRIISetLagPreconditioner()
371: @*/
372: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
373: {
376: PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
377: PetscFunctionReturn(0);
378: }
380: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
381: {
382: NEP_RII *ctx = (NEP_RII*)nep->data;
384: ctx->cctol = cct;
385: PetscFunctionReturn(0);
386: }
388: /*@
389: NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
390: in the linear solver constant.
392: Logically Collective on nep
394: Input Parameters:
395: + nep - nonlinear eigenvalue solver
396: - cct - a boolean value
398: Options Database Keys:
399: . -nep_rii_const_correction_tol <bool> - set the boolean flag
401: Notes:
402: By default, an exponentially decreasing tolerance is set in the KSP used
403: within the nonlinear iteration, so that each Newton iteration requests
404: better accuracy than the previous one. The constant correction tolerance
405: flag stops this behaviour.
407: Level: intermediate
409: .seealso: NEPRIIGetConstCorrectionTol()
410: @*/
411: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
412: {
415: PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
416: PetscFunctionReturn(0);
417: }
419: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
420: {
421: NEP_RII *ctx = (NEP_RII*)nep->data;
423: *cct = ctx->cctol;
424: PetscFunctionReturn(0);
425: }
427: /*@
428: NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.
430: Not Collective
432: Input Parameter:
433: . nep - nonlinear eigenvalue solver
435: Output Parameter:
436: . cct - the value of the constant tolerance flag
438: Level: intermediate
440: .seealso: NEPRIISetConstCorrectionTol()
441: @*/
442: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
443: {
446: PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
447: PetscFunctionReturn(0);
448: }
450: static PetscErrorCode NEPRIISetHermitian_RII(NEP nep,PetscBool herm)
451: {
452: NEP_RII *ctx = (NEP_RII*)nep->data;
454: ctx->herm = herm;
455: PetscFunctionReturn(0);
456: }
458: /*@
459: NEPRIISetHermitian - Sets a flag to indicate if the Hermitian version of the
460: scalar nonlinear equation must be used by the solver.
462: Logically Collective on nep
464: Input Parameters:
465: + nep - nonlinear eigenvalue solver
466: - herm - a boolean value
468: Options Database Keys:
469: . -nep_rii_hermitian <bool> - set the boolean flag
471: Notes:
472: By default, the scalar nonlinear equation x'*inv(T(sigma))*T(z)*x=0 is solved
473: at each step of the nonlinear iteration. When this flag is set the simpler
474: form x'*T(z)*x=0 is used, which is supposed to be valid only for Hermitian
475: problems.
477: Level: intermediate
479: .seealso: NEPRIIGetHermitian()
480: @*/
481: PetscErrorCode NEPRIISetHermitian(NEP nep,PetscBool herm)
482: {
485: PetscTryMethod(nep,"NEPRIISetHermitian_C",(NEP,PetscBool),(nep,herm));
486: PetscFunctionReturn(0);
487: }
489: static PetscErrorCode NEPRIIGetHermitian_RII(NEP nep,PetscBool *herm)
490: {
491: NEP_RII *ctx = (NEP_RII*)nep->data;
493: *herm = ctx->herm;
494: PetscFunctionReturn(0);
495: }
497: /*@
498: NEPRIIGetHermitian - Returns the flag about using the Hermitian version of
499: the scalar nonlinear equation.
501: Not Collective
503: Input Parameter:
504: . nep - nonlinear eigenvalue solver
506: Output Parameter:
507: . herm - the value of the hermitian flag
509: Level: intermediate
511: .seealso: NEPRIISetHermitian()
512: @*/
513: PetscErrorCode NEPRIIGetHermitian(NEP nep,PetscBool *herm)
514: {
517: PetscUseMethod(nep,"NEPRIIGetHermitian_C",(NEP,PetscBool*),(nep,herm));
518: PetscFunctionReturn(0);
519: }
521: static PetscErrorCode NEPRIISetDeflationThreshold_RII(NEP nep,PetscReal deftol)
522: {
523: NEP_RII *ctx = (NEP_RII*)nep->data;
525: ctx->deftol = deftol;
526: PetscFunctionReturn(0);
527: }
529: /*@
530: NEPRIISetDeflationThreshold - Sets the threshold value used to switch between
531: deflated and non-deflated iteration.
533: Logically Collective on nep
535: Input Parameters:
536: + nep - nonlinear eigenvalue solver
537: - deftol - the threshold value
539: Options Database Keys:
540: . -nep_rii_deflation_threshold <deftol> - set the threshold
542: Notes:
543: Normally, the solver iterates on the extended problem in order to deflate
544: previously converged eigenpairs. If this threshold is set to a nonzero value,
545: then once the residual error is below this threshold the solver will
546: continue the iteration without deflation. The intention is to be able to
547: improve the current eigenpair further, despite having previous eigenpairs
548: with somewhat bad precision.
550: Level: advanced
552: .seealso: NEPRIIGetDeflationThreshold()
553: @*/
554: PetscErrorCode NEPRIISetDeflationThreshold(NEP nep,PetscReal deftol)
555: {
558: PetscTryMethod(nep,"NEPRIISetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
559: PetscFunctionReturn(0);
560: }
562: static PetscErrorCode NEPRIIGetDeflationThreshold_RII(NEP nep,PetscReal *deftol)
563: {
564: NEP_RII *ctx = (NEP_RII*)nep->data;
566: *deftol = ctx->deftol;
567: PetscFunctionReturn(0);
568: }
570: /*@
571: NEPRIIGetDeflationThreshold - Returns the threshold value that controls deflation.
573: Not Collective
575: Input Parameter:
576: . nep - nonlinear eigenvalue solver
578: Output Parameter:
579: . deftol - the threshold
581: Level: advanced
583: .seealso: NEPRIISetDeflationThreshold()
584: @*/
585: PetscErrorCode NEPRIIGetDeflationThreshold(NEP nep,PetscReal *deftol)
586: {
589: PetscUseMethod(nep,"NEPRIIGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
590: PetscFunctionReturn(0);
591: }
593: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
594: {
595: NEP_RII *ctx = (NEP_RII*)nep->data;
597: PetscObjectReference((PetscObject)ksp);
598: KSPDestroy(&ctx->ksp);
599: ctx->ksp = ksp;
600: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
601: nep->state = NEP_STATE_INITIAL;
602: PetscFunctionReturn(0);
603: }
605: /*@
606: NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
607: eigenvalue solver.
609: Collective on nep
611: Input Parameters:
612: + nep - eigenvalue solver
613: - ksp - the linear solver object
615: Level: advanced
617: .seealso: NEPRIIGetKSP()
618: @*/
619: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
620: {
624: PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
625: PetscFunctionReturn(0);
626: }
628: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
629: {
630: NEP_RII *ctx = (NEP_RII*)nep->data;
632: if (!ctx->ksp) {
633: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
634: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
635: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
636: KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_");
637: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
638: PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
639: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
640: KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
641: }
642: *ksp = ctx->ksp;
643: PetscFunctionReturn(0);
644: }
646: /*@
647: NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
648: the nonlinear eigenvalue solver.
650: Not Collective
652: Input Parameter:
653: . nep - nonlinear eigenvalue solver
655: Output Parameter:
656: . ksp - the linear solver object
658: Level: advanced
660: .seealso: NEPRIISetKSP()
661: @*/
662: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
663: {
666: PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
667: PetscFunctionReturn(0);
668: }
670: PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
671: {
672: NEP_RII *ctx = (NEP_RII*)nep->data;
673: PetscBool isascii;
675: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
676: if (isascii) {
677: PetscViewerASCIIPrintf(viewer," maximum number of inner iterations: %" PetscInt_FMT "\n",ctx->max_inner_it);
678: if (ctx->cctol) PetscViewerASCIIPrintf(viewer," using a constant tolerance for the linear solver\n");
679: if (ctx->herm) PetscViewerASCIIPrintf(viewer," using the Hermitian version of the scalar nonlinear equation\n");
680: if (ctx->lag) PetscViewerASCIIPrintf(viewer," updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag);
681: if (ctx->deftol) PetscViewerASCIIPrintf(viewer," deflation threshold: %g\n",(double)ctx->deftol);
682: if (!ctx->ksp) NEPRIIGetKSP(nep,&ctx->ksp);
683: PetscViewerASCIIPushTab(viewer);
684: KSPView(ctx->ksp,viewer);
685: PetscViewerASCIIPopTab(viewer);
686: }
687: PetscFunctionReturn(0);
688: }
690: PetscErrorCode NEPReset_RII(NEP nep)
691: {
692: NEP_RII *ctx = (NEP_RII*)nep->data;
694: KSPReset(ctx->ksp);
695: PetscFunctionReturn(0);
696: }
698: PetscErrorCode NEPDestroy_RII(NEP nep)
699: {
700: NEP_RII *ctx = (NEP_RII*)nep->data;
702: KSPDestroy(&ctx->ksp);
703: PetscFree(nep->data);
704: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL);
705: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL);
706: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL);
707: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL);
708: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL);
709: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL);
710: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NULL);
711: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NULL);
712: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NULL);
713: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NULL);
714: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL);
715: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL);
716: PetscFunctionReturn(0);
717: }
719: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
720: {
721: NEP_RII *ctx;
723: PetscNewLog(nep,&ctx);
724: nep->data = (void*)ctx;
725: ctx->max_inner_it = 10;
726: ctx->lag = 1;
727: ctx->cctol = PETSC_FALSE;
728: ctx->herm = PETSC_FALSE;
729: ctx->deftol = 0.0;
731: nep->useds = PETSC_TRUE;
733: nep->ops->solve = NEPSolve_RII;
734: nep->ops->setup = NEPSetUp_RII;
735: nep->ops->setfromoptions = NEPSetFromOptions_RII;
736: nep->ops->reset = NEPReset_RII;
737: nep->ops->destroy = NEPDestroy_RII;
738: nep->ops->view = NEPView_RII;
739: nep->ops->computevectors = NEPComputeVectors_Schur;
741: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII);
742: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII);
743: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII);
744: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII);
745: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII);
746: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII);
747: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NEPRIISetHermitian_RII);
748: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NEPRIIGetHermitian_RII);
749: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NEPRIISetDeflationThreshold_RII);
750: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NEPRIIGetDeflationThreshold_RII);
751: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII);
752: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII);
753: PetscFunctionReturn(0);
754: }