Actual source code: filter.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }