Actual source code: filter.c

slepc-3.14.2 2021-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: PetscErrorCode STComputeOperator_Filter(ST st)
 23: {
 25:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

 28:   if (st->nmat>1) SETERRQ(PetscObjectComm((PetscObject)st),1,"Only implemented for standard eigenvalue problem");
 29:   if (ctx->intb >= PETSC_MAX_REAL && ctx->inta <= PETSC_MIN_REAL) SETERRQ(PetscObjectComm((PetscObject)st),1,"Must pass an interval with STFilterSetInterval()");
 30:   if (ctx->right == 0.0 && ctx->left == 0.0) SETERRQ(PetscObjectComm((PetscObject)st),1,"Must pass an approximate numerical range with STFilterSetRange()");
 31:   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);
 32:   if (!ctx->polyDegree) ctx->polyDegree = 100;
 33:   ctx->frame[0] = ctx->left;
 34:   ctx->frame[1] = ctx->inta;
 35:   ctx->frame[2] = ctx->intb;
 36:   ctx->frame[3] = ctx->right;
 37:   STFilter_FILTLAN_setFilter(st,&st->T[0]);
 38:   st->M = st->T[0];
 39:   MatDestroy(&st->P);
 40:   return(0);
 41: }

 43: PetscErrorCode STSetUp_Filter(ST st)
 44: {

 48:   STSetWorkVecs(st,4);
 49:   return(0);
 50: }

 52: PetscErrorCode STSetFromOptions_Filter(PetscOptionItems *PetscOptionsObject,ST st)
 53: {
 55:   PetscReal      array[2]={0,0};
 56:   PetscInt       k;
 57:   PetscBool      flg;

 60:   PetscOptionsHead(PetscOptionsObject,"ST Filter Options");

 62:     k = 2;
 63:     PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",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_interval (comma-separated without spaces)");
 66:       STFilterSetInterval(st,array[0],array[1]);
 67:     }
 68:     k = 2;
 69:     PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg);
 70:     if (flg) {
 71:       if (k<2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_range (comma-separated without spaces)");
 72:       STFilterSetRange(st,array[0],array[1]);
 73:     }
 74:     PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg);
 75:     if (flg) { STFilterSetDegree(st,k); }

 77:   PetscOptionsTail();
 78:   return(0);
 79: }

 81: static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
 82: {
 83:   ST_FILTER *ctx = (ST_FILTER*)st->data;

 86:   if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
 87:   if (ctx->inta != inta || ctx->intb != intb) {
 88:     ctx->inta   = inta;
 89:     ctx->intb   = intb;
 90:     st->state   = ST_STATE_INITIAL;
 91:     st->opready = PETSC_FALSE;
 92:     ctx->filtch = PETSC_TRUE;
 93:   }
 94:   return(0);
 95: }

 97: /*@
 98:    STFilterSetInterval - Defines the interval containing the desired eigenvalues.

100:    Logically Collective on st

102:    Input Parameters:
103: +  st   - the spectral transformation context
104: .  inta - left end of the interval
105: -  intb - right end of the interval

107:    Options Database Key:
108: .  -st_filter_interval <a,b> - set [a,b] as the interval of interest

110:    Level: intermediate

112:    Notes:
113:    The filter will be configured to emphasize eigenvalues contained in the given
114:    interval, and damp out eigenvalues outside it. If the interval is open, then
115:    the filter is low- or high-pass, otherwise it is mid-pass.

117:    Common usage is to set the interval in EPS with EPSSetInterval().

119:    The interval must be contained within the numerical range of the matrix, see
120:    STFilterSetRange().

122: .seealso: STFilterGetInterval(), STFilterSetRange(), EPSSetInterval()
123: @*/
124: PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
125: {

132:   PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
133:   return(0);
134: }

136: static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
137: {
138:   ST_FILTER *ctx = (ST_FILTER*)st->data;

141:   if (inta) *inta = ctx->inta;
142:   if (intb) *intb = ctx->intb;
143:   return(0);
144: }

146: /*@
147:    STFilterGetInterval - Gets the interval containing the desired eigenvalues.

149:    Not Collective

151:    Input Parameter:
152: .  st  - the spectral transformation context

154:    Output Parameter:
155: +  inta - left end of the interval
156: -  intb - right end of the interval

158:    Level: intermediate

160: .seealso: STFilterSetInterval()
161: @*/
162: PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
163: {

168:   PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
169:   return(0);
170: }

