Actual source code: mfnbasic.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:    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 = NULL;
 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: /*@
 31:    MFNView - Prints the `MFN` data structure.

 33:    Collective

 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:    Notes:
 43:    The available visualization contexts include
 44: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
 45: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
 46:          first process opens the file; all other processes send their data to the
 47:          first one to print

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

 52:    Level: beginner

 54: .seealso: [](ch:mfn), `MFNCreate()`, `MFNViewFromOptions()`
 55: @*/
 56: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
 57: {
 58:   PetscBool      isascii;

 60:   PetscFunctionBegin;
 62:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer));
 64:   PetscCheckSameComm(mfn,1,viewer,2);

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

 85: /*@
 86:    MFNViewFromOptions - View (print) an `MFN` object based on values in the options database.

 88:    Collective

 90:    Input Parameters:
 91: +  mfn  - the matrix function solver context
 92: .  obj  - optional object that provides the options prefix used to query the options database
 93: -  name - command line option

 95:    Level: intermediate

 97: .seealso: [](ch:mfn), `MFNView()`, `MFNCreate()`, `PetscObjectViewFromOptions()`
 98: @*/
 99: PetscErrorCode MFNViewFromOptions(MFN mfn,PetscObject obj,const char name[])
100: {
101:   PetscFunctionBegin;
103:   PetscCall(PetscObjectViewFromOptions((PetscObject)mfn,obj,name));
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }
106: /*@
107:    MFNConvergedReasonView - Displays the reason an `MFN` solve converged or diverged.

109:    Collective

111:    Input Parameters:
112: +  mfn - the matrix function solver context
113: -  viewer - the viewer to display the reason

115:    Options Database Key:
116: .  -mfn_converged_reason - print reason for convergence/divergence, and number of iterations

118:    Notes:
119:    Use `MFNConvergedReasonViewFromOptions()` to display the reason based on values
120:    in the options database.

122:    To change the format of the output call `PetscViewerPushFormat()` before this
123:    call. Use `PETSC_VIEWER_DEFAULT` for the default, or `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: [](ch:mfn), `MFNSetTolerances()`, `MFNGetIterationNumber()`, `MFNConvergedReasonViewFromOptions()`
130: @*/
131: PetscErrorCode MFNConvergedReasonView(MFN mfn,PetscViewer viewer)
132: {
133:   PetscBool         isAscii;
134:   PetscViewerFormat format;

136:   PetscFunctionBegin;
137:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)mfn));
138:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
139:   if (isAscii) {
140:     PetscCall(PetscViewerGetFormat(viewer,&format));
141:     PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel));
142:     if (mfn->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(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));
143:     else if (mfn->reason <= 0) PetscCall(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));
144:     PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel));
145:   }
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }

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

153:    Collective

155:    Input Parameter:
156: .  mfn - the matrix function solver context

158:    Level: intermediate

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

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

183: /*@
184:    MFNCreate - Creates the default `MFN` context.

186:    Collective

188:    Input Parameter:
189: .  comm - MPI communicator

191:    Output Parameter:
192: .  outmfn - location to put the `MFN` context

194:    Note:
195:    The default `MFN` type is `MFNKRYLOV`.

197:    Level: beginner

199: .seealso: [](ch:mfn), `MFNSetUp()`, `MFNSolve()`, `MFNDestroy()`, `MFN`
200: @*/
201: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
202: {
203:   MFN            mfn;

205:   PetscFunctionBegin;
206:   PetscAssertPointer(outmfn,2);
207:   PetscCall(MFNInitializePackage());
208:   PetscCall(SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView));

210:   mfn->A               = NULL;
211:   mfn->fn              = NULL;
212:   mfn->max_it          = PETSC_DETERMINE;
213:   mfn->ncv             = PETSC_DETERMINE;
214:   mfn->tol             = PETSC_DETERMINE;
215:   mfn->errorifnotconverged = PETSC_FALSE;

217:   mfn->numbermonitors  = 0;

219:   mfn->V               = NULL;
220:   mfn->nwork           = 0;
221:   mfn->work            = NULL;
222:   mfn->data            = NULL;

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

