Actual source code: mfnbasic.c

slepc-3.17.2 2022-08-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:    Basic MFN routines
 12: */

 14: #include <slepc/private/mfnimpl.h>

 16: /* Logging support */
 17: PetscClassId      MFN_CLASSID = 0;
 18: PetscLogEvent     MFN_SetUp = 0,MFN_Solve = 0;

 20: /* List of registered MFN routines */
 21: PetscFunctionList MFNList = 0;
 22: PetscBool         MFNRegisterAllCalled = PETSC_FALSE;

 24: /* List of registered MFN monitors */
 25: PetscFunctionList MFNMonitorList              = NULL;
 26: PetscFunctionList MFNMonitorCreateList        = NULL;
 27: PetscFunctionList MFNMonitorDestroyList       = NULL;
 28: PetscBool         MFNMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@C
 31:    MFNView - Prints the MFN data structure.

 33:    Collective on mfn

 35:    Input Parameters:
 36: +  mfn - the matrix function solver context
 37: -  viewer - optional visualization context

 39:    Options Database Key:
 40: .  -mfn_view -  Calls MFNView() at end of MFNSolve()

 42:    Note:
 43:    The available visualization contexts include
 44: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 45: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 46:          output where only the first processor opens
 47:          the file.  All other processors send their
 48:          data to the first processor to print.

 50:    The user can open an alternative visualization context with
 51:    PetscViewerASCIIOpen() - output to a specified file.

 53:    Level: beginner

 55: .seealso: MFNCreate()
 56: @*/
 57: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
 58: {
 59:   PetscBool      isascii;

 62:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer);

 66:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 67:   if (isascii) {
 68:     PetscObjectPrintClassNamePrefixType((PetscObject)mfn,viewer);
 69:     if (mfn->ops->view) {
 70:       PetscViewerASCIIPushTab(viewer);
 71:       (*mfn->ops->view)(mfn,viewer);
 72:       PetscViewerASCIIPopTab(viewer);
 73:     }
 74:     PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",mfn->ncv);
 75:     PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",mfn->max_it);
 76:     PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)mfn->tol);
 77:   } else {
 78:     if (mfn->ops->view) (*mfn->ops->view)(mfn,viewer);
 79:   }
 80:   PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
 81:   if (!mfn->V) MFNGetFN(mfn,&mfn->fn);
 82:   FNView(mfn->fn,viewer);
 83:   if (!mfn->V) MFNGetBV(mfn,&mfn->V);
 84:   BVView(mfn->V,viewer);
 85:   PetscViewerPopFormat(viewer);
 86:   PetscFunctionReturn(0);
 87: }

 89: /*@C
 90:    MFNViewFromOptions - View from options

 92:    Collective on MFN

 94:    Input Parameters:
 95: +  mfn  - the matrix function context
 96: .  obj  - optional object
 97: -  name - command line option

 99:    Level: intermediate

101: .seealso: MFNView(), MFNCreate()
102: @*/
103: PetscErrorCode MFNViewFromOptions(MFN mfn,PetscObject obj,const char name[])
104: {
106:   PetscObjectViewFromOptions((PetscObject)mfn,obj,name);
107:   PetscFunctionReturn(0);
108: }
109: /*@C
110:    MFNConvergedReasonView - Displays the reason an MFN solve converged or diverged.

112:    Collective on mfn

114:    Input Parameters:
115: +  mfn - the matrix function context
116: -  viewer - the viewer to display the reason

118:    Options Database Keys:
119: .  -mfn_converged_reason - print reason for convergence, and number of iterations

121:    Note:
122:    To change the format of the output call PetscViewerPushFormat(viewer,format) before
123:    this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
124:    display a reason if it fails. The latter can be set in the command line with
125:    -mfn_converged_reason ::failed

127:    Level: intermediate

129: .seealso: MFNSetTolerances(), MFNGetIterationNumber(), MFNConvergedReasonViewFromOptions()
130: @*/
131: PetscErrorCode MFNConvergedReasonView(MFN mfn,PetscViewer viewer)
132: {
133:   PetscBool         isAscii;
134:   PetscViewerFormat format;

136:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)mfn));
137:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
138:   if (isAscii) {
139:     PetscViewerGetFormat(viewer,&format);
140:     PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel);
141:     if (mfn->reason > 0 && format != PETSC_VIEWER_FAILED) PetscViewerASCIIPrintf(viewer,"%s Matrix function solve converged due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
142:     else if (mfn->reason <= 0) PetscViewerASCIIPrintf(viewer,"%s Matrix function solve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
143:     PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel);
144:   }
145:   PetscFunctionReturn(0);
146: }

