Actual source code: filter.c

  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: const char *STFilterTypes[] = {"","FILTLAN","CHEBYSHEV","STFilterType","ST_FILTER_",NULL};
 18: const char *STFilterDampings[] = {"NONE","JACKSON","LANCZOS","FEJER","STFilterDamping","ST_FILTER_DAMPING_",NULL};

 20: static PetscErrorCode STFilterSetType_Private(ST st,STFilterType type)
 21: {
 22:   ST_FILTER    *ctx = (ST_FILTER*)st->data;

 24:   PetscFunctionBegin;
 25:   ctx->type = type;
 26:   switch(type) {
 27:     case ST_FILTER_FILTLAN:
 28:       PetscCall(STCreate_Filter_FILTLAN(st));
 29:       break;
 30:     case ST_FILTER_CHEBYSHEV:
 31:       PetscCall(STCreate_Filter_Chebyshev(st));
 32:       break;
 33:     default:
 34:       SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Invalid type");
 35:   }
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: /*
 40:    Operator (filter):
 41:                Op               P         M
 42:    if nmat=1:  p(A)             NULL      p(A)
 43: */
 44: static PetscErrorCode STComputeOperator_Filter(ST st)
 45: {
 46:   ST_FILTER *ctx = (ST_FILTER*)st->data;

 48:   PetscFunctionBegin;
 49:   PetscCheck(st->nmat==1,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Only implemented for standard eigenvalue problem");
 50:   PetscCheck(ctx->intb<PETSC_MAX_REAL || ctx->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an interval with STFilterSetInterval()");
 51:   PetscCheck(ctx->right!=0.0 || ctx->left!=0.0,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"Must pass an approximate numerical range with STFilterSetRange()");
 52:   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);

 54:   if (!ctx->type) PetscCall(STFilterSetType_Private(st,ST_FILTER_FILTLAN)); /* default type */
 55:   if (!ctx->polyDegree) ctx->polyDegree = 100;
 56:   PetscCall(ctx->computeoperator(st,&st->T[0]));
 57:   st->M = st->T[0];
 58:   PetscCall(MatDestroy(&st->P));
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: static PetscErrorCode STSetFromOptions_Filter(ST st,PetscOptionItems PetscOptionsObject)
 63: {
 64:   ST_FILTER        *ctx = (ST_FILTER*)st->data;
 65:   PetscReal        array[2]={0,0};
 66:   PetscInt         k;
 67:   PetscBool        flg;
 68:   STFilterType     type;
 69:   STFilterDamping  damping;

 71:   PetscFunctionBegin;
 72:   PetscOptionsHeadBegin(PetscOptionsObject,"ST Filter Options");

 74:     PetscCall(PetscOptionsEnum("-st_filter_type","How to construct the filter","STFilterSetType",STFilterTypes,(PetscEnum)ctx->type,(PetscEnum*)&type,&flg));
 75:     if (flg) PetscCall(STFilterSetType(st,type));

 77:     k = 2;
 78:     PetscCall(PetscOptionsRealArray("-st_filter_interval","Interval containing the desired eigenvalues (two real values separated with a comma without spaces)","STFilterSetInterval",array,&k,&flg));
 79:     if (flg) {
 80:       PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_interval (comma-separated without spaces)");
 81:       PetscCall(STFilterSetInterval(st,array[0],array[1]));
 82:     }
 83:     k = 2;
 84:     PetscCall(PetscOptionsRealArray("-st_filter_range","Interval containing all eigenvalues (two real values separated with a comma without spaces)","STFilterSetRange",array,&k,&flg));
 85:     if (flg) {
 86:       PetscCheck(k>1,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Must pass two values in -st_filter_range (comma-separated without spaces)");
 87:       PetscCall(STFilterSetRange(st,array[0],array[1]));
 88:     }
 89:     PetscCall(PetscOptionsInt("-st_filter_degree","Degree of filter polynomial","STFilterSetDegree",100,&k,&flg));
 90:     if (flg) PetscCall(STFilterSetDegree(st,k));

 92:     PetscCall(PetscOptionsEnum("-st_filter_damping","Type of damping","STFilterSetDamping",STFilterDampings,(PetscEnum)ctx->damping,(PetscEnum*)&damping,&flg));
 93:     if (flg) PetscCall(STFilterSetDamping(st,damping));

 95:   PetscOptionsHeadEnd();
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: static PetscErrorCode STReset_Filter(ST st)
100: {
101:   ST_FILTER *ctx = (ST_FILTER*)st->data;

103:   PetscFunctionBegin;
104:   ctx->left  = 0.0;
105:   ctx->right = 0.0;
106:   PetscCall(MatDestroy(&ctx->T));
107:   PetscCall(MatDestroyMatrices(ctx->nW,&ctx->W));
108:   ctx->nW = 0;
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }

112: static PetscErrorCode STFilterSetType_Filter(ST st,STFilterType type)
113: {
114:   ST_FILTER *ctx = (ST_FILTER*)st->data;

116:   PetscFunctionBegin;
117:   if (ctx->type != type) {
118:     PetscCall(STReset_Filter(st));
119:     PetscCall(STFilterSetType_Private(st,type));
120:     st->state   = ST_STATE_INITIAL;
121:     st->opready = PETSC_FALSE;
122:     ctx->filtch = PETSC_TRUE;
123:   }
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: /*@
128:    STFilterSetType - Sets the method to be used to build the polynomial filter.

130:    Logically Collective

132:    Input Parameters:
133: +  st   - the spectral transformation context
134: -  type - the type of filter

136:    Options Database Key:
137: .  -st_filter_type <type> - set the type of filter

139:    Level: intermediate

141: .seealso: [](ch:st), `STFILTER`, `STFilterGetType()`
142: @*/
143: PetscErrorCode STFilterSetType(ST st,STFilterType type)
144: {
145:   PetscFunctionBegin;
148:   PetscTryMethod(st,"STFilterSetType_C",(ST,STFilterType),(st,type));
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: static PetscErrorCode STFilterGetType_Filter(ST st,STFilterType *type)
153: {
154:   ST_FILTER *ctx = (ST_FILTER*)st->data;

156:   PetscFunctionBegin;
157:   *type = ctx->type ;
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: /*@
162:    STFilterGetType - Gets the method to be used to build the polynomial filter.

164:    Not Collective

166:    Input Parameter:
167: .  st  - the spectral transformation context

169:    Output Parameter:
170: .  type - the type of filter

172:    Level: intermediate

174: .seealso: [](ch:st), `STFilterSetType()`
175: @*/
176: PetscErrorCode STFilterGetType(ST st,STFilterType *type)
177: {
178:   PetscFunctionBegin;
180:   PetscAssertPointer(type,2);
181:   PetscUseMethod(st,"STFilterGetType_C",(ST,STFilterType*),(st,type));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: static PetscErrorCode STFilterSetInterval_Filter(ST st,PetscReal inta,PetscReal intb)
186: {
187:   ST_FILTER *ctx = (ST_FILTER*)st->data;

189:   PetscFunctionBegin;
190:   PetscCheck(inta<intb,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
191:   if (ctx->inta != inta || ctx->intb != intb) {
192:     ctx->inta   = inta;
193:     ctx->intb   = intb;
194:     st->state   = ST_STATE_INITIAL;
195:     st->opready = PETSC_FALSE;
196:     ctx->filtch = PETSC_TRUE;
197:   }
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: /*@
202:    STFilterSetInterval - Defines the interval containing the desired eigenvalues.

204:    Logically Collective

206:    Input Parameters:
207: +  st   - the spectral transformation context
208: .  inta - left end of the interval
209: -  intb - right end of the interval

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

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

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

221:    The interval must be contained within the numerical range of the matrix, see
222:    `STFilterSetRange()`.

224:    Level: intermediate

226: .seealso: [](ch:st), `STFILTER`, `STFilterGetInterval()`, `STFilterSetRange()`, `EPSSetInterval()`
227: @*/
228: PetscErrorCode STFilterSetInterval(ST st,PetscReal inta,PetscReal intb)
229: {
230:   PetscFunctionBegin;
234:   PetscTryMethod(st,"STFilterSetInterval_C",(ST,PetscReal,PetscReal),(st,inta,intb));
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: static PetscErrorCode STFilterGetInterval_Filter(ST st,PetscReal *inta,PetscReal *intb)
239: {
240:   ST_FILTER *ctx = (ST_FILTER*)st->data;

242:   PetscFunctionBegin;
243:   if (inta) *inta = ctx->inta;
244:   if (intb) *intb = ctx->intb;
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: /*@
249:    STFilterGetInterval - Gets the interval containing the desired eigenvalues.

251:    Not Collective

253:    Input Parameter:
254: .  st  - the spectral transformation context

256:    Output Parameters:
257: +  inta - left end of the interval
258: -  intb - right end of the interval

260:    Level: intermediate

262: .seealso: [](ch:st), `STFilterSetInterval()`
263: @*/
264: PetscErrorCode STFilterGetInterval(ST st,PetscReal *inta,PetscReal *intb)
265: {
266:   PetscFunctionBegin;
268:   PetscUseMethod(st,"STFilterGetInterval_C",(ST,PetscReal*,PetscReal*),(st,inta,intb));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: static PetscErrorCode STFilterSetRange_Filter(ST st,PetscReal left,PetscReal right)
273: {
274:   ST_FILTER *ctx = (ST_FILTER*)st->data;

276:   PetscFunctionBegin;
277:   PetscCheck(left<right,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be left<right");
278:   if (ctx->left != left || ctx->right != right) {
279:     ctx->left   = left;
280:     ctx->right  = right;
281:     st->state   = ST_STATE_INITIAL;
282:     st->opready = PETSC_FALSE;
283:     ctx->filtch = PETSC_TRUE;
284:   }
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

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

292:    Logically Collective

294:    Input Parameters:
295: +  st    - the spectral transformation context
296: .  left  - left end of the spectral range
297: -  right - right end of the spectral range

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

302:    Notes:
303:    The filter will be most effective if the numerical range is tight, that is,
304:    `left` and `right` are good approximations to the leftmost and rightmost
305:    eigenvalues, respectively.

307:    Level: intermediate

309: .seealso: [](ch:st), `STFILTER`, `STFilterGetRange()`, `STFilterSetInterval()`
310: @*/
311: PetscErrorCode STFilterSetRange(ST st,PetscReal left,PetscReal right)
312: {
313:   PetscFunctionBegin;
317:   PetscTryMethod(st,"STFilterSetRange_C",(ST,PetscReal,PetscReal),(st,left,right));
318:   PetscFunctionReturn(PETSC_SUCCESS);
319: }

321: static PetscErrorCode STFilterGetRange_Filter(ST st,PetscReal *left,PetscReal *right)
322: {
323:   ST_FILTER *ctx = (ST_FILTER*)st->data;

325:   PetscFunctionBegin;
326:   if (left)  *left  = ctx->left;
327:   if (right) *right = ctx->right;
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: /*@
332:    STFilterGetRange - Gets the interval containing all eigenvalues.

334:    Not Collective

336:    Input Parameter:
337: .  st  - the spectral transformation context

339:    Output Parameters:
340: +  left  - left end of the spectral range
341: -  right - right end of the spectral range

343:    Level: intermediate

345: .seealso: [](ch:st), `STFilterSetRange()`
346: @*/
347: PetscErrorCode STFilterGetRange(ST st,PetscReal *left,PetscReal *right)
348: {
349:   PetscFunctionBegin;
351:   PetscUseMethod(st,"STFilterGetRange_C",(ST,PetscReal*,PetscReal*),(st,left,right));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: static PetscErrorCode STFilterSetDegree_Filter(ST st,PetscInt deg)
356: {
357:   ST_FILTER *ctx = (ST_FILTER*)st->data;

359:   PetscFunctionBegin;
360:   if (deg == PETSC_DEFAULT || deg == PETSC_DECIDE) {
361:     ctx->polyDegree = 0;
362:     st->state       = ST_STATE_INITIAL;
363:     st->opready     = PETSC_FALSE;
364:     ctx->filtch     = PETSC_TRUE;
365:   } else {
366:     PetscCheck(deg>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
367:     if (ctx->polyDegree != deg) {
368:       ctx->polyDegree = deg;
369:       st->state       = ST_STATE_INITIAL;
370:       st->opready     = PETSC_FALSE;
371:       ctx->filtch     = PETSC_TRUE;
372:     }
373:   }
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }

377: /*@
378:    STFilterSetDegree - Sets the degree of the filter polynomial.

380:    Logically Collective

382:    Input Parameters:
383: +  st  - the spectral transformation context
384: -  deg - polynomial degree

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

389:    Level: intermediate

391: .seealso: [](ch:st), `STFILTER`, `STFilterGetDegree()`
392: @*/
393: PetscErrorCode STFilterSetDegree(ST st,PetscInt deg)
394: {
395:   PetscFunctionBegin;
398:   PetscTryMethod(st,"STFilterSetDegree_C",(ST,PetscInt),(st,deg));
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: static PetscErrorCode STFilterGetDegree_Filter(ST st,PetscInt *deg)
403: {
404:   ST_FILTER *ctx = (ST_FILTER*)st->data;

406:   PetscFunctionBegin;
407:   *deg = ctx->polyDegree;
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }

411: /*@
412:    STFilterGetDegree - Gets the degree of the filter polynomial.

414:    Not Collective

416:    Input Parameter:
417: .  st  - the spectral transformation context

419:    Output Parameter:
420: .  deg - polynomial degree

422:    Level: intermediate

424: .seealso: [](ch:st), `STFilterSetDegree()`
425: @*/
426: PetscErrorCode STFilterGetDegree(ST st,PetscInt *deg)
427: {
428:   PetscFunctionBegin;
430:   PetscAssertPointer(deg,2);
431:   PetscUseMethod(st,"STFilterGetDegree_C",(ST,PetscInt*),(st,deg));
432:   PetscFunctionReturn(PETSC_SUCCESS);
433: }

435: static PetscErrorCode STFilterGetThreshold_Filter(ST st,PetscReal *gamma)
436: {
437:   ST_FILTER *ctx = (ST_FILTER*)st->data;

439:   PetscFunctionBegin;
440:   if (ctx->getthreshold) PetscCall(ctx->getthreshold(st,gamma));
441:   else *gamma = 0.5;
442:   PetscFunctionReturn(PETSC_SUCCESS);
443: }

445: /*@
446:    STFilterGetThreshold - Gets the threshold value $\gamma$ used to decide
447:    whether eigenvalue approximations are inside or outside the wanted interval.

449:    Not Collective

451:    Input Parameter:
452: .  st  - the spectral transformation context

454:    Output Parameter:
455: .  gamma - the threshold value

457:    Note:
458:    An eigenvalue $\lambda$ is considered to be inside the wanted interval
459:    if $|p(\lambda)|\geq\gamma$, and outside otherwise.

461:    Level: developer

463: .seealso: [](ch:st), `STFILTER`, `STFilterGetRange()`
464: @*/
465: PetscErrorCode STFilterGetThreshold(ST st,PetscReal *gamma)
466: {
467:   PetscFunctionBegin;
469:   PetscAssertPointer(gamma,2);
470:   PetscUseMethod(st,"STFilterGetThreshold_C",(ST,PetscReal*),(st,gamma));
471:   PetscFunctionReturn(PETSC_SUCCESS);
472: }

474: static PetscErrorCode STFilterSetDamping_Filter(ST st,STFilterDamping damping)
475: {
476:   ST_FILTER *ctx = (ST_FILTER*)st->data;

478:   PetscFunctionBegin;
479:   if (ctx->damping != damping) {
480:     ctx->damping = damping;
481:     st->state    = ST_STATE_INITIAL;
482:     st->opready  = PETSC_FALSE;
483:     ctx->filtch  = PETSC_TRUE;
484:   }
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: /*@
489:    STFilterSetDamping - Sets the type of damping to be used in the polynomial filter.

491:    Logically Collective

493:    Input Parameters:
494: +  st      - the spectral transformation context
495: -  damping - the type of damping

497:    Options Database Key:
498: .  -st_filter_damping <damping> - sets the type of damping

500:    Note:
501:    Only used in `ST_FILTER_CHEBYSHEV` filters.

503:    Level: advanced

505: .seealso: [](ch:st), `STFILTER`, `STFilterGetDamping()`
506: @*/
507: PetscErrorCode STFilterSetDamping(ST st,STFilterDamping damping)
508: {
509:   PetscFunctionBegin;
512:   PetscTryMethod(st,"STFilterSetDamping_C",(ST,STFilterDamping),(st,damping));
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: static PetscErrorCode STFilterGetDamping_Filter(ST st,STFilterDamping *damping)
517: {
518:   ST_FILTER *ctx = (ST_FILTER*)st->data;

520:   PetscFunctionBegin;
521:   *damping = ctx->damping;
522:   PetscFunctionReturn(PETSC_SUCCESS);
523: }

525: /*@
526:    STFilterGetDamping - Gets the type of damping used in the polynomial filter.

528:    Not Collective

530:    Input Parameter:
531: .  st  - the spectral transformation context

533:    Output Parameter:
534: .  damping - the type of damping

536:    Level: advanced

538: .seealso: [](ch:st), `STFilterSetDamping()`
539: @*/
540: PetscErrorCode STFilterGetDamping(ST st,STFilterDamping *damping)
541: {
542:   PetscFunctionBegin;
544:   PetscAssertPointer(damping,2);
545:   PetscUseMethod(st,"STFilterGetDamping_C",(ST,STFilterDamping*),(st,damping));
546:   PetscFunctionReturn(PETSC_SUCCESS);
547: }

549: static PetscErrorCode STView_Filter(ST st,PetscViewer viewer)
550: {
551:   ST_FILTER *ctx = (ST_FILTER*)st->data;
552:   PetscReal gamma;
553:   PetscBool isascii;

555:   PetscFunctionBegin;
556:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
557:   if (isascii) {
558:     PetscCall(PetscViewerASCIIPrintf(viewer,"  filter type: %s\n",STFilterTypes[ctx->type]));
559:     PetscCall(PetscViewerASCIIPrintf(viewer,"  interval of desired eigenvalues: [%g,%g]\n",(double)ctx->inta,(double)ctx->intb));
560:     PetscCall(PetscViewerASCIIPrintf(viewer,"  numerical range: [%g,%g]\n",(double)ctx->left,(double)ctx->right));
561:     PetscCall(PetscViewerASCIIPrintf(viewer,"  degree of filter polynomial: %" PetscInt_FMT "\n",ctx->polyDegree));
562:     if (ctx->damping && ctx->type==ST_FILTER_CHEBYSHEV) PetscCall(PetscViewerASCIIPrintf(viewer,"  type of damping = %s\n",STFilterDampings[ctx->damping]));
563:     if (st->state>=ST_STATE_SETUP && ctx->getthreshold) {
564:       PetscCall(STFilterGetThreshold(st,&gamma));
565:       PetscCall(PetscViewerASCIIPrintf(viewer,"  limit to accept eigenvalues: gamma=%g\n",(double)gamma));
566:     }
567:   }
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

571: static PetscErrorCode STDestroy_Filter(ST st)
572: {
573:   ST_FILTER *ctx = (ST_FILTER*)st->data;

575:   PetscFunctionBegin;
576:   if (ctx->destroy) PetscCall(ctx->destroy(st));
577:   PetscCall(PetscFree(st->data));
578:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetType_C",NULL));
579:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetType_C",NULL));
580:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",NULL));
581:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",NULL));
582:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",NULL));
583:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",NULL));
584:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",NULL));
585:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",NULL));
586:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",NULL));
587:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDamping_C",NULL));
588:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDamping_C",NULL));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: /*MC
593:    STFILTER - STFILTER = "filter" - A special type of `ST` that represents
594:    a polynomial filter.

596:    Level: beginner

598:    Notes:
599:    In standard eigenvalue problems, when the eigenvalues of interest are interior
600:    to the spectrum and we want to avoid the high cost associated with the matrix
601:    factorization of the shift-and-invert spectral transformation, an alternative
602:    is to build a high-order polynomial such that $p(A)$ enhances the wanted
603:    eigenvalues and filters out the unwanted ones.

605:    The definition of the filter is done with functions such as `STFilterSetDegree()`,
606:    `STFilterSetInterval()` or `STFilterSetType()`.

608: .seealso: [](ch:st), `ST`, `STType`, `STSetType()`, `STSetMatrices()`, `STFilterSetDegree()`, `STFilterSetInterval()`, `STFilterSetType()`
609: M*/

611: SLEPC_EXTERN PetscErrorCode STCreate_Filter(ST st)
612: {
613:   ST_FILTER   *ctx;

615:   PetscFunctionBegin;
616:   PetscCall(PetscNew(&ctx));
617:   st->data = (void*)ctx;

619:   st->usesksp = PETSC_FALSE;

621:   ctx->type               = (STFilterType)0;
622:   ctx->inta               = PETSC_MIN_REAL;
623:   ctx->intb               = PETSC_MAX_REAL;
624:   ctx->left               = 0.0;
625:   ctx->right              = 0.0;
626:   ctx->polyDegree         = 0;

628:   st->ops->apply           = STApply_Generic;
629:   st->ops->computeoperator = STComputeOperator_Filter;
630:   st->ops->setfromoptions  = STSetFromOptions_Filter;
631:   st->ops->destroy         = STDestroy_Filter;
632:   st->ops->reset           = STReset_Filter;
633:   st->ops->view            = STView_Filter;

635:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetType_C",STFilterSetType_Filter));
636:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetType_C",STFilterGetType_Filter));
637:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetInterval_C",STFilterSetInterval_Filter));
638:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetInterval_C",STFilterGetInterval_Filter));
639:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetRange_C",STFilterSetRange_Filter));
640:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetRange_C",STFilterGetRange_Filter));
641:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDegree_C",STFilterSetDegree_Filter));
642:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDegree_C",STFilterGetDegree_Filter));
643:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetThreshold_C",STFilterGetThreshold_Filter));
644:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterSetDamping_C",STFilterSetDamping_Filter));
645:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STFilterGetDamping_C",STFilterGetDamping_Filter));
646:   PetscFunctionReturn(PETSC_SUCCESS);
647: }