Actual source code: mfnbasic.c

slepc-3.14.0 2020-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: PetscFunctionList MFNList = 0;
 17: PetscBool         MFNRegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      MFN_CLASSID = 0;
 19: PetscLogEvent     MFN_SetUp = 0,MFN_Solve = 0;

 21: /*@C
 22:    MFNView - Prints the MFN data structure.

 24:    Collective on mfn

 26:    Input Parameters:
 27: +  mfn - the matrix function solver context
 28: -  viewer - optional visualization context

 30:    Options Database Key:
 31: .  -mfn_view -  Calls MFNView() at end of MFNSolve()

 33:    Note:
 34:    The available visualization contexts include
 35: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 36: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 37:          output where only the first processor opens
 38:          the file.  All other processors send their
 39:          data to the first processor to print.

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

 44:    Level: beginner
 45: @*/
 46: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
 47: {
 49:   PetscBool      isascii;

 53:   if (!viewer) {
 54:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer);
 55:   }

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

 84: /*@C
 85:    MFNViewFromOptions - View from options

 87:    Collective on MFN

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

 94:    Level: intermediate

 96: .seealso: MFNView(), MFNCreate()
 97: @*/
 98: PetscErrorCode MFNViewFromOptions(MFN mfn,PetscObject obj,const char name[])
 99: {

104:   PetscObjectViewFromOptions((PetscObject)mfn,obj,name);
105:   return(0);
106: }
107: /*@C
108:    MFNConvergedReasonView - Displays the reason an MFN solve converged or diverged.

110:    Collective on mfn

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:   PetscErrorCode    ierr;
132:   PetscBool         isAscii;
133:   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) {
142:       PetscViewerASCIIPrintf(viewer,"%s Matrix function solve converged due to %s; iterations %D\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
143:     } else if (mfn->reason <= 0) {
144:       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);
145:     }
146:     PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel);
147:   }
148:   return(0);
149: }

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

155:    Collective on mfn

157:    Input Parameter:
158: .  mfn - the matrix function context

160:    Level: developer

162: .seealso: MFNConvergedReasonView()
163: @*/
164: PetscErrorCode MFNConvergedReasonViewFromOptions(MFN mfn)
165: {
166:   PetscErrorCode    ierr;
167:   PetscViewer       viewer;
168:   PetscBool         flg;
169:   static PetscBool  incall = PETSC_FALSE;
170:   PetscViewerFormat format;

173:   if (incall) return(0);
174:   incall = PETSC_TRUE;
175:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg);
176:   if (flg) {
177:     PetscViewerPushFormat(viewer,format);
178:     MFNConvergedReasonView(mfn,viewer);
179:     PetscViewerPopFormat(viewer);
180:     PetscViewerDestroy(&viewer);
181:   }
182:   incall = PETSC_FALSE;
183:   return(0);
184: }

186: /*@
187:    MFNCreate - Creates the default MFN context.

189:    Collective

191:    Input Parameter:
192: .  comm - MPI communicator

194:    Output Parameter:
195: .  mfn - location to put the MFN context

197:    Note:
198:    The default MFN type is MFNKRYLOV

200:    Level: beginner

202: .seealso: MFNSetUp(), MFNSolve(), MFNDestroy(), MFN
203: @*/
204: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
205: {
207:   MFN            mfn;

211:   *outmfn = 0;
212:   MFNInitializePackage();
213:   SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView);

215:   mfn->A               = NULL;
216:   mfn->fn              = NULL;
217:   mfn->max_it          = PETSC_DEFAULT;
218:   mfn->ncv             = PETSC_DEFAULT;
219:   mfn->tol             = PETSC_DEFAULT;
220:   mfn->errorifnotconverged = PETSC_FALSE;

222:   mfn->numbermonitors  = 0;

224:   mfn->V               = NULL;
225:   mfn->nwork           = 0;
226:   mfn->work            = NULL;
227:   mfn->data            = NULL;

229:   mfn->its             = 0;
230:   mfn->nv              = 0;
231:   mfn->errest          = 0;
232:   mfn->setupcalled     = 0;
233:   mfn->reason          = MFN_CONVERGED_ITERATING;

235:   *outmfn = mfn;
236:   return(0);
237: }

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

242:    Logically Collective on mfn

244:    Input Parameters:
245: +  mfn  - the matrix function context
246: -  type - a known method

248:    Options Database Key:
249: .  -mfn_type <method> - Sets the method; use -help for a list
250:     of available methods

252:    Notes:
253:    See "slepc/include/slepcmfn.h" for available methods. The default
254:    is MFNKRYLOV

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

264:    Level: intermediate

266: .seealso: MFNType
267: @*/
268: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
269: {
270:   PetscErrorCode ierr,(*r)(MFN);
271:   PetscBool      match;


277:   PetscObjectTypeCompare((PetscObject)mfn,type,&match);
278:   if (match) return(0);

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

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

286:   mfn->setupcalled = 0;
287:   PetscObjectChangeTypeName((PetscObject)mfn,type);
288:   (*r)(mfn);
289:   return(0);
290: }

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