148: /*@
149:    MFNConvergedReasonViewFromOptions - Processes command line options to determine if/how
150:    the MFN converged reason is to be viewed.

152:    Collective on mfn

154:    Input Parameter:
155: .  mfn - the matrix function context

157:    Level: developer

159: .seealso: MFNConvergedReasonView()
160: @*/
161: PetscErrorCode MFNConvergedReasonViewFromOptions(MFN mfn)
162: {
163:   PetscViewer       viewer;
164:   PetscBool         flg;
165:   static PetscBool  incall = PETSC_FALSE;
166:   PetscViewerFormat format;

168:   if (incall) PetscFunctionReturn(0);
169:   incall = PETSC_TRUE;
170:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg);
171:   if (flg) {
172:     PetscViewerPushFormat(viewer,format);
173:     MFNConvergedReasonView(mfn,viewer);
174:     PetscViewerPopFormat(viewer);
175:     PetscViewerDestroy(&viewer);
176:   }
177:   incall = PETSC_FALSE;
178:   PetscFunctionReturn(0);
179: }

181: /*@
182:    MFNCreate - Creates the default MFN context.

184:    Collective

186:    Input Parameter:
187: .  comm - MPI communicator

189:    Output Parameter:
190: .  outmfn - location to put the MFN context

192:    Note:
193:    The default MFN type is MFNKRYLOV

195:    Level: beginner

197: .seealso: MFNSetUp(), MFNSolve(), MFNDestroy(), MFN
198: @*/
199: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
200: {
201:   MFN            mfn;

204:   *outmfn = 0;
205:   MFNInitializePackage();
206:   SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView);

208:   mfn->A               = NULL;
209:   mfn->fn              = NULL;
210:   mfn->max_it          = PETSC_DEFAULT;
211:   mfn->ncv             = PETSC_DEFAULT;
212:   mfn->tol             = PETSC_DEFAULT;
213:   mfn->errorifnotconverged = PETSC_FALSE;

215:   mfn->numbermonitors  = 0;

217:   mfn->V               = NULL;
218:   mfn->nwork           = 0;
219:   mfn->work            = NULL;
220:   mfn->data            = NULL;

222:   mfn->its             = 0;
223:   mfn->nv              = 0;
224:   mfn->errest          = 0;
225:   mfn->setupcalled     = 0;
226:   mfn->reason          = MFN_CONVERGED_ITERATING;

228:   *outmfn = mfn;
229:   PetscFunctionReturn(0);
230: }

232: /*@C
233:    MFNSetType - Selects the particular solver to be used in the MFN object.

235:    Logically Collective on mfn

237:    Input Parameters:
238: +  mfn  - the matrix function context
239: -  type - a known method

241:    Options Database Key:
242: .  -mfn_type <method> - Sets the method; use -help for a list
243:     of available methods

245:    Notes:
246:    See "slepc/include/slepcmfn.h" for available methods. The default
247:    is MFNKRYLOV

249:    Normally, it is best to use the MFNSetFromOptions() command and
250:    then set the MFN type from the options database rather than by using
251:    this routine.  Using the options database provides the user with
252:    maximum flexibility in evaluating the different available methods.
253:    The MFNSetType() routine is provided for those situations where it
254:    is necessary to set the iterative solver independently of the command
255:    line or options database.

257:    Level: intermediate

259: .seealso: MFNType
260: @*/
261: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
262: {
263:   PetscErrorCode (*r)(MFN);
264:   PetscBool      match;


269:   PetscObjectTypeCompare((PetscObject)mfn,type,&match);
270:   if (match) PetscFunctionReturn(0);

272:   PetscFunctionListFind(MFNList,type,&r);

275:   if (mfn->ops->destroy) (*mfn->ops->destroy)(mfn);
276:   PetscMemzero(mfn->ops,sizeof(struct _MFNOps));

278:   mfn->setupcalled = 0;
279:   PetscObjectChangeTypeName((PetscObject)mfn,type);
280:   (*r)(mfn);
281:   PetscFunctionReturn(0);
282: }

284: /*@C
285:    MFNGetType - Gets the MFN type as a string from the MFN object.

287:    Not Collective

289:    Input Parameter:
290: .  mfn - the matrix function context

292:    Output Parameter:
293: .  type - name of MFN method

295:    Level: intermediate

297: .seealso: MFNSetType()
298: @*/
299: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
300: {
303:   *type = ((PetscObject)mfn)->type_name;
304:   PetscFunctionReturn(0);
305: }

307: /*@C
308:    MFNRegister - Adds a method to the matrix function solver package.

310:    Not Collective

312:    Input Parameters:
313: +  name - name of a new user-defined solver
314: -  function - routine to create the solver context

316:    Notes:
317:    MFNRegister() may be called multiple times to add several user-defined solvers.

319:    Sample usage:
320: .vb
321:     MFNRegister("my_solver",MySolverCreate);
322: .ve

324:    Then, your solver can be chosen with the procedural interface via
325: $     MFNSetType(mfn,"my_solver")
326:    or at runtime via the option
327: $     -mfn_type my_solver

329:    Level: advanced

331: .seealso: MFNRegisterAll()
332: @*/
333: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
334: {
335:   MFNInitializePackage();
336:   PetscFunctionListAdd(&MFNList,name,function);
337:   PetscFunctionReturn(0);
338: }

340: /*@C
341:    MFNMonitorRegister - Adds MFN monitor routine.

343:    Not Collective

345:    Input Parameters:
346: +  name    - name of a new monitor routine
347: .  vtype   - a PetscViewerType for the output
348: .  format  - a PetscViewerFormat for the output
349: .  monitor - monitor routine
350: .  create  - creation routine, or NULL
351: -  destroy - destruction routine, or NULL

353:    Notes:
354:    MFNMonitorRegister() may be called multiple times to add several user-defined monitors.

356:    Sample usage:
357: .vb
358:    MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
359: .ve

361:    Then, your monitor can be chosen with the procedural interface via
362: $      MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL)
363:    or at runtime via the option
364: $      -mfn_monitor_my_monitor

366:    Level: advanced

368: .seealso: MFNMonitorRegisterAll()
369: @*/
370: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(MFN,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
371: {
372:   char           key[PETSC_MAX_PATH_LEN];

374:   MFNInitializePackage();
375:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
376:   PetscFunctionListAdd(&MFNMonitorList,key,monitor);
377:   if (create)  PetscFunctionListAdd(&MFNMonitorCreateList,key,create);
378:   if (destroy) PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy);
379:   PetscFunctionReturn(0);
380: }

382: /*@
383:    MFNReset - Resets the MFN context to the initial state (prior to setup)
384:    and destroys any allocated Vecs and Mats.

386:    Collective on mfn

388:    Input Parameter:
389: .  mfn - matrix function context obtained from MFNCreate()

391:    Level: advanced

393: .seealso: MFNDestroy()
394: @*/
395: PetscErrorCode MFNReset(MFN mfn)
396: {
398:   if (!mfn) PetscFunctionReturn(0);
399:   if (mfn->ops->reset) (mfn->ops->reset)(mfn);
400:   MatDestroy(&mfn->A);
401:   BVDestroy(&mfn->V);
402:   VecDestroyVecs(mfn->nwork,&mfn->work);
403:   mfn->nwork = 0;
404:   mfn->setupcalled = 0;
405:   PetscFunctionReturn(0);
406: }

408: /*@C
409:    MFNDestroy - Destroys the MFN context.

411:    Collective on mfn

413:    Input Parameter:
414: .  mfn - matrix function context obtained from MFNCreate()

416:    Level: beginner

418: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
419: @*/
420: PetscErrorCode MFNDestroy(MFN *mfn)
421: {
422:   if (!*mfn) PetscFunctionReturn(0);
424:   if (--((PetscObject)(*mfn))->refct > 0) { *mfn = 0; PetscFunctionReturn(0); }
425:   MFNReset(*mfn);
426:   if ((*mfn)->ops->destroy) (*(*mfn)->ops->destroy)(*mfn);
427:   FNDestroy(&(*mfn)->fn);
428:   MatDestroy(&(*mfn)->AT);
429:   MFNMonitorCancel(*mfn);
430:   PetscHeaderDestroy(mfn);
431:   PetscFunctionReturn(0);
432: }

434: /*@
435:    MFNSetBV - Associates a basis vectors object to the matrix function solver.

437:    Collective on mfn

439:    Input Parameters:
440: +  mfn - matrix function context obtained from MFNCreate()
441: -  bv  - the basis vectors object

443:    Note:
444:    Use MFNGetBV() to retrieve the basis vectors context (for example,
445:    to free it at the end of the computations).

447:    Level: advanced

449: .seealso: MFNGetBV()
450: @*/
451: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
452: {
456:   PetscObjectReference((PetscObject)bv);
457:   BVDestroy(&mfn->V);
458:   mfn->V = bv;
459:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
460:   PetscFunctionReturn(0);
461: }

463: /*@
464:    MFNGetBV - Obtain the basis vectors object associated to the matrix
465:    function solver.

467:    Not Collective

469:    Input Parameters:
470: .  mfn - matrix function context obtained from MFNCreate()

472:    Output Parameter:
473: .  bv - basis vectors context

475:    Level: advanced

477: .seealso: MFNSetBV()
478: @*/
479: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
480: {
483:   if (!mfn->V) {
484:     BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V);
485:     PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0);
486:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
487:     PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options);
488:   }
489:   *bv = mfn->V;
490:   PetscFunctionReturn(0);
491: }

493: /*@
494:    MFNSetFN - Specifies the function to be computed.

496:    Collective on mfn

498:    Input Parameters:
499: +  mfn - matrix function context obtained from MFNCreate()
500: -  fn  - the math function object

502:    Note:
503:    Use MFNGetFN() to retrieve the math function context (for example,
504:    to free it at the end of the computations).

506:    Level: beginner

508: .seealso: MFNGetFN()
509: @*/
510: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
511: {
515:   PetscObjectReference((PetscObject)fn);
516:   FNDestroy(&mfn->fn);
517:   mfn->fn = fn;
518:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
519:   PetscFunctionReturn(0);
520: }

522: /*@
523:    MFNGetFN - Obtain the math function object associated to the MFN object.

525:    Not Collective

527:    Input Parameters:
528: .  mfn - matrix function context obtained from MFNCreate()

530:    Output Parameter:
531: .  fn - math function context

533:    Level: beginner

535: .seealso: MFNSetFN()
536: @*/
537: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
538: {
541:   if (!mfn->fn) {
542:     FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn);
543:     PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0);
544:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
545:     PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options);
546:   }
547:   *fn = mfn->fn;
548:   PetscFunctionReturn(0);
549: }