Actual source code: filter.c

slepc-main 2024-11-09
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: 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: }