Actual source code: interpol.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: "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> /*I "slepcnep.h" I*/
27: #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
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: {
39: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
40: ST st;
41: RG rg;
42: PetscReal a,b,c,d,s,tol;
43: PetscScalar zero=0.0;
44: PetscBool flg,istrivial,trackall;
45: PetscInt its,in;
48: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
49: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
50: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
51: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPINTERPOL only available for split operator");
52: if (nep->stopping!=NEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support user-defined stopping test");
54: /* transfer PEP options */
55: if (!ctx->pep) { NEPInterpolGetPEP(nep,&ctx->pep); }
56: PEPSetBasis(ctx->pep,PEP_BASIS_CHEBYSHEV1);
57: PEPSetWhichEigenpairs(ctx->pep,PEP_TARGET_MAGNITUDE);
58: PetscObjectTypeCompare((PetscObject)ctx->pep,PEPJD,&flg);
59: if (!flg) {
60: PEPGetST(ctx->pep,&st);
61: STSetType(st,STSINVERT);
62: }
63: PEPSetDimensions(ctx->pep,nep->nev,nep->ncv?nep->ncv:PETSC_DEFAULT,nep->mpd?nep->mpd:PETSC_DEFAULT);
64: PEPGetTolerances(ctx->pep,&tol,&its);
65: if (tol==PETSC_DEFAULT) tol = (nep->tol==PETSC_DEFAULT)?SLEPC_DEFAULT_TOL:nep->tol;
66: if (ctx->tol==PETSC_DEFAULT) ctx->tol = tol;
67: if (!its) its = nep->max_it?nep->max_it:PETSC_DEFAULT;
68: PEPSetTolerances(ctx->pep,tol,its);
69: NEPGetTrackAll(nep,&trackall);
70: PEPSetTrackAll(ctx->pep,trackall);
72: /* transfer region options */
73: RGIsTrivial(nep->rg,&istrivial);
74: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPINTERPOL requires a nontrivial region");
75: PetscObjectTypeCompare((PetscObject)nep->rg,RGINTERVAL,&flg);
76: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for interval regions");
77: RGIntervalGetEndpoints(nep->rg,&a,&b,&c,&d);
78: if (a<=-PETSC_MAX_REAL || b>=PETSC_MAX_REAL) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for bounded intervals");
79: PEPGetRG(ctx->pep,&rg);
80: RGSetType(rg,RGINTERVAL);
81: if (a==b) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for intervals on the real axis");
82: s = 2.0/(b-a);
83: c = c*s;
84: d = d*s;
85: RGIntervalSetEndpoints(rg,-1.0,1.0,c,d);
86: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
87: if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
88: PEPSetTarget(ctx->pep,(nep->target-(a+b)/2)*s);
90: NEPAllocateSolution(nep,0);
91: return(0);
92: }
94: /*
95: Input:
96: d, number of nodes to compute
97: a,b, interval extrems
98: Output:
99: *x, array containing the d Chebyshev nodes of the interval [a,b]
100: *dct2, coefficients to compute a discrete cosine transformation (DCT-II)
101: */
102: static PetscErrorCode ChebyshevNodes(PetscInt d,PetscReal a,PetscReal b,PetscScalar *x,PetscReal *dct2)
103: {
104: PetscInt j,i;
105: PetscReal t;
108: for (j=0;j<d+1;j++) {
109: t = ((2*j+1)*PETSC_PI)/(2*(d+1));
110: x[j] = (a+b)/2.0+((b-a)/2.0)*PetscCosReal(t);
111: for (i=0;i<d+1;i++) dct2[j*(d+1)+i] = PetscCosReal(i*t);
112: }
113: return(0);
114: }
116: PetscErrorCode NEPSolve_Interpol(NEP nep)
117: {
119: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
120: Mat *A;
121: PetscScalar *x,*fx,t;
122: PetscReal *cs,a,b,s,aprox,aprox0=1.0,*matnorm;
123: PetscInt i,j,k,deg=ctx->maxdeg;
124: PetscBool hasmnorm=PETSC_FALSE;
125: Vec vr,vi=NULL;
128: PetscMalloc5(deg+1,&A,(deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx,nep->nt,&matnorm);
129: for (j=0;j<nep->nt;j++) {
130: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
131: if (!hasmnorm) break;
132: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
133: }
134: if (!hasmnorm) for (j=0;j<nep->nt;j++) matnorm[j] = 1.0;
135: RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL);
136: ChebyshevNodes(deg,a,b,x,cs);
137: for (j=0;j<nep->nt;j++) {
138: for (i=0;i<=deg;i++) {
139: FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]);
140: }
141: }
142: /* Polynomial coefficients */
143: ctx->deg = deg;
144: for (k=0;k<=deg;k++) {
145: MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]);
146: t = 0.0;
147: for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
148: t *= 2.0/(deg+1);
149: if (k==0) t /= 2.0;
150: aprox = matnorm[0]*PetscAbsScalar(t);
151: MatScale(A[k],t);
152: for (j=1;j<nep->nt;j++) {
153: t = 0.0;
154: for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
155: t *= 2.0/(deg+1);
156: if (k==0) t /= 2.0;
157: aprox += matnorm[j]*PetscAbsScalar(t);
158: MatAXPY(A[k],t,nep->A[j],nep->mstr);
159: }
160: if (k==0) aprox0 = aprox;
161: if (aprox/aprox0<ctx->tol) { ctx->deg = k; deg = k; break; }
162: }
164: PEPSetOperators(ctx->pep,deg+1,A);
165: for (k=0;k<=deg;k++) {
166: MatDestroy(&A[k]);
167: }
168: PetscFree5(A,cs,x,fx,matnorm);
170: /* Solve polynomial eigenproblem */
171: PEPSolve(ctx->pep);
172: PEPGetConverged(ctx->pep,&nep->nconv);
173: PEPGetIterationNumber(ctx->pep,&nep->its);
174: PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason);
175: BVSetActiveColumns(nep->V,0,nep->nconv);
176: BVCreateVec(nep->V,&vr);
177: #if !defined(PETSC_USE_COMPLEX)
178: VecDuplicate(vr,&vi);
179: #endif
180: s = 2.0/(b-a);
181: for (i=0;i<nep->nconv;i++) {
182: PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],vr,vi);
183: nep->eigr[i] /= s;
184: nep->eigr[i] += (a+b)/2.0;
185: nep->eigi[i] /= s;
186: BVInsertVec(nep->V,i,vr);
187: #if !defined(PETSC_USE_COMPLEX)
188: if (nep->eigi[i]!=0.0) {
189: BVInsertVec(nep->V,++i,vi);
190: }
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;
209: n = PetscMin(nest,nep->ncv);
210: for (i=0;i<n;i++) {
211: nep->eigr[i] = eigr[i];
212: nep->eigi[i] = eigi[i];
213: nep->errest[i] = errest[i];
214: }
215: PEPGetST(pep,&st);
216: STBackTransform(st,n,nep->eigr,nep->eigi);
217: RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL);
218: s = 2.0/(b-a);
219: for (i=0;i<n;i++) {
220: nep->eigr[i] /= s;
221: nep->eigr[i] += (a+b)/2.0;
222: nep->eigi[i] /= s;
223: }
224: NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest);
225: return(0);
226: }
228: PetscErrorCode NEPSetFromOptions_Interpol(PetscOptionItems *PetscOptionsObject,NEP nep)
229: {
231: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
232: PetscInt i;
233: PetscBool flg1,flg2;
234: PetscReal r;
237: PetscOptionsHead(PetscOptionsObject,"NEP Interpol Options");
239: NEPInterpolGetInterpolation(nep,&r,&i);
240: if (!i) i = PETSC_DEFAULT;
241: PetscOptionsInt("-nep_interpol_interpolation_degree","Maximum degree of polynomial interpolation","NEPInterpolSetInterpolation",i,&i,&flg1);
242: PetscOptionsReal("-nep_interpol_interpolation_tol","Tolerance for interpolation coefficients","NEPInterpolSetInterpolation",r,&r,&flg2);
243: if (flg1 || flg2) { NEPInterpolSetInterpolation(nep,r,i); }
245: PetscOptionsTail();
247: if (!ctx->pep) { NEPInterpolGetPEP(nep,&ctx->pep); }
248: PEPSetFromOptions(ctx->pep);
249: return(0);
250: }
252: static PetscErrorCode NEPInterpolSetInterpolation_Interpol(NEP nep,PetscReal tol,PetscInt degree)
253: {
255: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
258: if (tol == PETSC_DEFAULT) {
259: ctx->tol = PETSC_DEFAULT;
260: nep->state = NEP_STATE_INITIAL;
261: } else {
262: if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
263: ctx->tol = tol;
264: }
265: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
266: ctx->maxdeg = 0;
267: if (nep->state) { NEPReset(nep); }
268: nep->state = NEP_STATE_INITIAL;
269: } else {
270: if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
271: if (ctx->maxdeg != degree) {
272: ctx->maxdeg = degree;
273: if (nep->state) { NEPReset(nep); }
274: nep->state = NEP_STATE_INITIAL;
275: }
276: }
277: return(0);
278: }
280: /*@
281: NEPInterpolSetInterpolation - Sets the tolerance and maximum degree when building
282: the interpolation polynomial.
284: Collective on NEP
286: Input Parameters:
287: + nep - nonlinear eigenvalue solver
288: . tol - tolerance to stop computing polynomial coefficients
289: - deg - maximum degree of interpolation
291: Options Database Key:
292: + -nep_interpol_interpolation_tol <tol> - Sets the tolerance to stop computing polynomial coefficients
293: - -nep_interpol_interpolation_degree <degree> - Sets the maximum degree of interpolation
295: Notes:
296: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
298: Level: advanced
300: .seealso: NEPInterpolGetInterpolation()
301: @*/
302: PetscErrorCode NEPInterpolSetInterpolation(NEP nep,PetscReal tol,PetscInt deg)
303: {
310: PetscTryMethod(nep,"NEPInterpolSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,deg));
311: return(0);
312: }
314: static PetscErrorCode NEPInterpolGetInterpolation_Interpol(NEP nep,PetscReal *tol,PetscInt *deg)
315: {
316: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
319: if (tol) *tol = ctx->tol;
320: if (deg) *deg = ctx->maxdeg;
321: return(0);
322: }
324: /*@
325: NEPInterpolGetInterpolation - Gets the tolerance and maximum degree when building
326: the interpolation polynomial.
328: Not Collective
330: Input Parameter:
331: . nep - nonlinear eigenvalue solver
333: Output Parameter:
334: + tol - tolerance to stop computing polynomial coefficients
335: - deg - maximum degree of interpolation
337: Level: advanced
339: .seealso: NEPInterpolSetInterpolation()
340: @*/
341: PetscErrorCode NEPInterpolGetInterpolation(NEP nep,PetscReal *tol,PetscInt *deg)
342: {
347: PetscUseMethod(nep,"NEPInterpolGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,deg));
348: return(0);
349: }
351: static PetscErrorCode NEPInterpolSetPEP_Interpol(NEP nep,PEP pep)
352: {
354: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
357: PetscObjectReference((PetscObject)pep);
358: PEPDestroy(&ctx->pep);
359: ctx->pep = pep;
360: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->pep);
361: nep->state = NEP_STATE_INITIAL;
362: return(0);
363: }
365: /*@
366: NEPInterpolSetPEP - Associate a polynomial eigensolver object (PEP) to the
367: nonlinear eigenvalue solver.
369: Collective on NEP
371: Input Parameters:
372: + nep - nonlinear eigenvalue solver
373: - pep - the polynomial eigensolver object
375: Level: advanced
377: .seealso: NEPInterpolGetPEP()
378: @*/
379: PetscErrorCode NEPInterpolSetPEP(NEP nep,PEP pep)
380: {
387: PetscTryMethod(nep,"NEPInterpolSetPEP_C",(NEP,PEP),(nep,pep));
388: return(0);
389: }
391: static PetscErrorCode NEPInterpolGetPEP_Interpol(NEP nep,PEP *pep)
392: {
394: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
397: if (!ctx->pep) {
398: PEPCreate(PetscObjectComm((PetscObject)nep),&ctx->pep);
399: PetscObjectIncrementTabLevel((PetscObject)ctx->pep,(PetscObject)nep,1);
400: PEPSetOptionsPrefix(ctx->pep,((PetscObject)nep)->prefix);
401: PEPAppendOptionsPrefix(ctx->pep,"nep_interpol_");
402: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->pep);
403: PetscObjectSetOptions((PetscObject)ctx->pep,((PetscObject)nep)->options);
404: PEPMonitorSet(ctx->pep,PEPMonitor_Interpol,nep,NULL);
405: }
406: *pep = ctx->pep;
407: return(0);
408: }
410: /*@
411: NEPInterpolGetPEP - Retrieve the polynomial eigensolver object (PEP)
412: associated with the nonlinear eigenvalue solver.
414: Not Collective
416: Input Parameter:
417: . nep - nonlinear eigenvalue solver
419: Output Parameter:
420: . pep - the polynomial eigensolver object
422: Level: advanced
424: .seealso: NEPInterpolSetPEP()
425: @*/
426: PetscErrorCode NEPInterpolGetPEP(NEP nep,PEP *pep)
427: {
433: PetscUseMethod(nep,"NEPInterpolGetPEP_C",(NEP,PEP*),(nep,pep));
434: return(0);
435: }
437: PetscErrorCode NEPView_Interpol(NEP nep,PetscViewer viewer)
438: {
440: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
441: PetscBool isascii;
444: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
445: if (isascii) {
446: if (!ctx->pep) { NEPInterpolGetPEP(nep,&ctx->pep); }
447: PetscViewerASCIIPrintf(viewer," polynomial degree %D, max=%D\n",ctx->deg,ctx->maxdeg);
448: PetscViewerASCIIPrintf(viewer," tolerance for norm of polynomial coefficients %g\n",ctx->tol);
449: PetscViewerASCIIPushTab(viewer);
450: PEPView(ctx->pep,viewer);
451: PetscViewerASCIIPopTab(viewer);
452: }
453: return(0);
454: }
456: PetscErrorCode NEPReset_Interpol(NEP nep)
457: {
459: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
462: PEPReset(ctx->pep);
463: return(0);
464: }
466: PetscErrorCode NEPDestroy_Interpol(NEP nep)
467: {
469: NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
472: PEPDestroy(&ctx->pep);
473: PetscFree(nep->data);
474: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NULL);
475: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NULL);
476: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NULL);
477: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NULL);
478: return(0);
479: }
481: SLEPC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
482: {
484: NEP_INTERPOL *ctx;
487: PetscNewLog(nep,&ctx);
488: nep->data = (void*)ctx;
489: ctx->maxdeg = 5;
490: ctx->tol = PETSC_DEFAULT;
492: nep->ops->solve = NEPSolve_Interpol;
493: nep->ops->setup = NEPSetUp_Interpol;
494: nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
495: nep->ops->reset = NEPReset_Interpol;
496: nep->ops->destroy = NEPDestroy_Interpol;
497: nep->ops->view = NEPView_Interpol;
499: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NEPInterpolSetInterpolation_Interpol);
500: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NEPInterpolGetInterpolation_Interpol);
501: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol);
502: PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol);
503: return(0);
504: }