Actual source code: mfnbasic.c

slepc-main 2024-11-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 = 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 - Adds MFN monitor routine.

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
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:    Example Usage:
360: .vb
361:    MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
362: .ve

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

369:    Level: advanced

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

377:   PetscFunctionBegin;
378:   PetscCall(MFNInitializePackage());
379:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
380:   PetscCall(PetscFunctionListAdd(&MFNMonitorList,key,monitor));
381:   if (create)  PetscCall(PetscFunctionListAdd(&MFNMonitorCreateList,key,create));
382:   if (destroy) PetscCall(PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy));
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

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

390:    Collective

392:    Input Parameter:
393: .  mfn - matrix function context obtained from MFNCreate()

395:    Level: advanced

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

413: /*@
414:    MFNDestroy - Destroys the MFN context.

416:    Collective

418:    Input Parameter:
419: .  mfn - matrix function context obtained from MFNCreate()

421:    Level: beginner

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

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

443:    Collective

445:    Input Parameters:
446: +  mfn - matrix function context obtained from MFNCreate()
447: -  bv  - the basis vectors object

449:    Note:
450:    Use MFNGetBV() to retrieve the basis vectors context (for example,
451:    to free it at the end of the computations).

453:    Level: advanced

455: .seealso: MFNGetBV()
456: @*/
457: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
458: {
459:   PetscFunctionBegin;
462:   PetscCheckSameComm(mfn,1,bv,2);
463:   PetscCall(PetscObjectReference((PetscObject)bv));
464:   PetscCall(BVDestroy(&mfn->V));
465:   mfn->V = bv;
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

469: /*@
470:    MFNGetBV - Obtain the basis vectors object associated to the matrix
471:    function solver.

473:    Not Collective

475:    Input Parameters:
476: .  mfn - matrix function context obtained from MFNCreate()

478:    Output Parameter:
479: .  bv - basis vectors context

481:    Level: advanced

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

499: /*@
500:    MFNSetFN - Specifies the function to be computed.

502:    Collective

504:    Input Parameters:
505: +  mfn - matrix function context obtained from MFNCreate()
506: -  fn  - the math function object

508:    Note:
509:    Use MFNGetFN() to retrieve the math function context (for example,
510:    to free it at the end of the computations).

512:    Level: beginner

514: .seealso: MFNGetFN()
515: @*/
516: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
517: {
518:   PetscFunctionBegin;
521:   PetscCheckSameComm(mfn,1,fn,2);
522:   PetscCall(PetscObjectReference((PetscObject)fn));
523:   PetscCall(FNDestroy(&mfn->fn));
524:   mfn->fn = fn;
525:   PetscFunctionReturn(PETSC_SUCCESS);
526: }

528: /*@
529:    MFNGetFN - Obtain the math function object associated to the MFN object.

531:    Not Collective

533:    Input Parameters:
534: .  mfn - matrix function context obtained from MFNCreate()

536:    Output Parameter:
537: .  fn - math function context

539:    Level: beginner

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