230:   *outmfn = mfn;
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: /*@
235:    MFNSetType - Selects the particular solver to be used in the `MFN` object.

237:    Logically Collective

239:    Input Parameters:
240: +  mfn  - the matrix function solver context
241: -  type - a known method

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

247:    Notes:
248:    See `MFNType` for available methods. The default is `MFNKRYLOV`.

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

258:    Level: intermediate

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

267:   PetscFunctionBegin;
269:   PetscAssertPointer(type,2);

271:   PetscCall(PetscObjectTypeCompare((PetscObject)mfn,type,&match));
272:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

274:   PetscCall(PetscFunctionListFind(MFNList,type,&r));
275:   PetscCheck(r,PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MFN type given: %s",type);

277:   PetscTryTypeMethod(mfn,destroy);
278:   PetscCall(PetscMemzero(mfn->ops,sizeof(struct _MFNOps)));

280:   mfn->setupcalled = 0;
281:   PetscCall(PetscObjectChangeTypeName((PetscObject)mfn,type));
282:   PetscCall((*r)(mfn));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: /*@
287:    MFNGetType - Gets the `MFN` type as a string from the `MFN` object.

289:    Not Collective

291:    Input Parameter:
292: .  mfn - the matrix function solver context

294:    Output Parameter:
295: .  type - name of `MFN` method

297:    Level: intermediate

299: .seealso: [](ch:mfn), `MFNSetType()`
300: @*/
301: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
302: {
303:   PetscFunctionBegin;
305:   PetscAssertPointer(type,2);
306:   *type = ((PetscObject)mfn)->type_name;
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

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

313:    Not Collective

315:    Input Parameters:
316: +  name - name of a new user-defined solver
317: -  function - routine to create the solver context

319:    Note:
320:    `MFNRegister()` may be called multiple times to add several user-defined solvers.

322:    Example Usage:
323: .vb
324:    MFNRegister("my_solver",MySolverCreate);
325: .ve

327:    Then, your solver can be chosen with the procedural interface via
328: .vb
329:    MFNSetType(mfn,"my_solver")
330: .ve
331:    or at runtime via the option `-mfn_type my_solver`.

333:    Level: advanced

335: .seealso: [](ch:mfn), `MFNRegisterAll()`
336: @*/
337: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
338: {
339:   PetscFunctionBegin;
340:   PetscCall(MFNInitializePackage());
341:   PetscCall(PetscFunctionListAdd(&MFNList,name,function));
342:   PetscFunctionReturn(PETSC_SUCCESS);
343: }

345: /*@C
346:    MFNMonitorRegister - Registers an `MFN` monitor routine that may be accessed with
347:    `MFNMonitorSetFromOptions()`.

349:    Not Collective

351:    Input Parameters:
352: +  name    - name of a new monitor routine
353: .  vtype   - a `PetscViewerType` for the output
354: .  format  - a `PetscViewerFormat` for the output
355: .  monitor - monitor routine, see `MFNMonitorRegisterFn`
356: .  create  - creation routine, or `NULL`
357: -  destroy - destruction routine, or `NULL`

359:    Notes:
360:    `MFNMonitorRegister()` may be called multiple times to add several user-defined monitors.

362:    The calling sequence for the given function matches the calling sequence of `MFNMonitorFn`
363:    functions passed to `MFNMonitorSet()` with the additional requirement that its final argument
364:    be a `PetscViewerAndFormat`.

366:    Example Usage:
367: .vb
368:    MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
369: .ve

371:    Then, your monitor can be chosen with the procedural interface via
372: .vb
373:    MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL);
374: .ve
375:    or at runtime via the option `-mfn_monitor_my_monitor`.

377:    Level: advanced

379: .seealso: [](ch:mfn), `MFNMonitorSet()`, `MFNMonitorRegisterAll()`, `MFNMonitorSetFromOptions()`
380: @*/
381: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,MFNMonitorRegisterFn *monitor,MFNMonitorRegisterCreateFn *create,MFNMonitorRegisterDestroyFn *destroy)
382: {
383:   char           key[PETSC_MAX_PATH_LEN];

385:   PetscFunctionBegin;
386:   PetscCall(MFNInitializePackage());
387:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
388:   PetscCall(PetscFunctionListAdd(&MFNMonitorList,key,monitor));
389:   if (create)  PetscCall(PetscFunctionListAdd(&MFNMonitorCreateList,key,create));
390:   if (destroy) PetscCall(PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy));
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

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

398:    Collective

400:    Input Parameter:
401: .  mfn - the matrix function solver context

403:    Level: advanced

