Actual source code: filter.c
slepc-3.21.0 2024-03-30
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: Filter spectral transformation, to encapsulate polynomial filters
12: */
14: #include <slepc/private/stimpl.h>
15: #include "filter.h"
17: /*
18: Operator (filter):
19: Op P M
20: if nmat=1: p(A) NULL p(A)
21: */
22: static PetscErrorCode STComputeOperator_Filter(ST st)
23: {
24: ST_FILTER *ctx = (ST_FILTER*)st->data;
26: PetscFunctionBegin;
27: PetscCheck(st->nmat==1,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Only implemented for standard eigenvalue problem");
28: PetscCheck(ctx->intb<PETSC_MAX_REAL || ctx->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an interval with STFilterSetInterval()");
29: PetscCheck(ctx->right!=0.0 || ctx->left!=0.0,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an approximate numerical range with STFilterSetRange()");
30: PetscCheck(ctx->left<=ctx->inta && ctx->right>=ctx->intb,PetscObjectComm((PetscObject)st),PETSC_ERR_USER_INPUT,"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);
31: if (!ctx->polyDegree) ctx->polyDegree = 100;
32: ctx->frame[0] = ctx->left;
33: ctx->frame[1] = ctx->inta;
34: ctx->frame[2] = ctx->intb;
35: ctx->frame[3] = ctx->right;
36: PetscCall(STFilter_FILTLAN_setFilter(st,&st->T[0]));
37: st->M = st->T[0];
38: PetscCall(MatDestroy(&st->P));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode STSetUp_Filter(ST st)
43: {
44: PetscFunctionBegin;
45: PetscCall(STSetWorkVecs(st,4));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode STSetFromOptions_Filter(ST st,PetscOptionItems *PetscOptionsObject)
50: {
51: PetscReal array[2]={0,0};
52: PetscInt k;
53: PetscBool flg;
55: PetscFunctionBegin;
56: PetscOptionsHeadBegin(PetscOptionsObject,"ST Filter Options");
58: k = 2;
59: PetscCall(PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",array,&k,&flg));
60: if (flg) {
61: PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_interval (comma-separated without spaces)");
62: PetscCall(STFilterSetInterval(st,array[0],array[1]));
63: }
64: k = 2;
65: PetscCall(PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg));
66: if (flg) {
67: PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_range (comma-separated without spaces)");
68: PetscCall(STFilterSetRange(st,array[0],array[1]));
69: }
70: PetscCall(PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg));
71: if (flg) PetscCall(STFilterSetDegree(st,k));
73: PetscOptionsHeadEnd();
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
78: {
79: ST_FILTER *ctx = (ST_FILTER*)st->data;
81: PetscFunctionBegin;
82: PetscCheck(inta<intb,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
83: if (ctx->inta != inta || ctx->intb != intb) {
84: ctx->inta = inta;
85: ctx->intb = intb;
86: st->state = ST_STATE_INITIAL;
87: st->opready = PETSC_FALSE;
88: ctx->filtch = PETSC_TRUE;
89: }
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*@
94: STFilterSetInterval - Defines the interval containing the desired eigenvalues.
96: Logically Collective
98: Input Parameters:
99: + st - the spectral transformation context
100: . inta - left end of the interval
101: - intb - right end of the interval
103: Options Database Key:
104: . -st_filter_interval <a,b> - set [a,b] as the interval of interest
106: Notes:
107: The filter will be configured to emphasize eigenvalues contained in the given
108: interval, and damp out eigenvalues outside it. If the interval is open, then
109: the filter is low- or high-pass, otherwise it is mid-pass.
111: Common usage is to set the interval in EPS with EPSSetInterval().
113: The interval must be contained within the numerical range of the matrix, see
114: STFilterSetRange().
116: Level: intermediate
118: .seealso: STFilterGetInterval(), STFilterSetRange(), EPSSetInterval()
119: @*/
120: PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
121: {
122: PetscFunctionBegin;
126: PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
131: {
132: ST_FILTER *ctx = (ST_FILTER*)st->data;
134: PetscFunctionBegin;
135: if (inta) *inta = ctx->inta;
136: if (intb) *intb = ctx->intb;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@
141: STFilterGetInterval - Gets the interval containing the desired eigenvalues.
143: Not Collective
145: Input Parameter:
146: . st - the spectral transformation context
148: Output Parameters:
149: + inta - left end of the interval
150: - intb - right end of the interval
152: Level: intermediate
154: .seealso: STFilterSetInterval()
155: @*/
156: PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
157: {
158: PetscFunctionBegin;
160: PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
165: {
166: ST_FILTER *ctx = (ST_FILTER*)st->data;
168: PetscFunctionBegin;
169: PetscCheck(left<right,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: st->opready = PETSC_FALSE;
175: ctx->filtch = PETSC_TRUE;
176: }
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: STFilterSetRange - Defines the numerical range (or field of values) of the matrix, that is,
182: the interval containing all eigenvalues.
184: Logically Collective
186: Input Parameters:
187: + st - the spectral transformation context
188: . left - left end of the interval
189: - right - right end of the interval
191: Options Database Key:
192: . -st_filter_range <a,b> - set [a,b] as the numerical range
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: Level: intermediate
200: .seealso: STFilterGetRange(), STFilterSetInterval()
201: @*/
202: PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
203: {
204: PetscFunctionBegin;
208: PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
213: {
214: ST_FILTER *ctx = (ST_FILTER*)st->data;
216: PetscFunctionBegin;
217: if (left) *left = ctx->left;
218: if (right) *right = ctx->right;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: STFilterGetRange - Gets the interval containing all eigenvalues.
225: Not Collective
227: Input Parameter:
228: . st - the spectral transformation context
230: Output Parameters:
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: {
240: PetscFunctionBegin;
242: PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
247: {
248: ST_FILTER *ctx = (ST_FILTER*)st->data;
250: PetscFunctionBegin;
251: if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
252: ctx->polyDegree = 0;
253: st->state = ST_STATE_INITIAL;
254: st->opready = PETSC_FALSE;
255: ctx->filtch = PETSC_TRUE;
256: } else {
257: PetscCheck(deg>0,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: st->opready = PETSC_FALSE;
262: ctx->filtch = PETSC_TRUE;
263: }
264: }
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: STFilterSetDegree - Sets the degree of the filter polynomial.
271: Logically Collective
273: Input Parameters:
274: + st - the spectral transformation context
275: - deg - polynomial degree
277: Options Database Key:
278: . -st_filter_degree <deg> - sets the degree of the filter polynomial
280: Level: intermediate
282: .seealso: STFilterGetDegree()
283: @*/
284: PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
285: {
286: PetscFunctionBegin;
289: PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
294: {
295: ST_FILTER *ctx = (ST_FILTER*)st->data;
297: PetscFunctionBegin;
298: *deg = ctx->polyDegree;
299: PetscFunctionReturn(PETSC_SUCCESS);
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: {
319: PetscFunctionBegin;
321: PetscAssertPointer(deg,2);
322: PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
327: {
328: ST_FILTER *ctx = (ST_FILTER*)st->data;
330: PetscFunctionBegin;
331: *gamma = ctx->filterInfo->yLimit;
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: /*@
336: STFilterGetThreshold - Gets the threshold value gamma such that rho(lambda)>=gamma for
337: eigenvalues lambda inside the wanted interval and rho(lambda)<gamma for those outside.
339: Not Collective
341: Input Parameter:
342: . st - the spectral transformation context
344: Output Parameter:
345: . gamma - the threshold value
347: Level: developer
349: .seealso: STFilterGetRange()
350: @*/
351: PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
352: {
353: PetscFunctionBegin;
355: PetscAssertPointer(gamma,2);
356: PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode STReset_Filter(ST st)
361: {
362: ST_FILTER *ctx = (ST_FILTER*)st->data;
364: PetscFunctionBegin;
365: ctx->left = 0.0;
366: ctx->right = 0.0;
367: PetscCall(MatDestroy(&ctx->T));
368: PetscCall(MatDestroyMatrices(ctx->nW,&ctx->W));
369: ctx->nW = 0;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
374: {
375: ST_FILTER *ctx = (ST_FILTER*)st->data;
376: PetscBool isascii;
378: PetscFunctionBegin;
379: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
380: if (isascii) {
381: PetscCall(PetscViewerASCIIPrintf(viewer," interval of desired eigenvalues: [%g,%g]\n",(double)ctx->inta,(double)ctx->intb));
382: PetscCall(PetscViewerASCIIPrintf(viewer," numerical range: [%g,%g]\n",(double)ctx->left,(double)ctx->right));
383: PetscCall(PetscViewerASCIIPrintf(viewer," degree of filter polynomial: %" PetscInt_FMT "\n",ctx->polyDegree));
384: if (st->state>=ST_STATE_SETUP) PetscCall(PetscViewerASCIIPrintf(viewer," limit to accept eigenvalues: theta=%g\n",(double)ctx->filterInfo->yLimit));
385: }
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode STDestroy_Filter(ST st)
390: {
391: ST_FILTER *ctx = (ST_FILTER*)st->data;
393: PetscFunctionBegin;
394: PetscCall(PetscFree(ctx->opts));
395: PetscCall(PetscFree(ctx->filterInfo));
396: PetscCall(PetscFree(ctx->baseFilter));
397: PetscCall(PetscFree(st->data));
398: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL));
399: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL));
400: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL));
401: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL));
402: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL));
403: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL));
404: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
409: {
410: ST_FILTER *ctx;
411: FILTLAN_IOP iop;
412: FILTLAN_PFI pfi;
414: PetscFunctionBegin;
415: PetscCall(PetscNew(&ctx));
416: st->data = (void*)ctx;
418: st->usesksp = PETSC_FALSE;
420: ctx->inta = PETSC_MIN_REAL;
421: ctx->intb = PETSC_MAX_REAL;
422: ctx->left = 0.0;
423: ctx->right = 0.0;
424: ctx->polyDegree = 0;
425: ctx->baseDegree = 10;
427: PetscCall(PetscNew(&iop));
428: ctx->opts = iop;
429: iop->intervalWeights[0] = 100.0;
430: iop->intervalWeights[1] = 1.0;
431: iop->intervalWeights[2] = 1.0;
432: iop->intervalWeights[3] = 1.0;
433: iop->intervalWeights[4] = 100.0;
434: iop->transIntervalRatio = 0.6;
435: iop->reverseInterval = PETSC_FALSE;
436: iop->initialPlateau = 0.1;
437: iop->plateauShrinkRate = 1.5;
438: iop->initialShiftStep = 0.01;
439: iop->shiftStepExpanRate = 1.5;
440: iop->maxInnerIter = 30;
441: iop->yLimitTol = 0.0001;
442: iop->maxOuterIter = 50;
443: iop->numGridPoints = 200;
444: iop->yBottomLine = 0.001;
445: iop->yRippleLimit = 0.8;
447: PetscCall(PetscNew(&pfi));
448: ctx->filterInfo = pfi;
450: st->ops->apply = STApply_Generic;
451: st->ops->setup = STSetUp_Filter;
452: st->ops->computeoperator = STComputeOperator_Filter;
453: st->ops->setfromoptions = STSetFromOptions_Filter;
454: st->ops->destroy = STDestroy_Filter;
455: st->ops->reset = STReset_Filter;
456: st->ops->view = STView_Filter;
458: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter));
459: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter));
460: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter));
461: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter));
462: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter));
463: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter));
464: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }