Actual source code: rii.c
slepc-3.9.2 2018-07-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc nonlinear eigensolver: "rii"
13: Method: Residual inverse iteration
15: Algorithm:
17: Simple residual inverse iteration with varying shift.
19: References:
21: [1] A. Neumaier, "Residual inverse iteration for the nonlinear
22: eigenvalue problem", SIAM J. Numer. Anal. 22(5):914-923, 1985.
23: */
25: #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
27: typedef struct {
28: PetscInt max_inner_it; /* maximum number of Newton iterations */
29: PetscInt lag; /* interval to rebuild preconditioner */
30: PetscBool cctol; /* constant correction tolerance */
31: KSP ksp; /* linear solver object */
32: } NEP_RII;
34: PETSC_STATIC_INLINE PetscErrorCode NEPRII_KSPSolve(NEP nep,Vec b,Vec x)
35: {
37: PetscInt lits;
38: NEP_RII *ctx = (NEP_RII*)nep->data;
41: KSPSolve(ctx->ksp,b,x);
42: KSPGetIterationNumber(ctx->ksp,&lits);
43: PetscInfo2(nep,"iter=%D, linear solve iterations=%D\n",nep->its,lits);
44: return(0);
45: }
47: PetscErrorCode NEPSetUp_RII(NEP nep)
48: {
50: PetscBool istrivial;
53: if (nep->nev>1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Requested several eigenpairs but this solver can compute only one");
54: if (nep->ncv) { PetscInfo(nep,"Setting ncv = 1, ignoring user-provided value\n"); }
55: nep->ncv = 1;
56: if (nep->mpd) { PetscInfo(nep,"Setting mpd = 1, ignoring user-provided value\n"); }
57: nep->mpd = 1;
58: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
59: if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");
61: RGIsTrivial(nep->rg,&istrivial);
62: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");
64: NEPAllocateSolution(nep,0);
65: NEPSetWorkVecs(nep,2);
66: return(0);
67: }
69: PetscErrorCode NEPSolve_RII(NEP nep)
70: {
71: PetscErrorCode ierr;
72: NEP_RII *ctx = (NEP_RII*)nep->data;
73: Mat T=nep->function,Tp=nep->jacobian,Tsigma;
74: Vec u,r=nep->work[0],delta=nep->work[1];
75: PetscScalar lambda,a1,a2,corr;
76: PetscReal resnorm=1.0,ktol=0.1;
77: PetscBool hascopy;
78: PetscInt inner_its;
79: KSPConvergedReason kspreason;
82: /* get initial approximation of eigenvalue and eigenvector */
83: NEPGetDefaultShift(nep,&lambda);
84: if (!nep->nini) {
85: BVSetRandomColumn(nep->V,0);
86: }
87: BVGetColumn(nep->V,0,&u);
88: NEPComputeFunction(nep,lambda,T,T);
90: /* prepare linear solver */
91: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
92: MatDuplicate(T,MAT_COPY_VALUES,&Tsigma);
93: KSPSetOperators(ctx->ksp,Tsigma,Tsigma);
95: /* Restart loop */
96: while (nep->reason == NEP_CONVERGED_ITERATING) {
97: nep->its++;
99: /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
100: estimate as starting value. */
101: inner_its=0;
102: do {
103: NEPApplyFunction(nep,lambda,u,delta,r,T,T);
104: VecDot(r,u,&a1);
105: NEPApplyJacobian(nep,lambda,u,delta,r,Tp);
106: VecDot(r,u,&a2);
107: corr = a1/a2;
108: lambda = lambda - corr;
109: inner_its++;
110: } while (PetscAbsScalar(corr)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);
112: /* update preconditioner and set adaptive tolerance */
113: if (ctx->lag && !(nep->its%ctx->lag) && nep->its>2*ctx->lag && resnorm<1e-2) {
114: MatHasOperation(T,MATOP_COPY,&hascopy);
115: if (hascopy) {
116: MatCopy(T,Tsigma,DIFFERENT_NONZERO_PATTERN);
117: } else {
118: MatDestroy(&Tsigma);
119: MatDuplicate(T,MAT_COPY_VALUES,&Tsigma);
120: }
121: KSPSetOperators(ctx->ksp,Tsigma,Tsigma);
122: }
123: if (!ctx->cctol) {
124: ktol = PetscMax(ktol/2.0,PETSC_MACHINE_EPSILON*10.0);
125: KSPSetTolerances(ctx->ksp,ktol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
126: }
128: /* form residual, r = T(lambda)*u */
129: NEPApplyFunction(nep,lambda,u,delta,r,T,T);
131: /* convergence test */
132: VecNorm(r,NORM_2,&resnorm);
133: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
134: nep->eigr[nep->nconv] = lambda;
135: if (nep->errest[nep->nconv]<=nep->tol) {
136: nep->nconv = nep->nconv + 1;
137: }
138: (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
139: NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,1);
141: if (nep->reason == NEP_CONVERGED_ITERATING) {
142: /* eigenvector correction: delta = T(sigma)\r */
143: NEPRII_KSPSolve(nep,r,delta);
144: KSPGetConvergedReason(ctx->ksp,&kspreason);
145: if (kspreason<0) {
146: PetscInfo1(nep,"iter=%D, linear solve failed, stopping solve\n",nep->its);
147: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
148: break;
149: }
151: /* update eigenvector: u = u - delta */
152: VecAXPY(u,-1.0,delta);
154: /* normalize eigenvector */
155: VecNormalize(u,NULL);
156: }
157: }
158: MatDestroy(&Tsigma);
159: BVRestoreColumn(nep->V,0,&u);
160: return(0);
161: }
163: PetscErrorCode NEPSetFromOptions_RII(PetscOptionItems *PetscOptionsObject,NEP nep)
164: {
166: NEP_RII *ctx = (NEP_RII*)nep->data;
167: PetscBool flg;
168: PetscInt i;
171: PetscOptionsHead(PetscOptionsObject,"NEP RII Options");
173: PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&ctx->max_inner_it,NULL);
175: PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL);
177: i = 0;
178: PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg);
179: if (flg) { NEPRIISetLagPreconditioner(nep,i); }
181: PetscOptionsTail();
183: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
184: KSPSetFromOptions(ctx->ksp);
185: return(0);
186: }
188: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
189: {
190: NEP_RII *ctx = (NEP_RII*)nep->data;
193: ctx->max_inner_it = its;
194: return(0);
195: }
197: /*@
198: NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
199: used in the RII solver. These are the Newton iterations related to the computation
200: of the nonlinear Rayleigh functional.
202: Logically Collective on NEP
204: Input Parameters:
205: + nep - nonlinear eigenvalue solver
206: - its - maximum inner iterations
208: Level: advanced
210: .seealso: NEPRIIGetMaximumIterations()
211: @*/
212: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
213: {
219: PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
220: return(0);
221: }
223: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
224: {
225: NEP_RII *ctx = (NEP_RII*)nep->data;
228: *its = ctx->max_inner_it;
229: return(0);
230: }
232: /*@
233: NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.
235: Not Collective
237: Input Parameter:
238: . nep - nonlinear eigenvalue solver
240: Output Parameter:
241: . its - maximum inner iterations
243: Level: advanced
245: .seealso: NEPRIISetMaximumIterations()
246: @*/
247: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
248: {
254: PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
255: return(0);
256: }
258: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
259: {
260: NEP_RII *ctx = (NEP_RII*)nep->data;
263: if (lag<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
264: ctx->lag = lag;
265: return(0);
266: }
268: /*@
269: NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
270: nonlinear solve.
272: Logically Collective on NEP
274: Input Parameters:
275: + nep - nonlinear eigenvalue solver
276: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
277: computed within the nonlinear iteration, 2 means every second time
278: the Jacobian is built, etc.
280: Options Database Keys:
281: . -nep_rii_lag_preconditioner <lag>
283: Notes:
284: The default is 1.
285: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
287: Level: intermediate
289: .seealso: NEPRIIGetLagPreconditioner()
290: @*/
291: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
292: {
298: PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
299: return(0);
300: }
302: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
303: {
304: NEP_RII *ctx = (NEP_RII*)nep->data;
307: *lag = ctx->lag;
308: return(0);
309: }
311: /*@
312: NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
314: Not Collective
316: Input Parameter:
317: . nep - nonlinear eigenvalue solver
319: Output Parameter:
320: . lag - the lag parameter
322: Level: intermediate
324: .seealso: NEPRIISetLagPreconditioner()
325: @*/
326: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
327: {
333: PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
334: return(0);
335: }
337: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
338: {
339: NEP_RII *ctx = (NEP_RII*)nep->data;
342: ctx->cctol = cct;
343: return(0);
344: }
346: /*@
347: NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
348: in the linear solver constant.
350: Logically Collective on NEP
352: Input Parameters:
353: + nep - nonlinear eigenvalue solver
354: - cct - a boolean value
356: Options Database Keys:
357: . -nep_rii_const_correction_tol <bool> - set the boolean flag
359: Notes:
360: By default, an exponentially decreasing tolerance is set in the KSP used
361: within the nonlinear iteration, so that each Newton iteration requests
362: better accuracy than the previous one. The constant correction tolerance
363: flag stops this behaviour.
365: Level: intermediate
367: .seealso: NEPRIIGetConstCorrectionTol()
368: @*/
369: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
370: {
376: PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
377: return(0);
378: }
380: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
381: {
382: NEP_RII *ctx = (NEP_RII*)nep->data;
385: *cct = ctx->cctol;
386: return(0);
387: }
389: /*@
390: NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.
392: Not Collective
394: Input Parameter:
395: . nep - nonlinear eigenvalue solver
397: Output Parameter:
398: . cct - the value of the constant tolerance flag
400: Level: intermediate
402: .seealso: NEPRIISetConstCorrectionTol()
403: @*/
404: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
405: {
411: PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
412: return(0);
413: }
415: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
416: {
418: NEP_RII *ctx = (NEP_RII*)nep->data;
421: PetscObjectReference((PetscObject)ksp);
422: KSPDestroy(&ctx->ksp);
423: ctx->ksp = ksp;
424: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
425: nep->state = NEP_STATE_INITIAL;
426: return(0);
427: }
429: /*@
430: NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
431: eigenvalue solver.
433: Collective on NEP
435: Input Parameters:
436: + nep - eigenvalue solver
437: - ksp - the linear solver object
439: Level: advanced
441: .seealso: NEPRIIGetKSP()
442: @*/
443: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
444: {
451: PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
452: return(0);
453: }
455: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
456: {
458: NEP_RII *ctx = (NEP_RII*)nep->data;
461: if (!ctx->ksp) {
462: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
463: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
464: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
465: KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_");
466: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
467: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
468: KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
469: }
470: *ksp = ctx->ksp;
471: return(0);
472: }
474: /*@
475: NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
476: the nonlinear eigenvalue solver.
478: Not Collective
480: Input Parameter:
481: . nep - nonlinear eigenvalue solver
483: Output Parameter:
484: . ksp - the linear solver object
486: Level: advanced
488: .seealso: NEPRIISetKSP()
489: @*/
490: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
491: {
497: PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
498: return(0);
499: }
501: PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
502: {
504: NEP_RII *ctx = (NEP_RII*)nep->data;
505: PetscBool isascii;
508: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
509: if (isascii) {
510: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
511: PetscViewerASCIIPrintf(viewer," maximum number of inner iterations: %D\n",ctx->max_inner_it);
512: if (ctx->cctol) {
513: PetscViewerASCIIPrintf(viewer," using a constant tolerance for the linear solver\n");
514: }
515: if (ctx->lag) {
516: PetscViewerASCIIPrintf(viewer," updating the preconditioner every %D iterations\n",ctx->lag);
517: }
518: KSPView(ctx->ksp,viewer);
519: }
520: return(0);
521: }
523: PetscErrorCode NEPReset_RII(NEP nep)
524: {
526: NEP_RII *ctx = (NEP_RII*)nep->data;
529: KSPReset(ctx->ksp);
530: return(0);
531: }
533: PetscErrorCode NEPDestroy_RII(NEP nep)
534: {
536: NEP_RII *ctx = (NEP_RII*)nep->data;
539: KSPDestroy(&ctx->ksp);
540: PetscFree(nep->data);
541: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL);
542: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL);
543: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL);
544: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL);
545: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL);
546: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL);
547: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL);
548: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL);
549: return(0);
550: }
552: PETSC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
553: {
555: NEP_RII *ctx;
558: PetscNewLog(nep,&ctx);
559: nep->data = (void*)ctx;
560: ctx->max_inner_it = 10;
561: ctx->lag = 1;
562: ctx->cctol = PETSC_FALSE;
564: nep->ops->solve = NEPSolve_RII;
565: nep->ops->setup = NEPSetUp_RII;
566: nep->ops->setfromoptions = NEPSetFromOptions_RII;
567: nep->ops->reset = NEPReset_RII;
568: nep->ops->destroy = NEPDestroy_RII;
569: nep->ops->view = NEPView_RII;
571: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII);
572: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII);
573: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII);
574: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII);
575: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII);
576: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII);
577: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII);
578: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII);
579: return(0);
580: }