405: .seealso: [](ch:mfn), `MFNDestroy()`
406: @*/
407: PetscErrorCode MFNReset(MFN mfn)
408: {
409:   PetscFunctionBegin;
411:   if (!mfn) PetscFunctionReturn(PETSC_SUCCESS);
412:   PetscTryTypeMethod(mfn,reset);
413:   PetscCall(MatDestroy(&mfn->A));
414:   PetscCall(BVDestroy(&mfn->V));
415:   PetscCall(VecDestroyVecs(mfn->nwork,&mfn->work));
416:   mfn->nwork = 0;
417:   mfn->setupcalled = 0;
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: /*@
422:    MFNDestroy - Destroys the MFN context.

424:    Collective

426:    Input Parameter:
427: .  mfn - the matrix function solver context

429:    Level: beginner

431: .seealso: [](ch:mfn), `MFNCreate()`, `MFNSetUp()`, `MFNSolve()`
432: @*/
433: PetscErrorCode MFNDestroy(MFN *mfn)
434: {
435:   PetscFunctionBegin;
436:   if (!*mfn) PetscFunctionReturn(PETSC_SUCCESS);
438:   if (--((PetscObject)*mfn)->refct > 0) { *mfn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
439:   PetscCall(MFNReset(*mfn));
440:   PetscTryTypeMethod(*mfn,destroy);
441:   PetscCall(FNDestroy(&(*mfn)->fn));
442:   PetscCall(MatDestroy(&(*mfn)->AT));
443:   PetscCall(MFNMonitorCancel(*mfn));
444:   PetscCall(PetscHeaderDestroy(mfn));
445:   PetscFunctionReturn(PETSC_SUCCESS);
446: }

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

451:    Collective

453:    Input Parameters:
454: +  mfn - the matrix function solver context
455: -  bv  - the basis vectors object

457:    Note:
458:    Use MFNGetBV() to retrieve the basis vectors context (for example,
459:    to free it at the end of the computations).

461:    Level: advanced

463: .seealso: [](ch:mfn), `MFNGetBV()`
464: @*/
465: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
466: {
467:   PetscFunctionBegin;
470:   PetscCheckSameComm(mfn,1,bv,2);
471:   PetscCall(PetscObjectReference((PetscObject)bv));
472:   PetscCall(BVDestroy(&mfn->V));
473:   mfn->V = bv;
474:   PetscFunctionReturn(PETSC_SUCCESS);
475: }

477: /*@
478:    MFNGetBV - Obtain the basis vectors object associated to the matrix
479:    function solver.

481:    Not Collective

483:    Input Parameter:
484: .  mfn - the matrix function solver context

486:    Output Parameter:
487: .  bv - basis vectors context

489:    Level: advanced

491: .seealso: [](ch:mfn), `MFNSetBV()`
492: @*/
493: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
494: {
495:   PetscFunctionBegin;
497:   PetscAssertPointer(bv,2);
498:   if (!mfn->V) {
499:     PetscCall(BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V));
500:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0));
501:     PetscCall(PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options));
502:   }
503:   *bv = mfn->V;
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

507: /*@
508:    MFNSetFN - Specifies the function to be computed.

510:    Collective

512:    Input Parameters:
513: +  mfn - the matrix function solver context
514: -  fn  - the math function object

516:    Note:
517:    Use MFNGetFN() to retrieve the math function context (for example,
518:    to free it at the end of the computations).

520:    Level: beginner

522: .seealso: [](ch:mfn), `MFNGetFN()`
523: @*/
524: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
525: {
526:   PetscFunctionBegin;
529:   PetscCheckSameComm(mfn,1,fn,2);
530:   PetscCall(PetscObjectReference((PetscObject)fn));
531:   PetscCall(FNDestroy(&mfn->fn));
532:   mfn->fn = fn;
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: /*@
537:    MFNGetFN - Obtain the math function object associated to the MFN object.

539:    Not Collective

541:    Input Parameter:
542: .  mfn - the matrix function solver context

544:    Output Parameter:
545: .  fn - math function context

547:    Level: beginner

549: .seealso: [](ch:mfn), `MFNSetFN()`
550: @*/
551: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
552: {
553:   PetscFunctionBegin;
555:   PetscAssertPointer(fn,2);
556:   if (!mfn->fn) {
557:     PetscCall(FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn));
558:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0));
559:     PetscCall(PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options));
560:   }
561:   *fn = mfn->fn;
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }