Actual source code: mfnbasic.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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
 54: @*/
 55: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
 56: {
 58:   PetscBool      isascii;

 62:   if (!viewer) {
 63:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer);
 64:   }

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

 93: /*@C
 94:    MFNViewFromOptions - View from options

 96:    Collective on MFN

 98:    Input Parameters:
 99: +  mfn  - the matrix function context
100: .  obj  - optional object
101: -  name - command line option

103:    Level: intermediate

105: .seealso: MFNView(), MFNCreate()
106: @*/
107: PetscErrorCode MFNViewFromOptions(MFN mfn,PetscObject obj,const char name[])
108: {

113:   PetscObjectViewFromOptions((PetscObject)mfn,obj,name);
114:   return(0);
115: }
116: /*@C
117:    MFNConvergedReasonView - Displays the reason an MFN solve converged or diverged.

119:    Collective on mfn

121:    Input Parameters:
122: +  mfn - the matrix function context
123: -  viewer - the viewer to display the reason

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

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

134:    Level: intermediate

136: .seealso: MFNSetTolerances(), MFNGetIterationNumber(), MFNConvergedReasonViewFromOptions()
137: @*/
138: PetscErrorCode MFNConvergedReasonView(MFN mfn,PetscViewer viewer)
139: {
140:   PetscErrorCode    ierr;
141:   PetscBool         isAscii;
142:   PetscViewerFormat format;

145:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)mfn));
146:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
147:   if (isAscii) {
148:     PetscViewerGetFormat(viewer,&format);
149:     PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel);
150:     if (mfn->reason > 0 && format != PETSC_VIEWER_FAILED) {
151:       PetscViewerASCIIPrintf(viewer,"%s Matrix function solve converged due to %s; iterations %D\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
152:     } else if (mfn->reason <= 0) {
153:       PetscViewerASCIIPrintf(viewer,"%s Matrix function solve did not converge due to %s; iterations %D\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
154:     }
155:     PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel);
156:   }
157:   return(0);
158: }

160: /*@
161:    MFNConvergedReasonViewFromOptions - Processes command line options to determine if/how
162:    the MFN converged reason is to be viewed.

164:    Collective on mfn

166:    Input Parameter:
167: .  mfn - the matrix function context

169:    Level: developer

171: .seealso: MFNConvergedReasonView()
172: @*/
173: PetscErrorCode MFNConvergedReasonViewFromOptions(MFN mfn)
174: {
175:   PetscErrorCode    ierr;
176:   PetscViewer       viewer;
177:   PetscBool         flg;
178:   static PetscBool  incall = PETSC_FALSE;
179:   PetscViewerFormat format;

182:   if (incall) return(0);
183:   incall = PETSC_TRUE;
184:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg);
185:   if (flg) {
186:     PetscViewerPushFormat(viewer,format);
187:     MFNConvergedReasonView(mfn,viewer);
188:     PetscViewerPopFormat(viewer);
189:     PetscViewerDestroy(&viewer);
190:   }
191:   incall = PETSC_FALSE;
192:   return(0);
193: }

195: /*@
196:    MFNCreate - Creates the default MFN context.

198:    Collective

200:    Input Parameter:
201: .  comm - MPI communicator

203:    Output Parameter:
204: .  mfn - location to put the MFN context

206:    Note:
207:    The default MFN type is MFNKRYLOV

209:    Level: beginner

211: .seealso: MFNSetUp(), MFNSolve(), MFNDestroy(), MFN
212: @*/
213: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
214: {
216:   MFN            mfn;

220:   *outmfn = 0;
221:   MFNInitializePackage();
222:   SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView);

224:   mfn->A               = NULL;
225:   mfn->fn              = NULL;
226:   mfn->max_it          = PETSC_DEFAULT;
227:   mfn->ncv             = PETSC_DEFAULT;
228:   mfn->tol             = PETSC_DEFAULT;
229:   mfn->errorifnotconverged = PETSC_FALSE;

231:   mfn->numbermonitors  = 0;

