Actual source code: filter.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  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: PetscErrorCode STComputeOperator_Filter(ST st)
 23: {
 24:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

 30:   if (!ctx->polyDegree) ctx->polyDegree = 100;
 31:   ctx->frame[0] = ctx->left;
 32:   ctx->frame[1] = ctx->inta;
 33:   ctx->frame[2] = ctx->intb;
 34:   ctx->frame[3] = ctx->right;
 35:   STFilter_FILTLAN_setFilter(st,&st->T[0]);
 36:   st->M = st->T[0];
 37:   MatDestroy(&st->P);
 38:   return 0;
 39: }

 41: PetscErrorCode STSetUp_Filter(ST st)
 42: {
 43:   STSetWorkVecs(st,4);
 44:   return 0;
 45: }

 47: PetscErrorCode STSetFromOptions_Filter(ST st,PetscOptionItems *PetscOptionsObject)
 48: {
 49:   PetscReal      array[2]={0,0};
 50:   PetscInt       k;
 51:   PetscBool      flg;

 53:   PetscOptionsHeadBegin(PetscOptionsObject,"ST Filter Options");

 55:     k = 2;
 56:     PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",array,&k,&flg);
 57:     if (flg) {
 59:       STFilterSetInterval(st,array[0],array[1]);
 60:     }
 61:     k = 2;
 62:     PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg);
 63:     if (flg) {
 65:       STFilterSetRange(st,array[0],array[1]);
 66:     }
 67:     PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg);
 68:     if (flg) STFilterSetDegree(st,k);

 70:   PetscOptionsHeadEnd();
 71:   return 0;
 72: }

 74: static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
 75: {
 76:   ST_FILTER *ctx = (ST_FILTER*)st->data;

 79:   if (ctx->inta != inta || ctx->intb != intb) {
 80:     ctx->inta   = inta;
 81:     ctx->intb   = intb;
 82:     st->state   = ST_STATE_INITIAL;
 83:     st->opready = PETSC_FALSE;
 84:     ctx->filtch = PETSC_TRUE;
 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:    Notes:
103:    The filter will be configured to emphasize eigenvalues contained in the given
104:    interval, and damp out eigenvalues outside it. If the interval is open, then
105:    the filter is low- or high-pass, otherwise it is mid-pass.

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

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

112:    Level: intermediate

114: .seealso: STFilterGetInterval(), STFilterSetRange(), EPSSetInterval()
115: @*/
116: PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
117: {
121:   PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
122:   return 0;
123: }

125: static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
126: {
127:   ST_FILTER *ctx = (ST_FILTER*)st->data;

129:   if (inta) *inta = ctx->inta;
130:   if (intb) *intb = ctx->intb;
131:   return 0;
132: }

134: /*@
135:    STFilterGetInterval - Gets the interval containing the desired eigenvalues.

137:    Not Collective

139:    Input Parameter:
140: .  st  - the spectral transformation context

142:    Output Parameters:
143: +  inta - left end of the interval
144: -  intb - right end of the interval

146:    Level: intermediate

148: .seealso: STFilterSetInterval()
149: @*/
150: PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
151: {
153:   PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
154:   return 0;
155: }

157: static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
158: {
159:   ST_FILTER *ctx = (ST_FILTER*)st->data;

162:   if (ctx->left != left || ctx->right != right) {
163:     ctx->left   = left;
164:     ctx->right  = right;
165:     st->state   = ST_STATE_INITIAL;
166:     st->opready = PETSC_FALSE;
167:     ctx->filtch = PETSC_TRUE;
168:   }
169:   return 0;
170: }

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

176:    Logically Collective on st

178:    Input Parameters:
179: +  st    - the spectral transformation context
180: .  left  - left end of the interval
181: -  right - right end of the interval

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

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

190:    Level: intermediate

192: .seealso: STFilterGetRange(), STFilterSetInterval()
193: @*/
194: PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
195: {
199:   PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
200:   return 0;
201: }

203: static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
204: {
205:   ST_FILTER *ctx = (ST_FILTER*)st->data;

207:   if (left)  *left  = ctx->left;
208:   if (right) *right = ctx->right;
209:   return 0;
210: }

212: /*@
213:    STFilterGetRange - Gets the interval containing all eigenvalues.

215:    Not Collective

217:    Input Parameter:
218: .  st  - the spectral transformation context

220:    Output Parameters:
221: +  left  - left end of the interval
222: -  right - right end of the interval

224:    Level: intermediate

226: .seealso: STFilterSetRange()
227: @*/
228: PetscErrorCode STFilterGetRange(ST st,PetscReal *left,PetscReal *right)
229: {
231:   PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
232:   return 0;
233: }

235: static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
236: {
237:   ST_FILTER *ctx = (ST_FILTER*)st->data;

239:   if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
240:     ctx->polyDegree = 0;
241:     st->state       = ST_STATE_INITIAL;
242:     st->opready     = PETSC_FALSE;
243:     ctx->filtch     = PETSC_TRUE;
244:   } else {
246:     if (ctx->polyDegree != deg) {
247:       ctx->polyDegree = deg;
248:       st->state       = ST_STATE_INITIAL;
249:       st->opready     = PETSC_FALSE;
250:       ctx->filtch     = PETSC_TRUE;
251:     }
252:   }
253:   return 0;
254: }

256: /*@
257:    STFilterSetDegree - Sets the degree of the filter polynomial.

259:    Logically Collective on st

261:    Input Parameters:
262: +  st  - the spectral transformation context
263: -  deg - polynomial degree

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

268:    Level: intermediate

270: .seealso: STFilterGetDegree()
271: @*/
272: PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
273: {
276:   PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
277:   return 0;
278: }

280: static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
281: {
282:   ST_FILTER *ctx = (ST_FILTER*)st->data;

284:   *deg = ctx->polyDegree;
285:   return 0;
286: }

288: /*@
289:    STFilterGetDegree - Gets the degree of the filter polynomial.

291:    Not Collective

293:    Input Parameter:
294: .  st  - the spectral transformation context

296:    Output Parameter:
297: .  deg - polynomial degree

299:    Level: intermediate

301: .seealso: STFilterSetDegree()
302: @*/
303: PetscErrorCode STFilterGetDegree(ST st,PetscInt *deg)
304: {
307:   PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
308:   return 0;
309: }

311: static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
312: {
313:   ST_FILTER *ctx = (ST_FILTER*)st->data;

315:   *gamma = ctx->filterInfo->yLimit;
316:   return 0;
317: }

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

323:    Not Collective

325:    Input Parameter:
326: .  st  - the spectral transformation context

328:    Output Parameter:
329: .  gamma - the threshold value

331:    Level: developer

333: .seealso: STFilterGetRange()
334: @*/
335: PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
336: {
339:   PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
340:   return 0;
341: }

343: PetscErrorCode STReset_Filter(ST st)
344: {
345:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

347:   ctx->left  = 0.0;
348:   ctx->right = 0.0;
349:   MatDestroy(&ctx->T);
350:   MatDestroyMatrices(ctx->nW,&ctx->W);
351:   ctx->nW = 0;
352:   return 0;
353: }

355: PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
356: {
357:   ST_FILTER      *ctx = (ST_FILTER*)st->data;
358:   PetscBool      isascii;

360:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
361:   if (isascii) {
362:     PetscViewerASCIIPrintf(viewer,"  interval of desired eigenvalues: [%g,%g]\n",(double)ctx->inta,(double)ctx->intb);
363:     PetscViewerASCIIPrintf(viewer,"  numerical range: [%g,%g]\n",(double)ctx->left,(double)ctx->right);
364:     PetscViewerASCIIPrintf(viewer,"  degree of filter polynomial: %" PetscInt_FMT "\n",ctx->polyDegree);
365:     if (st->state>=ST_STATE_SETUP) PetscViewerASCIIPrintf(viewer,"  limit to accept eigenvalues: theta=%g\n",(double)ctx->filterInfo->yLimit);
366:   }
367:   return 0;
368: }

370: PetscErrorCode STDestroy_Filter(ST st)
371: {
372:   ST_FILTER      *ctx = (ST_FILTER*)st->data;

374:   PetscFree(ctx->opts);
375:   PetscFree(ctx->filterInfo);
376:   PetscFree(ctx->baseFilter);
377:   PetscFree(st->data);
378:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL);
379:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL);
380:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL);
381:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL);
382:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL);
383:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL);
384:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL);
385:   return 0;
386: }

