Actual source code: rii.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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*/
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: KSP ksp; /* linear solver object */
33: } NEP_RII;
35: PetscErrorCode NEPSetUp_RII(NEP nep)
36: {
38: PetscBool istrivial;
41: if (nep->ncv) { PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"); }
42: nep->ncv = nep->nev;
43: if (nep->mpd) { PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"); }
44: nep->mpd = nep->nev;
45: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
46: if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");
48: RGIsTrivial(nep->rg,&istrivial);
49: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");
51: NEPAllocateSolution(nep,0);
52: NEPSetWorkVecs(nep,2);
53: return(0);
54: }
56: PetscErrorCode NEPSolve_RII(NEP nep)
57: {
58: PetscErrorCode ierr;
59: NEP_RII *ctx = (NEP_RII*)nep->data;
60: Mat T,Tp,H;
61: Vec uu,u,r,delta;
62: PetscScalar lambda,lambda2,sigma,a1,a2,corr,*Hp,*Ap;
63: PetscReal nrm,resnorm=1.0,ktol=0.1,perr;
64: PetscBool skip=PETSC_FALSE;
65: PetscInt inner_its,its=0,ldh,ldds,i,j;
66: NEP_EXT_OP extop=NULL;
67: KSPConvergedReason kspreason;
70: /* get initial approximation of eigenvalue and eigenvector */
71: NEPGetDefaultShift(nep,&sigma);
72: lambda = sigma;
73: if (!nep->nini) {
74: BVSetRandomColumn(nep->V,0);
75: BVNormColumn(nep->V,0,NORM_2,&nrm);
76: BVScaleColumn(nep->V,0,1.0/nrm);
77: }
78: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
79: NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
80: NEPDeflationCreateVec(extop,&u);
81: VecDuplicate(u,&r);
82: VecDuplicate(u,&delta);
83: BVGetColumn(nep->V,0,&uu);
84: NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE);
85: BVRestoreColumn(nep->V,0,&uu);
87: /* prepare linear solver */
88: NEPDeflationSolveSetUp(extop,sigma);
90: VecCopy(u,r);
91: NEPDeflationFunctionSolve(extop,r,u);
92: VecNorm(u,NORM_2,&nrm);
93: VecScale(u,1.0/nrm);
95: /* Restart loop */
96: while (nep->reason == NEP_CONVERGED_ITERATING) {
97: its++;
99: /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
100: estimate as starting value. */
101: inner_its=0;
102: lambda2 = lambda;
103: do {
104: NEPDeflationComputeFunction(extop,lambda,&T);
105: MatMult(T,u,r);
106: VecDot(r,u,&a1);
107: if (inner_its && PetscAbsScalar(a1)/PetscAbsScalar(lambda)<=PETSC_SQRT_MACHINE_EPSILON) break;
108: NEPDeflationComputeJacobian(extop,lambda,&Tp);
109: MatMult(Tp,u,r);
110: VecDot(r,u,&a2);
111: corr = a1/a2;
112: lambda = lambda - corr;
113: inner_its++;
114: } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);
116: /* form residual, r = T(lambda)*u */
117: NEPDeflationComputeFunction(extop,lambda,&T);
118: MatMult(T,u,r);
120: /* convergence test */
121: perr = nep->errest[nep->nconv];
122: VecNorm(r,NORM_2,&resnorm);
123: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
124: nep->eigr[nep->nconv] = lambda;
125: if (nep->errest[nep->nconv]<=nep->tol) {
126: nep->nconv = nep->nconv + 1;
127: NEPDeflationLocking(extop,u,lambda);
128: nep->its += its;
129: skip = PETSC_TRUE;
130: }
131: (*nep->stopping)(nep,nep->its+its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
132: if (!skip || nep->reason>0) {
133: NEPMonitor(nep,nep->its+its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
134: }
136: if (nep->reason == NEP_CONVERGED_ITERATING) {
137: if (!skip) {
138: /* update preconditioner and set adaptive tolerance */
139: if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) {
140: NEPDeflationSolveSetUp(extop,lambda2);
141: }
142: if (!ctx->cctol) {
143: ktol = PetscMax(ktol/2.0,PETSC_MACHINE_EPSILON*10.0);
144: KSPSetTolerances(ctx->ksp,ktol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
145: }
147: /* eigenvector correction: delta = T(sigma)\r */
148: NEPDeflationFunctionSolve(extop,r,delta);
149: KSPGetConvergedReason(ctx->ksp,&kspreason);
150: if (kspreason<0) {
151: PetscInfo1(nep,"iter=%D, linear solve failed, stopping solve\n",nep->its);
152: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
153: break;
154: }
156: /* update eigenvector: u = u - delta */
157: VecAXPY(u,-1.0,delta);
159: /* normalize eigenvector */
160: VecNormalize(u,NULL);
161: } else {
162: its = -1;
163: NEPDeflationSetRandomVec(extop,u);
164: NEPDeflationSolveSetUp(extop,sigma);
165: VecCopy(u,r);
166: NEPDeflationFunctionSolve(extop,r,u);
167: VecNorm(u,NORM_2,&nrm);
168: VecScale(u,1.0/nrm);
169: lambda = sigma;
170: skip = PETSC_FALSE;
171: }
172: }
173: }
174: NEPDeflationGetInvariantPair(extop,NULL,&H);
175: MatGetSize(H,NULL,&ldh);
176: DSSetType(nep->ds,DSNHEP);
177: DSAllocate(nep->ds,PetscMax(nep->nconv,1));
178: DSGetLeadingDimension(nep->ds,&ldds);
179: MatDenseGetArray(H,&Hp);
180: DSGetArray(nep->ds,DS_MAT_A,&Ap);
181: for (j=0;j<nep->nconv;j++)
182: for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
183: DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
184: MatDenseRestoreArray(H,&Hp);
185: MatDestroy(&H);
186: DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
187: DSSolve(nep->ds,nep->eigr,nep->eigi);
188: NEPDeflationReset(extop);
189: VecDestroy(&u);
190: VecDestroy(&r);
191: VecDestroy(&delta);
192: return(0);
193: }
195: PetscErrorCode NEPSetFromOptions_RII(PetscOptionItems *PetscOptionsObject,NEP nep)
196: {
198: NEP_RII *ctx = (NEP_RII*)nep->data;
199: PetscBool flg;
200: PetscInt i;
203: PetscOptionsHead(PetscOptionsObject,"NEP RII Options");
205: PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&ctx->max_inner_it,NULL);
207: PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL);
209: i = 0;
210: PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg);
211: if (flg) { NEPRIISetLagPreconditioner(nep,i); }
213: PetscOptionsTail();
215: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
216: KSPSetFromOptions(ctx->ksp);
217: return(0);
218: }
220: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
221: {
222: NEP_RII *ctx = (NEP_RII*)nep->data;
225: ctx->max_inner_it = its;
226: return(0);
227: }
229: /*@
230: NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
231: used in the RII solver. These are the Newton iterations related to the computation
232: of the nonlinear Rayleigh functional.
234: Logically Collective on NEP
236: Input Parameters:
237: + nep - nonlinear eigenvalue solver
238: - its - maximum inner iterations
240: Level: advanced
242: .seealso: NEPRIIGetMaximumIterations()
243: @*/
244: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
245: {
251: PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
252: return(0);
253: }
255: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
256: {
257: NEP_RII *ctx = (NEP_RII*)nep->data;
260: *its = ctx->max_inner_it;
261: return(0);
262: }
264: /*@
265: NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.
267: Not Collective
269: Input Parameter:
270: . nep - nonlinear eigenvalue solver
272: Output Parameter:
273: . its - maximum inner iterations
275: Level: advanced
277: .seealso: NEPRIISetMaximumIterations()
278: @*/
279: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
280: {
286: PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
287: return(0);
288: }
290: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
291: {
292: NEP_RII *ctx = (NEP_RII*)nep->data;
295: if (lag<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
296: ctx->lag = lag;
297: return(0);
298: }
300: /*@
301: NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
302: nonlinear solve.
304: Logically Collective on NEP
306: Input Parameters:
307: + nep - nonlinear eigenvalue solver
308: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
309: computed within the nonlinear iteration, 2 means every second time
310: the Jacobian is built, etc.
312: Options Database Keys:
313: . -nep_rii_lag_preconditioner <lag>
315: Notes:
316: The default is 1.
317: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
319: Level: intermediate
321: .seealso: NEPRIIGetLagPreconditioner()
322: @*/
323: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
324: {
330: PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
331: return(0);
332: }
334: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
335: {
336: NEP_RII *ctx = (NEP_RII*)nep->data;
339: *lag = ctx->lag;
340: return(0);
341: }
343: /*@
344: NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
346: Not Collective
348: Input Parameter:
349: . nep - nonlinear eigenvalue solver
351: Output Parameter:
352: . lag - the lag parameter
354: Level: intermediate
356: .seealso: NEPRIISetLagPreconditioner()
357: @*/
358: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
359: {
365: PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
366: return(0);
367: }
369: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
370: {
371: NEP_RII *ctx = (NEP_RII*)nep->data;
374: ctx->cctol = cct;
375: return(0);
376: }
378: /*@
379: NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
380: in the linear solver constant.
382: Logically Collective on NEP
384: Input Parameters:
385: + nep - nonlinear eigenvalue solver
386: - cct - a boolean value
388: Options Database Keys:
389: . -nep_rii_const_correction_tol <bool> - set the boolean flag
391: Notes:
392: By default, an exponentially decreasing tolerance is set in the KSP used
393: within the nonlinear iteration, so that each Newton iteration requests
394: better accuracy than the previous one. The constant correction tolerance
395: flag stops this behaviour.
397: Level: intermediate
399: .seealso: NEPRIIGetConstCorrectionTol()
400: @*/
401: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
402: {
408: PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
409: return(0);
410: }
412: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
413: {
414: NEP_RII *ctx = (NEP_RII*)nep->data;
417: *cct = ctx->cctol;
418: return(0);
419: }
421: /*@
422: NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.
424: Not Collective
426: Input Parameter:
427: . nep - nonlinear eigenvalue solver
429: Output Parameter:
430: . cct - the value of the constant tolerance flag
432: Level: intermediate
434: .seealso: NEPRIISetConstCorrectionTol()
435: @*/
436: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
437: {
443: PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
444: return(0);
445: }
447: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
448: {
450: NEP_RII *ctx = (NEP_RII*)nep->data;
453: PetscObjectReference((PetscObject)ksp);
454: KSPDestroy(&ctx->ksp);
455: ctx->ksp = ksp;
456: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
457: nep->state = NEP_STATE_INITIAL;
458: return(0);
459: }
461: /*@
462: NEPRIISetKSP - Associate a linear solver object (KSP) to the nonlinear
463: eigenvalue solver.
465: Collective on NEP
467: Input Parameters:
468: + nep - eigenvalue solver
469: - ksp - the linear solver object
471: Level: advanced
473: .seealso: NEPRIIGetKSP()
474: @*/
475: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
476: {
483: PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
484: return(0);
485: }
487: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
488: {
490: NEP_RII *ctx = (NEP_RII*)nep->data;
493: if (!ctx->ksp) {
494: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
495: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
496: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
497: KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_");
498: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
499: PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
500: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
501: KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
502: }
503: *ksp = ctx->ksp;
504: return(0);
505: }
507: /*@
508: NEPRIIGetKSP - Retrieve the linear solver object (KSP) associated with
509: the nonlinear eigenvalue solver.
511: Not Collective
513: Input Parameter:
514: . nep - nonlinear eigenvalue solver
516: Output Parameter:
517: . ksp - the linear solver object
519: Level: advanced
521: .seealso: NEPRIISetKSP()
522: @*/
523: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
524: {
530: PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
531: return(0);
532: }
534: PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
535: {
537: NEP_RII *ctx = (NEP_RII*)nep->data;
538: PetscBool isascii;
541: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
542: if (isascii) {
543: if (!ctx->ksp) { NEPRIIGetKSP(nep,&ctx->ksp); }
544: PetscViewerASCIIPrintf(viewer," maximum number of inner iterations: %D\n",ctx->max_inner_it);
545: if (ctx->cctol) {
546: PetscViewerASCIIPrintf(viewer," using a constant tolerance for the linear solver\n");
547: }
548: if (ctx->lag) {
549: PetscViewerASCIIPrintf(viewer," updating the preconditioner every %D iterations\n",ctx->lag);
550: }
551: PetscViewerASCIIPushTab(viewer);
552: KSPView(ctx->ksp,viewer);
553: PetscViewerASCIIPopTab(viewer);
554: }
555: return(0);
556: }
558: PetscErrorCode NEPReset_RII(NEP nep)
559: {
561: NEP_RII *ctx = (NEP_RII*)nep->data;
564: KSPReset(ctx->ksp);
565: return(0);
566: }
568: PetscErrorCode NEPDestroy_RII(NEP nep)
569: {
571: NEP_RII *ctx = (NEP_RII*)nep->data;
574: KSPDestroy(&ctx->ksp);
575: PetscFree(nep->data);
576: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL);
577: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL);
578: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL);
579: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL);
580: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL);
581: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL);
582: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL);
583: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL);
584: return(0);
585: }
587: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
588: {
590: NEP_RII *ctx;
593: PetscNewLog(nep,&ctx);
594: nep->data = (void*)ctx;
595: ctx->max_inner_it = 10;
596: ctx->lag = 1;
597: ctx->cctol = PETSC_FALSE;
599: nep->useds = PETSC_TRUE;
601: nep->ops->solve = NEPSolve_RII;
602: nep->ops->setup = NEPSetUp_RII;
603: nep->ops->setfromoptions = NEPSetFromOptions_RII;
604: nep->ops->reset = NEPReset_RII;
605: nep->ops->destroy = NEPDestroy_RII;
606: nep->ops->view = NEPView_RII;
607: nep->ops->computevectors = NEPComputeVectors_Schur;
609: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII);
610: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII);
611: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII);
612: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII);
613: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII);
614: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII);
615: PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII);
616: PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII);
617: return(0);
618: }