Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/stimpl.h> /*I "slepcst.h" I*/
15 : #include "filter.h"
16 :
17 : /*
18 : Operator (filter):
19 : Op P M
20 : if nmat=1: p(A) NULL p(A)
21 : */
22 8 : static PetscErrorCode STComputeOperator_Filter(ST st)
23 : {
24 8 : ST_FILTER *ctx = (ST_FILTER*)st->data;
25 :
26 8 : PetscFunctionBegin;
27 8 : PetscCheck(st->nmat==1,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Only implemented for standard eigenvalue problem");
28 8 : PetscCheck(ctx->intb<PETSC_MAX_REAL || ctx->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an interval with STFilterSetInterval()");
29 8 : PetscCheck(ctx->right!=0.0 || ctx->left!=0.0,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an approximate numerical range with STFilterSetRange()");
30 8 : 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 8 : if (!ctx->polyDegree) ctx->polyDegree = 100;
32 8 : ctx->frame[0] = ctx->left;
33 8 : ctx->frame[1] = ctx->inta;
34 8 : ctx->frame[2] = ctx->intb;
35 8 : ctx->frame[3] = ctx->right;
36 8 : PetscCall(STFilter_FILTLAN_setFilter(st,&st->T[0]));
37 8 : st->M = st->T[0];
38 8 : PetscCall(MatDestroy(&st->P));
39 8 : PetscFunctionReturn(PETSC_SUCCESS);
40 : }
41 :
42 8 : static PetscErrorCode STSetUp_Filter(ST st)
43 : {
44 8 : PetscFunctionBegin;
45 8 : PetscCall(STSetWorkVecs(st,4));
46 8 : PetscFunctionReturn(PETSC_SUCCESS);
47 : }
48 :
49 11 : static PetscErrorCode STSetFromOptions_Filter(ST st,PetscOptionItems *PetscOptionsObject)
50 : {
51 11 : PetscReal array[2]={0,0};
52 11 : PetscInt k;
53 11 : PetscBool flg;
54 :
55 11 : PetscFunctionBegin;
56 11 : PetscOptionsHeadBegin(PetscOptionsObject,"ST Filter Options");
57 :
58 11 : k = 2;
59 11 : PetscCall(PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",array,&k,&flg));
60 11 : if (flg) {
61 0 : PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_interval (comma-separated without spaces)");
62 0 : PetscCall(STFilterSetInterval(st,array[0],array[1]));
63 : }
64 11 : k = 2;
65 11 : PetscCall(PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg));
66 11 : if (flg) {
67 2 : PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_range (comma-separated without spaces)");
68 2 : PetscCall(STFilterSetRange(st,array[0],array[1]));
69 : }
70 11 : PetscCall(PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg));
71 11 : if (flg) PetscCall(STFilterSetDegree(st,k));
72 :
73 11 : PetscOptionsHeadEnd();
74 11 : PetscFunctionReturn(PETSC_SUCCESS);
75 : }
76 :
77 9 : static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
78 : {
79 9 : ST_FILTER *ctx = (ST_FILTER*)st->data;
80 :
81 9 : PetscFunctionBegin;
82 9 : PetscCheck(inta<intb,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
83 9 : if (ctx->inta != inta || ctx->intb != intb) {
84 6 : ctx->inta = inta;
85 6 : ctx->intb = intb;
86 6 : st->state = ST_STATE_INITIAL;
87 6 : st->opready = PETSC_FALSE;
88 6 : ctx->filtch = PETSC_TRUE;
89 : }
90 9 : PetscFunctionReturn(PETSC_SUCCESS);
91 : }
92 :
93 : /*@
94 : STFilterSetInterval - Defines the interval containing the desired eigenvalues.
95 :
96 : Logically Collective
97 :
98 : Input Parameters:
99 : + st - the spectral transformation context
100 : . inta - left end of the interval
101 : - intb - right end of the interval
102 :
103 : Options Database Key:
104 : . -st_filter_interval <a,b> - set [a,b] as the interval of interest
105 :
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.
110 :
111 : Common usage is to set the interval in EPS with EPSSetInterval().
112 :
113 : The interval must be contained within the numerical range of the matrix, see
114 : STFilterSetRange().
115 :
116 : Level: intermediate
117 :
118 : .seealso: STFilterGetInterval(), STFilterSetRange(), EPSSetInterval()
119 : @*/
120 9 : PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
121 : {
122 9 : PetscFunctionBegin;
123 9 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
124 27 : PetscValidLogicalCollectiveReal(st,inta,2);
125 27 : PetscValidLogicalCollectiveReal(st,intb,3);
126 9 : PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
127 9 : PetscFunctionReturn(PETSC_SUCCESS);
128 : }
129 :
130 3 : static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
131 : {
132 3 : ST_FILTER *ctx = (ST_FILTER*)st->data;
133 :
134 3 : PetscFunctionBegin;
135 3 : if (inta) *inta = ctx->inta;
136 3 : if (intb) *intb = ctx->intb;
137 3 : PetscFunctionReturn(PETSC_SUCCESS);
138 : }
139 :
140 : /*@
141 : STFilterGetInterval - Gets the interval containing the desired eigenvalues.
142 :
143 : Not Collective
144 :
145 : Input Parameter:
146 : . st - the spectral transformation context
147 :
148 : Output Parameters:
149 : + inta - left end of the interval
150 : - intb - right end of the interval
151 :
152 : Level: intermediate
153 :
154 : .seealso: STFilterSetInterval()
155 : @*/
156 3 : PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
157 : {
158 3 : PetscFunctionBegin;
159 3 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
160 3 : PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
161 3 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 9 : static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
165 : {
166 9 : ST_FILTER *ctx = (ST_FILTER*)st->data;
167 :
168 9 : PetscFunctionBegin;
169 9 : PetscCheck(left<right,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be left<right");
170 9 : if (ctx->left != left || ctx->right != right) {
171 7 : ctx->left = left;
172 7 : ctx->right = right;
173 7 : st->state = ST_STATE_INITIAL;
174 7 : st->opready = PETSC_FALSE;
175 7 : ctx->filtch = PETSC_TRUE;
176 : }
177 9 : PetscFunctionReturn(PETSC_SUCCESS);
178 : }
179 :
180 : /*@
181 : STFilterSetRange - Defines the numerical range (or field of values) of the matrix, that is,
182 : the interval containing all eigenvalues.
183 :
184 : Logically Collective
185 :
186 : Input Parameters:
187 : + st - the spectral transformation context
188 : . left - left end of the interval
189 : - right - right end of the interval
190 :
191 : Options Database Key:
192 : . -st_filter_range <a,b> - set [a,b] as the numerical range
193 :
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.
197 :
198 : Level: intermediate
199 :
200 : .seealso: STFilterGetRange(), STFilterSetInterval()
201 : @*/
202 9 : PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
203 : {
204 9 : PetscFunctionBegin;
205 9 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
206 27 : PetscValidLogicalCollectiveReal(st,left,2);
207 27 : PetscValidLogicalCollectiveReal(st,right,3);
208 9 : PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
209 9 : PetscFunctionReturn(PETSC_SUCCESS);
210 : }
211 :
212 10 : static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
213 : {
214 10 : ST_FILTER *ctx = (ST_FILTER*)st->data;
215 :
216 10 : PetscFunctionBegin;
217 10 : if (left) *left = ctx->left;
218 10 : if (right) *right = ctx->right;
219 10 : PetscFunctionReturn(PETSC_SUCCESS);
220 : }
221 :
222 : /*@
223 : STFilterGetRange - Gets the interval containing all eigenvalues.
224 :
225 : Not Collective
226 :
227 : Input Parameter:
228 : . st - the spectral transformation context
229 :
230 : Output Parameters:
231 : + left - left end of the interval
232 : - right - right end of the interval
233 :
234 : Level: intermediate
235 :
236 : .seealso: STFilterSetRange()
237 : @*/
238 10 : PetscErrorCode STFilterGetRange(ST st,PetscReal *left,PetscReal *right)
239 : {
240 10 : PetscFunctionBegin;
241 10 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
242 10 : PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
243 10 : PetscFunctionReturn(PETSC_SUCCESS);
244 : }
245 :
246 5 : static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
247 : {
248 5 : ST_FILTER *ctx = (ST_FILTER*)st->data;
249 :
250 5 : PetscFunctionBegin;
251 5 : if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
252 0 : ctx->polyDegree = 0;
253 0 : st->state = ST_STATE_INITIAL;
254 0 : st->opready = PETSC_FALSE;
255 0 : ctx->filtch = PETSC_TRUE;
256 : } else {
257 5 : PetscCheck(deg>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
258 5 : if (ctx->polyDegree != deg) {
259 3 : ctx->polyDegree = deg;
260 3 : st->state = ST_STATE_INITIAL;
261 3 : st->opready = PETSC_FALSE;
262 3 : ctx->filtch = PETSC_TRUE;
263 : }
264 : }
265 5 : PetscFunctionReturn(PETSC_SUCCESS);
266 : }
267 :
268 : /*@
269 : STFilterSetDegree - Sets the degree of the filter polynomial.
270 :
271 : Logically Collective
272 :
273 : Input Parameters:
274 : + st - the spectral transformation context
275 : - deg - polynomial degree
276 :
277 : Options Database Key:
278 : . -st_filter_degree <deg> - sets the degree of the filter polynomial
279 :
280 : Level: intermediate
281 :
282 : .seealso: STFilterGetDegree()
283 : @*/
284 5 : PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
285 : {
286 5 : PetscFunctionBegin;
287 5 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
288 15 : PetscValidLogicalCollectiveInt(st,deg,2);
289 5 : PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
290 5 : PetscFunctionReturn(PETSC_SUCCESS);
291 : }
292 :
293 3 : static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
294 : {
295 3 : ST_FILTER *ctx = (ST_FILTER*)st->data;
296 :
297 3 : PetscFunctionBegin;
298 3 : *deg = ctx->polyDegree;
299 3 : PetscFunctionReturn(PETSC_SUCCESS);
300 : }
301 :
302 : /*@
303 : STFilterGetDegree - Gets the degree of the filter polynomial.
304 :
305 : Not Collective
306 :
307 : Input Parameter:
308 : . st - the spectral transformation context
309 :
310 : Output Parameter:
311 : . deg - polynomial degree
312 :
313 : Level: intermediate
314 :
315 : .seealso: STFilterSetDegree()
316 : @*/
317 3 : PetscErrorCode STFilterGetDegree(ST st,PetscInt *deg)
318 : {
319 3 : PetscFunctionBegin;
320 3 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
321 3 : PetscAssertPointer(deg,2);
322 3 : PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
323 3 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 23 : static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
327 : {
328 23 : ST_FILTER *ctx = (ST_FILTER*)st->data;
329 :
330 23 : PetscFunctionBegin;
331 23 : *gamma = ctx->filterInfo->yLimit;
332 23 : PetscFunctionReturn(PETSC_SUCCESS);
333 : }
334 :
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.
338 :
339 : Not Collective
340 :
341 : Input Parameter:
342 : . st - the spectral transformation context
343 :
344 : Output Parameter:
345 : . gamma - the threshold value
346 :
347 : Level: developer
348 :
349 : .seealso: STFilterGetRange()
350 : @*/
351 23 : PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
352 : {
353 23 : PetscFunctionBegin;
354 23 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
355 23 : PetscAssertPointer(gamma,2);
356 23 : PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
357 23 : PetscFunctionReturn(PETSC_SUCCESS);
358 : }
359 :
360 12 : static PetscErrorCode STReset_Filter(ST st)
361 : {
362 12 : ST_FILTER *ctx = (ST_FILTER*)st->data;
363 :
364 12 : PetscFunctionBegin;
365 12 : ctx->left = 0.0;
366 12 : ctx->right = 0.0;
367 12 : PetscCall(MatDestroy(&ctx->T));
368 12 : PetscCall(MatDestroyMatrices(ctx->nW,&ctx->W));
369 12 : ctx->nW = 0;
370 12 : PetscFunctionReturn(PETSC_SUCCESS);
371 : }
372 :
373 0 : static PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
374 : {
375 0 : ST_FILTER *ctx = (ST_FILTER*)st->data;
376 0 : PetscBool isascii;
377 :
378 0 : PetscFunctionBegin;
379 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
380 0 : if (isascii) {
381 0 : PetscCall(PetscViewerASCIIPrintf(viewer," interval of desired eigenvalues: [%g,%g]\n",(double)ctx->inta,(double)ctx->intb));
382 0 : PetscCall(PetscViewerASCIIPrintf(viewer," numerical range: [%g,%g]\n",(double)ctx->left,(double)ctx->right));
383 0 : PetscCall(PetscViewerASCIIPrintf(viewer," degree of filter polynomial: %" PetscInt_FMT "\n",ctx->polyDegree));
384 0 : if (st->state>=ST_STATE_SETUP) PetscCall(PetscViewerASCIIPrintf(viewer," limit to accept eigenvalues: theta=%g\n",(double)ctx->filterInfo->yLimit));
385 : }
386 0 : PetscFunctionReturn(PETSC_SUCCESS);
387 : }
388 :
389 6 : static PetscErrorCode STDestroy_Filter(ST st)
390 : {
391 6 : ST_FILTER *ctx = (ST_FILTER*)st->data;
392 :
393 6 : PetscFunctionBegin;
394 6 : PetscCall(PetscFree(ctx->opts));
395 6 : PetscCall(PetscFree(ctx->filterInfo));
396 6 : PetscCall(PetscFree(ctx->baseFilter));
397 6 : PetscCall(PetscFree(st->data));
398 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL));
399 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL));
400 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL));
401 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL));
402 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL));
403 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL));
404 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL));
405 6 : PetscFunctionReturn(PETSC_SUCCESS);
406 : }
407 :
408 6 : SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
409 : {
410 6 : ST_FILTER *ctx;
411 6 : FILTLAN_IOP iop;
412 6 : FILTLAN_PFI pfi;
413 :
414 6 : PetscFunctionBegin;
415 6 : PetscCall(PetscNew(&ctx));
416 6 : st->data = (void*)ctx;
417 :
418 6 : st->usesksp = PETSC_FALSE;
419 :
420 6 : ctx->inta = PETSC_MIN_REAL;
421 6 : ctx->intb = PETSC_MAX_REAL;
422 6 : ctx->left = 0.0;
423 6 : ctx->right = 0.0;
424 6 : ctx->polyDegree = 0;
425 6 : ctx->baseDegree = 10;
426 :
427 6 : PetscCall(PetscNew(&iop));
428 6 : ctx->opts = iop;
429 6 : iop->intervalWeights[0] = 100.0;
430 6 : iop->intervalWeights[1] = 1.0;
431 6 : iop->intervalWeights[2] = 1.0;
432 6 : iop->intervalWeights[3] = 1.0;
433 6 : iop->intervalWeights[4] = 100.0;
434 6 : iop->transIntervalRatio = 0.6;
435 6 : iop->reverseInterval = PETSC_FALSE;
436 6 : iop->initialPlateau = 0.1;
437 6 : iop->plateauShrinkRate = 1.5;
438 6 : iop->initialShiftStep = 0.01;
439 6 : iop->shiftStepExpanRate = 1.5;
440 6 : iop->maxInnerIter = 30;
441 6 : iop->yLimitTol = 0.0001;
442 6 : iop->maxOuterIter = 50;
443 6 : iop->numGridPoints = 200;
444 6 : iop->yBottomLine = 0.001;
445 6 : iop->yRippleLimit = 0.8;
446 :
447 6 : PetscCall(PetscNew(&pfi));
448 6 : ctx->filterInfo = pfi;
449 :
450 6 : st->ops->apply = STApply_Generic;
451 6 : st->ops->setup = STSetUp_Filter;
452 6 : st->ops->computeoperator = STComputeOperator_Filter;
453 6 : st->ops->setfromoptions = STSetFromOptions_Filter;
454 6 : st->ops->destroy = STDestroy_Filter;
455 6 : st->ops->reset = STReset_Filter;
456 6 : st->ops->view = STView_Filter;
457 :
458 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter));
459 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter));
460 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter));
461 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter));
462 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter));
463 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter));
464 6 : PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter));
465 6 : PetscFunctionReturn(PETSC_SUCCESS);
466 : }
|