295:    Not Collective

297:    Input Parameter:
298: .  mfn - the matrix function context

300:    Output Parameter:
301: .  name - name of MFN method

303:    Level: intermediate

305: .seealso: MFNSetType()
306: @*/
307: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
308: {
312:   *type = ((PetscObject)mfn)->type_name;
313:   return(0);
314: }

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

319:    Not Collective

321:    Input Parameters:
322: +  name - name of a new user-defined solver
323: -  function - routine to create the solver context

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

328:    Sample usage:
329: .vb
330:     MFNRegister("my_solver",MySolverCreate);
331: .ve

333:    Then, your solver can be chosen with the procedural interface via
334: $     MFNSetType(mfn,"my_solver")
335:    or at runtime via the option
336: $     -mfn_type my_solver

338:    Level: advanced

340: .seealso: MFNRegisterAll()
341: @*/
342: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
343: {

347:   MFNInitializePackage();
348:   PetscFunctionListAdd(&MFNList,name,function);
349:   return(0);
350: }

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

356:    Collective on mfn

358:    Input Parameter:
359: .  mfn - matrix function context obtained from MFNCreate()

361:    Level: advanced

363: .seealso: MFNDestroy()
364: @*/
365: PetscErrorCode MFNReset(MFN mfn)
366: {

371:   if (!mfn) return(0);
372:   if (mfn->ops->reset) { (mfn->ops->reset)(mfn); }
373:   MatDestroy(&mfn->A);
374:   BVDestroy(&mfn->V);
375:   VecDestroyVecs(mfn->nwork,&mfn->work);
376:   mfn->nwork = 0;
377:   mfn->setupcalled = 0;
378:   return(0);
379: }

381: /*@
382:    MFNDestroy - Destroys the MFN context.

384:    Collective on mfn

386:    Input Parameter:
387: .  mfn - matrix function context obtained from MFNCreate()

389:    Level: beginner

391: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
392: @*/
393: PetscErrorCode MFNDestroy(MFN *mfn)
394: {

398:   if (!*mfn) return(0);
400:   if (--((PetscObject)(*mfn))->refct > 0) { *mfn = 0; return(0); }
401:   MFNReset(*mfn);
402:   if ((*mfn)->ops->destroy) { (*(*mfn)->ops->destroy)(*mfn); }
403:   FNDestroy(&(*mfn)->fn);
404:   MatDestroy(&(*mfn)->AT);
405:   MFNMonitorCancel(*mfn);
406:   PetscHeaderDestroy(mfn);
407:   return(0);
408: }

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

413:    Collective on mfn

415:    Input Parameters:
416: +  mfn - matrix function context obtained from MFNCreate()
417: -  bv  - the basis vectors object

419:    Note:
420:    Use MFNGetBV() to retrieve the basis vectors context (for example,
421:    to free it at the end of the computations).

423:    Level: advanced

425: .seealso: MFNGetBV()
426: @*/
427: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
428: {

435:   PetscObjectReference((PetscObject)bv);
436:   BVDestroy(&mfn->V);
437:   mfn->V = bv;
438:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
439:   return(0);
440: }

442: /*@
443:    MFNGetBV - Obtain the basis vectors object associated to the matrix
444:    function solver.

446:    Not Collective

448:    Input Parameters:
449: .  mfn - matrix function context obtained from MFNCreate()

451:    Output Parameter:
452: .  bv - basis vectors context

454:    Level: advanced

456: .seealso: MFNSetBV()
457: @*/
458: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
459: {

465:   if (!mfn->V) {
466:     BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V);
467:     PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0);
468:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
469:     PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options);
470:   }
471:   *bv = mfn->V;
472:   return(0);
473: }

475: /*@
476:    MFNSetFN - Specifies the function to be computed.

478:    Collective on mfn

480:    Input Parameters:
481: +  mfn - matrix function context obtained from MFNCreate()
482: -  fn  - the math function object

484:    Note:
485:    Use MFNGetFN() to retrieve the math function context (for example,
486:    to free it at the end of the computations).

488:    Level: beginner

490: .seealso: MFNGetFN()
491: @*/
492: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
493: {

500:   PetscObjectReference((PetscObject)fn);
501:   FNDestroy(&mfn->fn);
502:   mfn->fn = fn;
503:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
504:   return(0);
505: }

507: /*@
508:    MFNGetFN - Obtain the math function object associated to the MFN object.

510:    Not Collective

512:    Input Parameters:
513: .  mfn - matrix function context obtained from MFNCreate()

515:    Output Parameter:
516: .  fn - math function context

518:    Level: beginner

520: .seealso: MFNSetFN()
521: @*/
522: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
523: {

529:   if (!mfn->fn) {
530:     FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn);
531:     PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0);
532:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
533:     PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options);
534:   }
535:   *fn = mfn->fn;
536:   return(0);
537: }