388: SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
389: {
390:   ST_FILTER      *ctx;
391:   FILTLAN_IOP    iop;
392:   FILTLAN_PFI    pfi;
393:   PetscNew(&ctx);
394:   st->data = (void*)ctx;

396:   st->usesksp = PETSC_FALSE;

398:   ctx->inta               = PETSC_MIN_REAL;
399:   ctx->intb               = PETSC_MAX_REAL;
400:   ctx->left               = 0.0;
401:   ctx->right              = 0.0;
402:   ctx->polyDegree         = 0;
403:   ctx->baseDegree         = 10;

405:   PetscNew(&iop);
406:   ctx->opts               = iop;
407:   iop->intervalWeights[0] = 100.0;
408:   iop->intervalWeights[1] = 1.0;
409:   iop->intervalWeights[2] = 1.0;
410:   iop->intervalWeights[3] = 1.0;
411:   iop->intervalWeights[4] = 100.0;
412:   iop->transIntervalRatio = 0.6;
413:   iop->reverseInterval    = PETSC_FALSE;
414:   iop->initialPlateau     = 0.1;
415:   iop->plateauShrinkRate  = 1.5;
416:   iop->initialShiftStep   = 0.01;
417:   iop->shiftStepExpanRate = 1.5;
418:   iop->maxInnerIter       = 30;
419:   iop->yLimitTol          = 0.0001;
420:   iop->maxOuterIter       = 50;
421:   iop->numGridPoints      = 200;
422:   iop->yBottomLine        = 0.001;
423:   iop->yRippleLimit       = 0.8;

425:   PetscNew(&pfi);
426:   ctx->filterInfo         = pfi;

428:   st->ops->apply           = STApply_Generic;
429:   st->ops->setup           = STSetUp_Filter;
430:   st->ops->computeoperator = STComputeOperator_Filter;
431:   st->ops->setfromoptions  = STSetFromOptions_Filter;
432:   st->ops->destroy         = STDestroy_Filter;
433:   st->ops->reset           = STReset_Filter;
434:   st->ops->view            = STView_Filter;

436:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter);
437:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter);
438:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter);
439:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter);
440:   PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter);
441:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter);
442:   PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter);
443:   return 0;
444: }