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:    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;

 61:   PetscFunctionBegin;
 63:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer));
 65:   PetscCheckSameComm(mfn,1,viewer,2);

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

 86: /*@
 87:    MFNViewFromOptions - View from options

 89:    Collective

 91:    Input Parameters:
 92: +  mfn  - the matrix function context
 93: .  obj  - optional object
 94: -  name - command line option

 96:    Level: intermediate

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

110:    Collective

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

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

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

125:    Level: intermediate

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

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

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

151:    Collective

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

156:    Level: developer

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

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

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

208:   mfn->A               = NULL;
209:   mfn->fn              = NULL;
210:   mfn->max_it          = PETSC_DETERMINE;
211:   mfn->ncv             = PETSC_DETERMINE;
212:   mfn->tol             = PETSC_DETERMINE;
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(PETSC_SUCCESS);
230: }

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

235:    Logically Collective

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;

266:   PetscFunctionBegin;
268:   PetscAssertPointer(type,2);

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

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

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

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

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

288:    Not Collective

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

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

296:    Level: intermediate

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

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

312:    Not Collective

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

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

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

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

331:    Level: advanced

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

343: /*@C
344:    MFNMonitorRegister - Registers an MFN monitor routine that may be accessed with MFNMonitorSetFromOptions().

346:    Not Collective

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

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

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

363:    Example Usage:
364: .vb
365:    MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
366: .ve

368:    Then, your monitor can be chosen with the procedural interface via
369: $      MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL)
370:    or at runtime via the option
371: $      -mfn_monitor_my_monitor

373:    Level: advanced

375: .seealso: `MFNMonitorSet()`, `MFNMonitorRegisterAll()`
376: @*/
377: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,MFNMonitorRegisterFn *monitor,MFNMonitorRegisterCreateFn *create,MFNMonitorRegisterDestroyFn *destroy)
378: {
379:   char           key[PETSC_MAX_PATH_LEN];

381:   PetscFunctionBegin;
382:   PetscCall(MFNInitializePackage());
383:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
384:   PetscCall(PetscFunctionListAdd(&MFNMonitorList,key,monitor));
385:   if (create)  PetscCall(PetscFunctionListAdd(&MFNMonitorCreateList,key,create));
386:   if (destroy) PetscCall(PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

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

394:    Collective

396:    Input Parameter:
397: .  mfn - matrix function context obtained from MFNCreate()

399:    Level: advanced

401: .seealso: `MFNDestroy()`
402: @*/
403: PetscErrorCode MFNReset(MFN mfn)
404: {
405:   PetscFunctionBegin;
407:   if (!mfn) PetscFunctionReturn(PETSC_SUCCESS);
408:   PetscTryTypeMethod(mfn,reset);
409:   PetscCall(MatDestroy(&mfn->A));
410:   PetscCall(BVDestroy(&mfn->V));
411:   PetscCall(VecDestroyVecs(mfn->nwork,&mfn->work));
412:   mfn->nwork = 0;
413:   mfn->setupcalled = 0;
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: /*@
418:    MFNDestroy - Destroys the MFN context.

420:    Collective

422:    Input Parameter:
423: .  mfn - matrix function context obtained from MFNCreate()

425:    Level: beginner

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

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

447:    Collective

449:    Input Parameters:
450: +  mfn - matrix function context obtained from MFNCreate()
451: -  bv  - the basis vectors object

453:    Note:
454:    Use MFNGetBV() to retrieve the basis vectors context (for example,
455:    to free it at the end of the computations).

457:    Level: advanced

459: .seealso: `MFNGetBV()`
460: @*/
461: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
462: {
463:   PetscFunctionBegin;
466:   PetscCheckSameComm(mfn,1,bv,2);
467:   PetscCall(PetscObjectReference((PetscObject)bv));
468:   PetscCall(BVDestroy(&mfn->V));
469:   mfn->V = bv;
470:   PetscFunctionReturn(PETSC_SUCCESS);
471: }

473: /*@
474:    MFNGetBV - Obtain the basis vectors object associated to the matrix
475:    function solver.

477:    Not Collective

479:    Input Parameters:
480: .  mfn - matrix function context obtained from MFNCreate()

482:    Output Parameter:
483: .  bv - basis vectors context

485:    Level: advanced

487: .seealso: `MFNSetBV()`
488: @*/
489: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
490: {
491:   PetscFunctionBegin;
493:   PetscAssertPointer(bv,2);
494:   if (!mfn->V) {
495:     PetscCall(BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V));
496:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0));
497:     PetscCall(PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options));
498:   }
499:   *bv = mfn->V;
500:   PetscFunctionReturn(PETSC_SUCCESS);
501: }

503: /*@
504:    MFNSetFN - Specifies the function to be computed.

506:    Collective

508:    Input Parameters:
509: +  mfn - matrix function context obtained from MFNCreate()
510: -  fn  - the math function object

512:    Note:
513:    Use MFNGetFN() to retrieve the math function context (for example,
514:    to free it at the end of the computations).

516:    Level: beginner

518: .seealso: `MFNGetFN()`
519: @*/
520: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
521: {
522:   PetscFunctionBegin;
525:   PetscCheckSameComm(mfn,1,fn,2);
526:   PetscCall(PetscObjectReference((PetscObject)fn));
527:   PetscCall(FNDestroy(&mfn->fn));
528:   mfn->fn = fn;
529:   PetscFunctionReturn(PETSC_SUCCESS);
530: }

532: /*@
533:    MFNGetFN - Obtain the math function object associated to the MFN object.

535:    Not Collective

537:    Input Parameters:
538: .  mfn - matrix function context obtained from MFNCreate()

540:    Output Parameter:
541: .  fn - math function context

543:    Level: beginner

545: .seealso: `MFNSetFN()`
546: @*/
547: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
548: {
549:   PetscFunctionBegin;
551:   PetscAssertPointer(fn,2);
552:   if (!mfn->fn) {
553:     PetscCall(FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn));
554:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0));
555:     PetscCall(PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options));
556:   }
557:   *fn = mfn->fn;
558:   PetscFunctionReturn(PETSC_SUCCESS);
559: }