172: static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
173: {
174:   ST_FILTER *ctx = (ST_FILTER*)st->data;

177:   if (left>right) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be left<right");
178:   if (ctx->left != left || ctx->right != right) {
179:     ctx->left   = left;
180:     ctx->right  = right;
181:     st->state   = ST_STATE_INITIAL;
182:     st->opready = PETSC_FALSE;
183:     ctx->filtch = PETSC_TRUE;
184:   }
185:   return(0);
186: }

188: /*@
189:    STFilterSetRange - Defines the numerical range (or field of values) of the matrix, that is,
190:    the interval containing all eigenvalues.

192:    Logically Collective on st

194:    Input Parameters:
195: +  st    - the spectral transformation context
196: .  left  - left end of the interval
197: -  right - right end of the interval

199:    Options Database Key:
200: .  -st_filter_range <a,b> - set [a,b] as the numerical range

202:    Level: intermediate

204:    Notes:
205:    The filter will be most effective if the numerical range is tight, that is, left and right
206:    are good approximations to the leftmost and rightmost eigenvalues, respectively.

208: .seealso: STFilterGetRange(), STFilterSetInterval()
209: @*/
210: PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
211: {

218:   PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
219:   return(0);
220: }

222: static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
223: {
224:   ST_FILTER *ctx = (ST_FILTER*)st->data;

227:   if (left)  *left  = ctx->left;
228:   if (right) *right = ctx->right;
229:   return(0);
230: }

232: /*@
233:    STFilterGetRange - Gets the interval containing all eigenvalues.

235:    Not Collective

237:    Input Parameter:
238: .  st  - the spectral transformation context

240:    Output Parameter:
241: +  left  - left end of the interval
242: -  right - right end of the interval

244:    Level: intermediate

246: .seealso: STFilterSetRange()
247: @*/
248: PetscErrorCode STFilterGetRange(ST st,PetscReal *left,PetscReal *right)
249: {

254:   PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
255:   return(0);
256: }

258: static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
259: {
260:   ST_FILTER *ctx = (ST_FILTER*)st->data;

263:   if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
264:     ctx->polyDegree = 0;
265:     st->state       = ST_STATE_INITIAL;
266:     st->opready     = PETSC_FALSE;
267:     ctx->filtch     = PETSC_TRUE;
268:   } else {
269:     if (deg<=0) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
270:     if (ctx->polyDegree != deg) {
271:       ctx->polyDegree = deg;
272:       st->state       = ST_STATE_INITIAL;
273:       st->opready     = PETSC_FALSE;
274:       ctx->filtch     = PETSC_TRUE;
275:     }
276:   }
277:   return(0);
278: }

280: /*@
281:    STFilterSetDegree - Sets the degree of the filter polynomial.

283:    Logically Collective on st

285:    Input Parameters:
286: +  st  - the spectral transformation context
287: -  deg - polynomial degree

289:    Options Database Key:
290: .  -st_filter_degree <deg> - sets the degree of the filter polynomial

292:    Level: intermediate

294: .seealso: STFilterGetDegree()
295: @*/
296: PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
297: {

303:   PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
304:   return(0);
305: }

307: static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
308: {
309:   ST_FILTER *ctx = (ST_FILTER*)st->data;

312:   *deg = ctx->polyDegree;
313:   return(0);
314: }

316: /*@
317:    STFilterGetDegree - Gets the degree of the filter polynomial.

319:    Not Collective

321:    Input Parameter:
322: .  st  - the spectral transformation context

324:    Output Parameter:
325: .  deg - polynomial degree

327:    Level: intermediate

329: .seealso: STFilterSetDegree()
330: @*/
331: PetscErrorCode STFilterGetDegree(ST st,PetscInt *deg)
332: {

338:   PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
339:   return(0);
340: }

342: static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
343: {
344:   ST_FILTER *ctx = (ST_FILTER*)st->data;

347:   *gamma = ctx->filterInfo->yLimit;
348:   return(0);
349: }