233:   mfn->V               = NULL;
234:   mfn->nwork           = 0;
235:   mfn->work            = NULL;
236:   mfn->data            = NULL;

238:   mfn->its             = 0;
239:   mfn->nv              = 0;
240:   mfn->errest          = 0;
241:   mfn->setupcalled     = 0;
242:   mfn->reason          = MFN_CONVERGED_ITERATING;

244:   *outmfn = mfn;
245:   return(0);
246: }

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

251:    Logically Collective on mfn

253:    Input Parameters:
254: +  mfn  - the matrix function context
255: -  type - a known method

257:    Options Database Key:
258: .  -mfn_type <method> - Sets the method; use -help for a list
259:     of available methods

261:    Notes:
262:    See "slepc/include/slepcmfn.h" for available methods. The default
263:    is MFNKRYLOV

265:    Normally, it is best to use the MFNSetFromOptions() command and
266:    then set the MFN type from the options database rather than by using
267:    this routine.  Using the options database provides the user with
268:    maximum flexibility in evaluating the different available methods.
269:    The MFNSetType() routine is provided for those situations where it
270:    is necessary to set the iterative solver independently of the command
271:    line or options database.

273:    Level: intermediate

275: .seealso: MFNType
276: @*/
277: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
278: {
279:   PetscErrorCode ierr,(*r)(MFN);
280:   PetscBool      match;


286:   PetscObjectTypeCompare((PetscObject)mfn,type,&match);
287:   if (match) return(0);

289:   PetscFunctionListFind(MFNList,type,&r);
290:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MFN type given: %s",type);

292:   if (mfn->ops->destroy) { (*mfn->ops->destroy)(mfn); }
293:   PetscMemzero(mfn->ops,sizeof(struct _MFNOps));

295:   mfn->setupcalled = 0;
296:   PetscObjectChangeTypeName((PetscObject)mfn,type);
297:   (*r)(mfn);
298:   return(0);
299: }

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

304:    Not Collective

306:    Input Parameter:
307: .  mfn - the matrix function context

309:    Output Parameter:
310: .  name - name of MFN method

312:    Level: intermediate

314: .seealso: MFNSetType()
315: @*/
316: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
317: {
321:   *type = ((PetscObject)mfn)->type_name;
322:   return(0);
323: }

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

328:    Not Collective

330:    Input Parameters:
331: +  name - name of a new user-defined solver
332: -  function - routine to create the solver context

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

337:    Sample usage:
338: .vb
339:     MFNRegister("my_solver",MySolverCreate);
340: .ve

342:    Then, your solver can be chosen with the procedural interface via
343: $     MFNSetType(mfn,"my_solver")
344:    or at runtime via the option
345: $     -mfn_type my_solver

347:    Level: advanced

349: .seealso: MFNRegisterAll()
350: @*/
351: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
352: {

356:   MFNInitializePackage();
357:   PetscFunctionListAdd(&MFNList,name,function);
358:   return(0);
359: }

361: /*@C
362:    MFNMonitorRegister - Adds MFN monitor routine.

364:    Not Collective

366:    Input Parameters:
367: +  name    - name of a new monitor routine
368: .  vtype   - a PetscViewerType for the output
369: .  format  - a PetscViewerFormat for the output
370: .  monitor - monitor routine
371: .  create  - creation routine, or NULL
372: -  destroy - destruction routine, or NULL

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

377:    Sample usage:
378: .vb
379:    MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
380: .ve

382:    Then, your monitor can be chosen with the procedural interface via
383: $      MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL)
384:    or at runtime via the option
385: $      -mfn_monitor_my_monitor

387:    Level: advanced

389: .seealso: MFNMonitorRegisterAll()
390: @*/
391: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(MFN,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
392: {
393:   char           key[PETSC_MAX_PATH_LEN];

397:   MFNInitializePackage();
398:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
399:   PetscFunctionListAdd(&MFNMonitorList,key,monitor);
400:   if (create)  { PetscFunctionListAdd(&MFNMonitorCreateList,key,create); }
401:   if (destroy) { PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy); }
402:   return(0);
403: }

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

