Actual source code: filter.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: Filter spectral transformation, to encapsulate polynomial filters
12: */
14: #include <slepc/private/stimpl.h> /*I "slepcst.h" I*/
15: #include "./filter.h"
17: PetscErrorCode STApply_Filter(ST st,Vec x,Vec y)
18: {
22: STFilter_FILTLAN_Apply(st,x,y);
23: return(0);
24: }
26: PetscErrorCode STSetUp_Filter(ST st)
27: {
29: ST_FILTER *ctx = (ST_FILTER*)st->data;
32: STSetWorkVecs(st,4);
33: if (st->nmat>1) SETERRQ(PetscObjectComm((PetscObject)st),1,"Only implemented for standard eigenvalue problem");
34: if (ctx->intb >= PETSC_MAX_REAL && ctx->inta <= PETSC_MIN_REAL) SETERRQ(PetscObjectComm((PetscObject)st),1,"Must pass an interval with STFilterSetInterval()");
35: if (ctx->right == 0.0 && ctx->left == 0.0) SETERRQ(PetscObjectComm((PetscObject)st),1,"Must pass an approximate numerical range with STFilterSetRange()");
36: if (ctx->left > ctx->inta || ctx->right < ctx->intb) SETERRQ4(PetscObjectComm((PetscObject)st),1,"The requested interval [%g,%g] must be contained in the numerical range [%g,%g]",(double)ctx->inta,(double)ctx->intb,(double)ctx->left,(double)ctx->right);
37: if (!ctx->polyDegree) ctx->polyDegree = 100;
38: ctx->frame[0] = ctx->left;
39: ctx->frame[1] = ctx->inta;
40: ctx->frame[2] = ctx->intb;
41: ctx->frame[3] = ctx->right;
42: STFilter_FILTLAN_setFilter(st);
43: return(0);
44: }
46: PetscErrorCode STSetFromOptions_Filter(PetscOptionItems *PetscOptionsObject,ST st)
47: {
49: PetscReal array[2]={0,0};
50: PetscInt k;
51: PetscBool flg;
54: PetscOptionsHead(PetscOptionsObject,"ST Filter Options");
56: k = 2;
57: PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",array,&k,&flg);
58: if (flg) {
59: if (k<2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_interval (comma-separated without spaces)");
60: STFilterSetInterval(st,array[0],array[1]);
61: }
62: k = 2;
63: PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg);
64: if (flg) {
65: if (k<2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_range (comma-separated without spaces)");
66: STFilterSetRange(st,array[0],array[1]);
67: }
68: PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg);
69: if (flg) { STFilterSetDegree(st,k); }
71: PetscOptionsTail();
72: return(0);
73: }
75: static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
76: {
77: ST_FILTER *ctx = (ST_FILTER*)st->data;
80: if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
81: if (ctx->inta != inta || ctx->intb != intb) {
82: ctx->inta = inta;
83: ctx->intb = intb;
84: st->state = ST_STATE_INITIAL;
85: }
86: return(0);
87: }
89: /*@
90: STFilterSetInterval - Defines the interval containing the desired eigenvalues.
92: Logically Collective on ST
94: Input Parameters:
95: + st - the spectral transformation context
96: . inta - left end of the interval
97: - intb - right end of the interval
99: Options Database Key:
100: . -st_filter_interval <a,b> - set [a,b] as the interval of interest
102: Level: intermediate
104: Notes:
105: The filter will be configured to emphasize eigenvalues contained in the given
106: interval, and damp out eigenvalues outside it. If the interval is open, then
107: the filter is low- or high-pass, otherwise it is mid-pass.
109: Common usage is to set the interval in EPS with EPSSetInterval().
111: The interval must be contained within the numerical range of the matrix, see
112: STFilterSetRange().
114: .seealso: STFilterGetInterval(), STFilterSetRange(), EPSSetInterval()
115: @*/
116: PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
117: {
124: PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
125: return(0);
126: }
128: static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
129: {
130: ST_FILTER *ctx = (ST_FILTER*)st->data;
133: if (inta) *inta = ctx->inta;
134: if (intb) *intb = ctx->intb;
135: return(0);
136: }
138: /*@
139: STFilterGetInterval - Gets the interval containing the desired eigenvalues.
141: Not Collective
143: Input Parameter:
144: . st - the spectral transformation context
146: Output Parameter:
147: + inta - left end of the interval
148: - intb - right end of the interval
150: Level: intermediate
152: .seealso: STFilterSetInterval()
153: @*/
154: PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
155: {
160: PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
161: return(0);
162: }
164: static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
165: {
166: ST_FILTER *ctx = (ST_FILTER*)st->data;
169: if (left>right) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be left<right");
170: if (ctx->left != left || ctx->right != right) {
171: ctx->left = left;
172: ctx->right = right;
173: st->state = ST_STATE_INITIAL;
174: }
175: return(0);
176: }
178: /*@
179: STFilterSetRange - Defines the numerical range (or field of values) of the matrix, that is,
180: the interval containing all eigenvalues.
182: Logically Collective on ST
184: Input Parameters:
185: + st - the spectral transformation context
186: . left - left end of the interval
187: - right - right end of the interval
189: Options Database Key:
190: . -st_filter_range <a,b> - set [a,b] as the numerical range
192: Level: intermediate
194: Notes:
195: The filter will be most effective if the numerical range is tight, that is, left and right
196: are good approximations to the leftmost and rightmost eigenvalues, respectively.
198: .seealso: STFilterGetRange(), STFilterSetInterval()
199: @*/
200: PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
201: {
208: PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
209: return(0);
210: }
212: static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
213: {
214: ST_FILTER *ctx = (ST_FILTER*)st->data;
217: if (left) *left = ctx->left;
218: if (right) *right = ctx->right;
219: return(0);
220: }
222: /*@
223: STFilterGetRange - Gets the interval containing all eigenvalues.
225: Not Collective
227: Input Parameter:
228: . st - the spectral transformation context
230: Output Parameter:
231: + left - left end of the interval
232: - right - right end of the interval
234: Level: intermediate
236: .seealso: STFilterSetRange()
237: @*/
238: PetscErrorCode STFilterGetRange(ST st,PetscReal *left,PetscReal *right)
239: {
244: PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
245: return(0);
246: }
248: static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
249: {
250: ST_FILTER *ctx = (ST_FILTER*)st->data;
253: if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
254: ctx->polyDegree = 0;
255: st->state = ST_STATE_INITIAL;
256: } else {
257: if (deg<=0) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
258: if (ctx->polyDegree != deg) {
259: ctx->polyDegree = deg;
260: st->state = ST_STATE_INITIAL;
261: }
262: }
263: return(0);
264: }
266: /*@
267: STFilterSetDegree - Sets the degree of the filter polynomial.
269: Logically Collective on ST
271: Input Parameters:
272: + st - the spectral transformation context
273: - deg - polynomial degree
275: Options Database Key:
276: . -st_filter_degree <deg> - sets the degree of the filter polynomial
278: Level: intermediate
280: .seealso: STFilterGetDegree()
281: @*/
282: PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
283: {
289: PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
290: return(0);
291: }
293: static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
294: {
295: ST_FILTER *ctx = (ST_FILTER*)st->data;
298: *deg = ctx->polyDegree;
299: return(0);
300: }
302: /*@
303: STFilterGetDegree - Gets the degree of the filter polynomial.
305: Not Collective
307: Input Parameter:
308: . st - the spectral transformation context
310: Output Parameter:
311: . deg - polynomial degree
313: Level: intermediate
315: .seealso: STFilterSetDegree()
316: @*/
317: PetscErrorCode STFilterGetDegree(ST st,PetscInt *deg)
318: {
324: PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
325: return(0);
326: }
328: static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
329: {
330: ST_FILTER *ctx = (ST_FILTER*)st->data;
333: *gamma = ctx->filterInfo->yLimit;
334: return(0);
335: }
337: /*@
338: STFilterGetThreshold - Gets the threshold value gamma such that rho(lambda)>=gamma for
339: eigenvalues lambda inside the wanted interval and rho(lambda)<gamma for those outside.
341: Not Collective
343: Input Parameter:
344: . st - the spectral transformation context
346: Output Parameter:
347: . gamma - the threshold value
349: Level: developer
350: @*/
351: PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
352: {
358: PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
359: return(0);
360: }
362: PetscErrorCode STReset_Filter(ST st)
363: {
364: ST_FILTER *ctx = (ST_FILTER*)st->data;
367: ctx->left = 0.0;
368: ctx->right = 0.0;
369: return(0);
370: }
372: PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
373: {
375: ST_FILTER *ctx = (ST_FILTER*)st->data;
376: PetscBool isascii;
379: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
380: if (isascii) {
381: PetscViewerASCIIPrintf(viewer," Filter: interval of desired eigenvalues is [%g,%g]\n",(double)ctx->inta,(double)ctx->intb);
382: PetscViewerASCIIPrintf(viewer," Filter: numerical range is [%g,%g]\n",(double)ctx->left,(double)ctx->right);
383: PetscViewerASCIIPrintf(viewer," Filter: degree of filter polynomial is %D\n",ctx->polyDegree);
384: if (st->state>=ST_STATE_SETUP) {
385: PetscViewerASCIIPrintf(viewer," Filter: limit to accept eigenvalues: theta=%g\n",ctx->filterInfo->yLimit);
386: }
387: }
388: return(0);
389: }
391: PetscErrorCode STDestroy_Filter(ST st)
392: {
394: ST_FILTER *ctx = (ST_FILTER*)st->data;
397: PetscFree(ctx->opts);
398: PetscFree(ctx->filterInfo);
399: PetscFree(ctx->baseFilter);
400: PetscFree(st->data);
401: PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL);
402: PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL);
403: PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL);
404: PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL);
405: PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL);
406: PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL);
407: PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL);
408: return(0);
409: }
411: SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
412: {
414: ST_FILTER *ctx;
415: FILTLAN_IOP iop;
416: FILTLAN_PFI pfi;
418: PetscNewLog(st,&ctx);
419: st->data = (void*)ctx;
421: ctx->inta = PETSC_MIN_REAL;
422: ctx->intb = PETSC_MAX_REAL;
423: ctx->left = 0.0;
424: ctx->right = 0.0;
425: ctx->polyDegree = 0;
426: ctx->baseDegree = 10;
428: PetscNewLog(st,&iop);
429: ctx->opts = iop;
430: iop->intervalWeights[0] = 100.0;
431: iop->intervalWeights[1] = 1.0;
432: iop->intervalWeights[2] = 1.0;
433: iop->intervalWeights[3] = 1.0;
434: iop->intervalWeights[4] = 100.0;
435: iop->transIntervalRatio = 0.6;
436: iop->reverseInterval = PETSC_FALSE;
437: iop->initialPlateau = 0.1;
438: iop->plateauShrinkRate = 1.5;
439: iop->initialShiftStep = 0.01;
440: iop->shiftStepExpanRate = 1.5;
441: iop->maxInnerIter = 30;
442: iop->yLimitTol = 0.0001;
443: iop->maxOuterIter = 50;
444: iop->numGridPoints = 200;
445: iop->yBottomLine = 0.001;
446: iop->yRippleLimit = 0.8;
448: PetscNewLog(st,&pfi);
449: ctx->filterInfo = pfi;
451: st->ops->apply = STApply_Filter;
452: st->ops->setfromoptions = STSetFromOptions_Filter;
453: st->ops->setup = STSetUp_Filter;
454: st->ops->destroy = STDestroy_Filter;
455: st->ops->reset = STReset_Filter;
456: st->ops->view = STView_Filter;
458: PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter);
459: PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter);
460: PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter);
461: PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter);
462: PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter);
463: PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter);
464: PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter);
465: return(0);
466: }