351: /*@
352:    STFilterGetThreshold - Gets the threshold value gamma such that rho(lambda)>=gamma for
353:    eigenvalues lambda inside the wanted interval and rho(lambda)<gamma for those outside.

355:    Not Collective

357:    Input Parameter:
358: .  st  - the spectral transformation context

360:    Output Parameter:
361: .  gamma - the threshold value

363:    Level: developer
364: @*/
365: PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
366: {

372:   PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
373:   return(0);
374: }

376: PetscErrorCode STReset_Filter(ST st)
377: {
379:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

382:   ctx->left  = 0.0;
383:   ctx->right = 0.0;
384:   MatDestroy(&ctx->T);
385:   return(0);
386: }

388: PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
389: {
391:   ST_FILTER      *ctx = (ST_FILTER*)st->data;
392:   PetscBool      isascii;

395:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
396:   if (isascii) {
397:     PetscViewerASCIIPrintf(viewer,"  Filter: interval of desired eigenvalues is [%g,%g]\n",(double)ctx->inta,(double)ctx->intb);
398:     PetscViewerASCIIPrintf(viewer,"  Filter: numerical range is [%g,%g]\n",(double)ctx->left,(double)ctx->right);
399:     PetscViewerASCIIPrintf(viewer,"  Filter: degree of filter polynomial is %D\n",ctx->polyDegree);
400:     if (st->state>=ST_STATE_SETUP) {
401:       PetscViewerASCIIPrintf(viewer,"  Filter: limit to accept eigenvalues: theta=%g\n",ctx->filterInfo->yLimit);
402:     }
403:   }
404:   return(0);
405: }

407: PetscErrorCode STDestroy_Filter(ST st)
408: {
410:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

413:   PetscFree(ctx->opts);
414:   PetscFree(ctx->filterInfo);
415:   PetscFree(ctx->baseFilter);
416:   PetscFree(st->data);
417:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL);
418:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL);
419:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL);
420:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL);
421:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL);
422:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL);
423:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL);
424:   return(0);
425: }

427: SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
428: {
430:   ST_FILTER      *ctx;
431:   FILTLAN_IOP    iop;
432:   FILTLAN_PFI    pfi;
434:   PetscNewLog(st,&ctx);
435:   st->data = (void*)ctx;

437:   st->usesksp = PETSC_FALSE;

439:   ctx->inta               = PETSC_MIN_REAL;
440:   ctx->intb               = PETSC_MAX_REAL;
441:   ctx->left               = 0.0;
442:   ctx->right              = 0.0;
443:   ctx->polyDegree         = 0;
444:   ctx->baseDegree         = 10;

446:   PetscNewLog(st,&iop);
447:   ctx->opts               = iop;
448:   iop->intervalWeights[0] = 100.0;
449:   iop->intervalWeights[1] = 1.0;
450:   iop->intervalWeights[2] = 1.0;
451:   iop->intervalWeights[3] = 1.0;
452:   iop->intervalWeights[4] = 100.0;
453:   iop->transIntervalRatio = 0.6;
454:   iop->reverseInterval    = PETSC_FALSE;
455:   iop->initialPlateau     = 0.1;
456:   iop->plateauShrinkRate  = 1.5;
457:   iop->initialShiftStep   = 0.01;
458:   iop->shiftStepExpanRate = 1.5;
459:   iop->maxInnerIter       = 30;
460:   iop->yLimitTol          = 0.0001;
461:   iop->maxOuterIter       = 50;
462:   iop->numGridPoints      = 200;
463:   iop->yBottomLine        = 0.001;
464:   iop->yRippleLimit       = 0.8;

466:   PetscNewLog(st,&pfi);
467:   ctx->filterInfo         = pfi;

469:   st->ops->apply           = STApply_Generic;
470:   st->ops->setup           = STSetUp_Filter;
471:   st->ops->computeoperator = STComputeOperator_Filter;
472:   st->ops->setfromoptions  = STSetFromOptions_Filter;
473:   st->ops->destroy         = STDestroy_Filter;
474:   st->ops->reset           = STReset_Filter;
475:   st->ops->view            = STView_Filter;

477:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter);
478:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter);
479:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter);
480:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter);
481:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter);
482:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter);
483:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter);
484:   return(0);
485: }