Actual source code: interpol.c
slepc-3.18.0 2022-10-01
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: "interpol"
13: Method: Polynomial interpolation
15: Algorithm:
17: Uses a PEP object to solve the interpolated NEP. Currently supports
18: only Chebyshev interpolation on an interval.
20: References:
22: [1] C. Effenberger and D. Kresser, "Chebyshev interpolation for
23: nonlinear eigenvalue problems", BIT 52:933-951, 2012.
24: */
26: #include <slepc/private/nepimpl.h>
27: #include <slepc/private/pepimpl.h>
29: typedef struct {
30: PEP pep;
31: PetscReal tol; /* tolerance for norm of polynomial coefficients */
32: PetscInt maxdeg; /* maximum degree of interpolation polynomial */
33: PetscInt deg; /* actual degree of interpolation polynomial */
34: } NEP_INTERPOL;
36: PetscErrorCode NEPSetUp_Interpol(NEP nep)
37: {
38: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
39: ST st;
40: RG rg;
41: PetscReal a,b,c,d,s,tol;
42: PetscScalar zero=0.0;
43: PetscBool flg,istrivial,trackall;
44: PetscInt its,in;
46: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
48: if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
49: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
51: NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);
53: /* transfer PEP options */
54: if (!ctx->pep) NEPInterpolGetPEP(nep,&ctx->pep);
55: PEPSetBasis(ctx->pep,PEP_BASIS_CHEBYSHEV1);
56: PEPSetWhichEigenpairs(ctx->pep,PEP_TARGET_MAGNITUDE);
57: PetscObjectTypeCompare((PetscObject)ctx->pep,PEPJD,&flg);
58: if (!flg) {
59: PEPGetST(ctx->pep,&st);
60: STSetType(st,STSINVERT);
61: }
62: PEPSetDimensions(ctx->pep,nep->nev,nep->ncv,nep->mpd);
63: PEPGetTolerances(ctx->pep,&tol,&its);
64: if (tol==PETSC_DEFAULT) tol = SlepcDefaultTol(nep->tol);
65: if (ctx->tol==PETSC_DEFAULT) ctx->tol = tol;
66: if (its==PETSC_DEFAULT) its = nep->max_it;
67: PEPSetTolerances(ctx->pep,tol,its);
68: NEPGetTrackAll(nep,&trackall);
69: PEPSetTrackAll(ctx->pep,trackall);
71: /* transfer region options */
72: RGIsTrivial(nep->rg,&istrivial);
74: PetscObjectTypeCompare((PetscObject)nep->rg,RGINTERVAL,&flg);
76: RGIntervalGetEndpoints(nep->rg,&a,&b,&c,&d);
78: PEPGetRG(ctx->pep,&rg);
79: RGSetType(rg,RGINTERVAL);
81: s = 2.0/(b-a);
82: c = c*s;
83: d = d*s;
84: RGIntervalSetEndpoints(rg,-1.0,1.0,c,d);
85: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
87: PEPSetTarget(ctx->pep,(nep->target-(a+b)/2)*s);
89: NEPAllocateSolution(nep,0);
90: return 0;
91: }
93: /*
94: Input:
95: d, number of nodes to compute
96: a,b, interval extremes
97: Output:
98: *x, array containing the d Chebyshev nodes of the interval [a,b]
99: *dct2, coefficients to compute a discrete cosine transformation (DCT-II)
100: */
101: static PetscErrorCode ChebyshevNodes(PetscInt d,PetscReal a,PetscReal b,PetscScalar *x,PetscReal *dct2)
102: {
103: PetscInt j,i;
104: PetscReal t;
106: for (j=0;j<d+1;j++) {
107: t = ((2*j+1)*PETSC_PI)/(2*(d+1));
108: x[j] = (a+b)/2.0+((b-a)/2.0)*PetscCosReal(t);
109: for (i=0;i<d+1;i++) dct2[j*(d+1)+i] = PetscCosReal(i*t);
110: }
111: return 0;
112: }
114: PetscErrorCode NEPSolve_Interpol(NEP nep)
115: {
116: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
117: Mat *A,*P;
118: PetscScalar *x,*fx,t;
119: PetscReal *cs,a,b,s,aprox,aprox0=1.0,*matnorm;
120: PetscInt i,j,k,deg=ctx->maxdeg;
121: PetscBool hasmnorm=PETSC_FALSE;
122: Vec vr,vi=NULL;
123: ST st;
125: PetscMalloc4((deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx,nep->nt,&matnorm);
126: for (j=0;j<nep->nt;j++) {
127: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
128: if (!hasmnorm) break;
129: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
130: }
131: if (!hasmnorm) for (j=0;j<nep->nt;j++) matnorm[j] = 1.0;
132: RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL);
133: ChebyshevNodes(deg,a,b,x,cs);
134: for (j=0;j<nep->nt;j++) {
135: for (i=0;i<=deg;i++) FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]);
136: }
137: /* Polynomial coefficients */
138: PetscMalloc1(deg+1,&A);
139: if (nep->P) PetscMalloc1(deg+1,&P);
140: ctx->deg = deg;
141: for (k=0;k<=deg;k++) {
142: MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]);
143: if (nep->P) MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P[k]);
144: t = 0.0;
145: for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
146: t *= 2.0/(deg+1);
147: if (k==0) t /= 2.0;
148: aprox = matnorm[0]*PetscAbsScalar(t);
149: MatScale(A[k],t);
150: if (nep->P) MatScale(P[k],t);
151: for (j=1;j<nep->nt;j++) {
152: t = 0.0;
153: for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
154: t *= 2.0/(deg+1);
155: if (k==0) t /= 2.0;
156: aprox += matnorm[j]*PetscAbsScalar(t);
157: MatAXPY(A[k],t,nep->A[j],nep->mstr);
158: if (nep->P) MatAXPY(P[k],t,nep->P[j],nep->mstrp);
159: }
160: if (k==0) aprox0 = aprox;
161: if (k>1 && aprox/aprox0<ctx->tol) { ctx->deg = k; deg = k; break; }
162: }
163: PEPSetOperators(ctx->pep,deg+1,A);
164: MatDestroyMatrices(deg+1,&A);
165: if (nep->P) {
166: PEPGetST(ctx->pep,&st);
167: STSetSplitPreconditioner(st,deg+1,P,nep->mstrp);
168: MatDestroyMatrices(deg+1,&P);
169: }
170: PetscFree4(cs,x,fx,matnorm);
172: /* Solve polynomial eigenproblem */
173: PEPSolve(ctx->pep);
174: PEPGetConverged(ctx->pep,&nep->nconv);
175: PEPGetIterationNumber(ctx->pep,&nep->its);
176: PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason);
177: BVSetActiveColumns(nep->V,0,nep->nconv);
178: BVCreateVec(nep->V,&vr);
179: #if !defined(PETSC_USE_COMPLEX)
180: VecDuplicate(vr,&vi);
181: #endif
182: s = 2.0/(b-a);
183: for (i=0;i<nep->nconv;i++) {
184: PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],vr,vi);
185: nep->eigr[i] /= s;
186: nep->eigr[i] += (a+b)/2.0;
187: nep->eigi[i] /= s;
188: BVInsertVec(nep->V,i,vr);
189: #if !defined(PETSC_USE_COMPLEX)
190: if (nep->eigi[i]!=0.0) BVInsertVec(nep->V,++i,vi);
191: #endif
192: }
193: VecDestroy(&vr);
194: VecDestroy(&vi);
196: nep->state = NEP_STATE_EIGENVECTORS;
197: return 0;
198: }
200: static PetscErrorCode PEPMonitor_Interpol(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
201: {
202: PetscInt i,n;
203: NEP nep = (NEP)ctx;
204: PetscReal a,b,s;
205: ST st;
207: n = PetscMin(nest,nep->ncv);
208: for (i=0;i<n;i++) {
209: nep->eigr[i] = eigr[i];
210: nep->eigi[i] = eigi[i];
211: nep->errest[i] = errest[i];
212: }
213: PEPGetST(pep,&st);
214: STBackTransform(st,n,nep->eigr,nep->eigi);
215: RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL);
216: s = 2.0/(b-a);
217: for (i=0;i<n;i++) {
218: nep->eigr[i] /= s;
219: nep->eigr[i] += (a+b)/2.0;
220: nep->eigi[i] /= s;
221: }
222: NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest);
223: return 0;
224: }
226: PetscErrorCode NEPSetFromOptions_Interpol(NEP nep,PetscOptionItems *PetscOptionsObject)
227: {
228: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
229: PetscInt i;
230: PetscBool flg1,flg2;
231: PetscReal r;
233: PetscOptionsHeadBegin(PetscOptionsObject,"NEP Interpol Options");
235: NEPInterpolGetInterpolation(nep,&r,&i);
236: if (!i) i = PETSC_DEFAULT;
237: PetscOptionsInt("-nep_interpol_interpolation_degree","Maximum degree of polynomial interpolation","NEPInterpolSetInterpolation",i,&i,&flg1);
238: PetscOptionsReal("-nep_interpol_interpolation_tol","Tolerance for interpolation coefficients","NEPInterpolSetInterpolation",r,&r,&flg2);
239: if (flg1 || flg2) NEPInterpolSetInterpolation(nep,r,i);
241: PetscOptionsHeadEnd();
243: if (!ctx->pep) NEPInterpolGetPEP(nep,&ctx->pep);
244: PEPSetFromOptions(ctx->pep);
245: return 0;
246: }
248: static PetscErrorCode NEPInterpolSetInterpolation_Interpol(NEP nep,PetscReal tol,PetscInt degree)
249: {
250: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
252: if (tol == PETSC_DEFAULT) {
253: ctx->tol = PETSC_DEFAULT;
254: nep->state = NEP_STATE_INITIAL;
255: } else {
257: ctx->tol = tol;
258: }
259: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
260: ctx->maxdeg = 0;
261: if (nep->state) NEPReset(nep);
262: nep->state = NEP_STATE_INITIAL;
263: } else {
265: if (ctx->maxdeg != degree) {
266: ctx->maxdeg = degree;
267: if (nep->state) NEPReset(nep);
268: nep->state = NEP_STATE_INITIAL;
269: }
270: }
271: return 0;
272: }
274: /*@
275: NEPInterpolSetInterpolation - Sets the tolerance and maximum degree when building
276: the interpolation polynomial.
278: Collective on nep
280: Input Parameters:
281: + nep - nonlinear eigenvalue solver
282: . tol - tolerance to stop computing polynomial coefficients
283: - deg - maximum degree of interpolation
285: Options Database Key:
286: + -nep_interpol_interpolation_tol <tol> - Sets the tolerance to stop computing polynomial coefficients
287: - -nep_interpol_interpolation_degree <degree> - Sets the maximum degree of interpolation
289: Notes:
290: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
292: Level: advanced
294: .seealso: NEPInterpolGetInterpolation()
295: @*/
296: PetscErrorCode NEPInterpolSetInterpolation(NEP nep,PetscReal tol,PetscInt deg)
297: {
301: PetscTryMethod(nep,"NEPInterpolSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,deg));
302: return 0;
303: }
305: static PetscErrorCode NEPInterpolGetInterpolation_Interpol(NEP nep,PetscReal *tol,PetscInt *deg)
306: {
307: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
309: if (tol) *tol = ctx->tol;
310: if (deg) *deg = ctx->maxdeg;
311: return 0;
312: }
314: /*@
315: NEPInterpolGetInterpolation - Gets the tolerance and maximum degree when building
316: the interpolation polynomial.
318: Not Collective
320: Input Parameter:
321: . nep - nonlinear eigenvalue solver
323: Output Parameters:
324: + tol - tolerance to stop computing polynomial coefficients
325: - deg - maximum degree of interpolation
327: Level: advanced
329: .seealso: NEPInterpolSetInterpolation()
330: @*/
331: PetscErrorCode NEPInterpolGetInterpolation(NEP nep,PetscReal *tol,PetscInt *deg)
332: {
334: PetscUseMethod(nep,"NEPInterpolGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,deg));
335: return 0;
336: }
338: static PetscErrorCode NEPInterpolSetPEP_Interpol(NEP nep,PEP pep)
339: {
340: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
342: PetscObjectReference((PetscObject)pep);
343: PEPDestroy(&ctx->pep);
344: ctx->pep = pep;
345: nep->state = NEP_STATE_INITIAL;
346: return 0;
347: }
349: /*@
350: NEPInterpolSetPEP - Associate a polynomial eigensolver object (PEP) to the
351: nonlinear eigenvalue solver.
353: Collective on nep
355: Input Parameters:
356: + nep - nonlinear eigenvalue solver
357: - pep - the polynomial eigensolver object
359: Level: advanced
361: .seealso: NEPInterpolGetPEP()
362: @*/
363: PetscErrorCode NEPInterpolSetPEP(NEP nep,PEP pep)
364: {
368: PetscTryMethod(nep,"NEPInterpolSetPEP_C",(NEP,PEP),(nep,pep));
369: return 0;
370: }
372: static PetscErrorCode NEPInterpolGetPEP_Interpol(NEP nep,PEP *pep)
373: {
374: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
376: if (!ctx->pep) {
377: PEPCreate(PetscObjectComm((PetscObject)nep),&ctx->pep);
378: PetscObjectIncrementTabLevel((PetscObject)ctx->pep,(PetscObject)nep,1);
379: PEPSetOptionsPrefix(ctx->pep,((PetscObject)nep)->prefix);
380: PEPAppendOptionsPrefix(ctx->pep,"nep_interpol_");
381: PetscObjectSetOptions((PetscObject)ctx->pep,((PetscObject)nep)->options);
382: PEPMonitorSet(ctx->pep,PEPMonitor_Interpol,nep,NULL);
383: }
384: *pep = ctx->pep;
385: return 0;
386: }
388: /*@
389: NEPInterpolGetPEP - Retrieve the polynomial eigensolver object (PEP)
390: associated with the nonlinear eigenvalue solver.
392: Not Collective
394: Input Parameter:
395: . nep - nonlinear eigenvalue solver
397: Output Parameter:
398: . pep - the polynomial eigensolver object
400: Level: advanced
402: .seealso: NEPInterpolSetPEP()
403: @*/
404: PetscErrorCode NEPInterpolGetPEP(NEP nep,PEP *pep)
405: {
408: PetscUseMethod(nep,"NEPInterpolGetPEP_C",(NEP,PEP*),(nep,pep));
409: return 0;
410: }
412: PetscErrorCode NEPView_Interpol(NEP nep,PetscViewer viewer)
413: {
414: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
415: PetscBool isascii;
417: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
418: if (isascii) {
419: if (!ctx->pep) NEPInterpolGetPEP(nep,&ctx->pep);
420: PetscViewerASCIIPrintf(viewer," polynomial degree %" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->deg,ctx->maxdeg);
421: PetscViewerASCIIPrintf(viewer," tolerance for norm of polynomial coefficients %g\n",(double)ctx->tol);
422: PetscViewerASCIIPushTab(viewer);
423: PEPView(ctx->pep,viewer);
424: PetscViewerASCIIPopTab(viewer);
425: }
426: return 0;
427: }
429: PetscErrorCode NEPReset_Interpol(NEP nep)
430: {
431: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
433: PEPReset(ctx->pep);
434: return 0;
435: }
437: PetscErrorCode NEPDestroy_Interpol(NEP nep)
438: {
439: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
441: PEPDestroy(&ctx->pep);
442: PetscFree(nep->data);
443: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NULL);
444: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NULL);
445: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NULL);
446: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NULL);
447: return 0;
448: }
450: SLEPC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
451: {
452: NEP_INTERPOL *ctx;
454: PetscNew(&ctx);
455: nep->data = (void*)ctx;
456: ctx->maxdeg = 5;
457: ctx->tol = PETSC_DEFAULT;
459: nep->ops->solve = NEPSolve_Interpol;
460: nep->ops->setup = NEPSetUp_Interpol;
461: nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
462: nep->ops->reset = NEPReset_Interpol;
463: nep->ops->destroy = NEPDestroy_Interpol;
464: nep->ops->view = NEPView_Interpol;
466: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NEPInterpolSetInterpolation_Interpol);
467: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NEPInterpolGetInterpolation_Interpol);
468: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol);
469: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol);
470: return 0;
471: }