409:    Collective on mfn

411:    Input Parameter:
412: .  mfn - matrix function context obtained from MFNCreate()

414:    Level: advanced

416: .seealso: MFNDestroy()
417: @*/
418: PetscErrorCode MFNReset(MFN mfn)
419: {

424:   if (!mfn) return(0);
425:   if (mfn->ops->reset) { (mfn->ops->reset)(mfn); }
426:   MatDestroy(&mfn->A);
427:   BVDestroy(&mfn->V);
428:   VecDestroyVecs(mfn->nwork,&mfn->work);
429:   mfn->nwork = 0;
430:   mfn->setupcalled = 0;
431:   return(0);
432: }

434: /*@C
435:    MFNDestroy - Destroys the MFN context.

437:    Collective on mfn

439:    Input Parameter:
440: .  mfn - matrix function context obtained from MFNCreate()

442:    Level: beginner

444: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
445: @*/
446: PetscErrorCode MFNDestroy(MFN *mfn)
447: {

451:   if (!*mfn) return(0);
453:   if (--((PetscObject)(*mfn))->refct > 0) { *mfn = 0; return(0); }
454:   MFNReset(*mfn);
455:   if ((*mfn)->ops->destroy) { (*(*mfn)->ops->destroy)(*mfn); }
456:   FNDestroy(&(*mfn)->fn);
457:   MatDestroy(&(*mfn)->AT);
458:   MFNMonitorCancel(*mfn);
459:   PetscHeaderDestroy(mfn);
460:   return(0);
461: }

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

466:    Collective on mfn

468:    Input Parameters:
469: +  mfn - matrix function context obtained from MFNCreate()
470: -  bv  - the basis vectors object

472:    Note:
473:    Use MFNGetBV() to retrieve the basis vectors context (for example,
474:    to free it at the end of the computations).

476:    Level: advanced

478: .seealso: MFNGetBV()
479: @*/
480: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
481: {

488:   PetscObjectReference((PetscObject)bv);
489:   BVDestroy(&mfn->V);
490:   mfn->V = bv;
491:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
492:   return(0);
493: }

495: /*@
496:    MFNGetBV - Obtain the basis vectors object associated to the matrix
497:    function solver.

499:    Not Collective

501:    Input Parameters:
502: .  mfn - matrix function context obtained from MFNCreate()

504:    Output Parameter:
505: .  bv - basis vectors context

507:    Level: advanced

509: .seealso: MFNSetBV()
510: @*/
511: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
512: {

518:   if (!mfn->V) {
519:     BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V);
520:     PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0);
521:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
522:     PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options);
523:   }
524:   *bv = mfn->V;
525:   return(0);
526: }

528: /*@
529:    MFNSetFN - Specifies the function to be computed.

531:    Collective on mfn

533:    Input Parameters:
534: +  mfn - matrix function context obtained from MFNCreate()
535: -  fn  - the math function object

537:    Note:
538:    Use MFNGetFN() to retrieve the math function context (for example,
539:    to free it at the end of the computations).

541:    Level: beginner

543: .seealso: MFNGetFN()
544: @*/
545: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
546: {

553:   PetscObjectReference((PetscObject)fn);
554:   FNDestroy(&mfn->fn);
555:   mfn->fn = fn;
556:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
557:   return(0);
558: }

560: /*@
561:    MFNGetFN - Obtain the math function object associated to the MFN object.

563:    Not Collective

565:    Input Parameters:
566: .  mfn - matrix function context obtained from MFNCreate()

568:    Output Parameter:
569: .  fn - math function context

571:    Level: beginner

573: .seealso: MFNSetFN()
574: @*/
575: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
576: {

582:   if (!mfn->fn) {
583:     FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn);
584:     PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0);
585:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
586:     PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options);
587:   }
588:   *fn = mfn->fn;
589:   return(0);
590: }