Actual source code: mfnbasic.c

slepc-3.13.2 2020-05-12
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:    MFNReasonView - 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:    Level: intermediate

121: .seealso: MFNSetTolerances(), MFNGetIterationNumber()
122: @*/
123: PetscErrorCode MFNReasonView(MFN mfn,PetscViewer viewer)
124: {
126:   PetscBool      isAscii;

129:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
130:   if (isAscii) {
131:     PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel);
132:     if (mfn->reason > 0) {
133:       PetscViewerASCIIPrintf(viewer,"%s Matrix function solve converged due to %s; iterations %D\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
134:     } else {
135:       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);
136:     }
137:     PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel);
138:   }
139:   return(0);
140: }

142: /*@
143:    MFNReasonViewFromOptions - Processes command line options to determine if/how
144:    the MFN converged reason is to be viewed.

146:    Collective on mfn

148:    Input Parameter:
149: .  mfn - the matrix function context

151:    Level: developer
152: @*/
153: PetscErrorCode MFNReasonViewFromOptions(MFN mfn)
154: {
155:   PetscErrorCode    ierr;
156:   PetscViewer       viewer;
157:   PetscBool         flg;
158:   static PetscBool  incall = PETSC_FALSE;
159:   PetscViewerFormat format;

162:   if (incall) return(0);
163:   incall = PETSC_TRUE;
164:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg);
165:   if (flg) {
166:     PetscViewerPushFormat(viewer,format);
167:     MFNReasonView(mfn,viewer);
168:     PetscViewerPopFormat(viewer);
169:     PetscViewerDestroy(&viewer);
170:   }
171:   incall = PETSC_FALSE;
172:   return(0);
173: }

175: /*@
176:    MFNCreate - Creates the default MFN context.

178:    Collective

180:    Input Parameter:
181: .  comm - MPI communicator

183:    Output Parameter:
184: .  mfn - location to put the MFN context

186:    Note:
187:    The default MFN type is MFNKRYLOV

189:    Level: beginner

191: .seealso: MFNSetUp(), MFNSolve(), MFNDestroy(), MFN
192: @*/
193: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
194: {
196:   MFN            mfn;

200:   *outmfn = 0;
201:   MFNInitializePackage();
202:   SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView);

204:   mfn->A               = NULL;
205:   mfn->fn              = NULL;
206:   mfn->max_it          = 0;
207:   mfn->ncv             = 0;
208:   mfn->tol             = PETSC_DEFAULT;
209:   mfn->errorifnotconverged = PETSC_FALSE;

211:   mfn->numbermonitors  = 0;

213:   mfn->V               = NULL;
214:   mfn->nwork           = 0;
215:   mfn->work            = NULL;
216:   mfn->data            = NULL;

218:   mfn->its             = 0;
219:   mfn->nv              = 0;
220:   mfn->errest          = 0;
221:   mfn->setupcalled     = 0;
222:   mfn->reason          = MFN_CONVERGED_ITERATING;

224:   *outmfn = mfn;
225:   return(0);
226: }

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

231:    Logically Collective on mfn

233:    Input Parameters:
234: +  mfn  - the matrix function context
235: -  type - a known method

237:    Options Database Key:
238: .  -mfn_type <method> - Sets the method; use -help for a list
239:     of available methods

241:    Notes:
242:    See "slepc/include/slepcmfn.h" for available methods. The default
243:    is MFNKRYLOV

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

253:    Level: intermediate

255: .seealso: MFNType
256: @*/
257: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
258: {
259:   PetscErrorCode ierr,(*r)(MFN);
260:   PetscBool      match;


266:   PetscObjectTypeCompare((PetscObject)mfn,type,&match);
267:   if (match) return(0);

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

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

275:   mfn->setupcalled = 0;
276:   PetscObjectChangeTypeName((PetscObject)mfn,type);
277:   (*r)(mfn);
278:   return(0);
279: }

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

284:    Not Collective

286:    Input Parameter:
287: .  mfn - the matrix function context

289:    Output Parameter:
290: .  name - name of MFN method

292:    Level: intermediate

294: .seealso: MFNSetType()
295: @*/
296: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
297: {
301:   *type = ((PetscObject)mfn)->type_name;
302:   return(0);
303: }

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

308:    Not Collective

310:    Input Parameters:
311: +  name - name of a new user-defined solver
312: -  function - routine to create the solver context

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

317:    Sample usage:
318: .vb
319:     MFNRegister("my_solver",MySolverCreate);
320: .ve

322:    Then, your solver can be chosen with the procedural interface via
323: $     MFNSetType(mfn,"my_solver")
324:    or at runtime via the option
325: $     -mfn_type my_solver

327:    Level: advanced

329: .seealso: MFNRegisterAll()
330: @*/
331: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
332: {

336:   MFNInitializePackage();
337:   PetscFunctionListAdd(&MFNList,name,function);
338:   return(0);
339: }

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

345:    Collective on mfn

347:    Input Parameter:
348: .  mfn - matrix function context obtained from MFNCreate()

350:    Level: advanced

352: .seealso: MFNDestroy()
353: @*/
354: PetscErrorCode MFNReset(MFN mfn)
355: {

360:   if (!mfn) return(0);
361:   if (mfn->ops->reset) { (mfn->ops->reset)(mfn); }
362:   MatDestroy(&mfn->A);
363:   BVDestroy(&mfn->V);
364:   VecDestroyVecs(mfn->nwork,&mfn->work);
365:   mfn->nwork = 0;
366:   mfn->setupcalled = 0;
367:   return(0);
368: }

370: /*@
371:    MFNDestroy - Destroys the MFN context.

373:    Collective on mfn

375:    Input Parameter:
376: .  mfn - matrix function context obtained from MFNCreate()

378:    Level: beginner

380: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
381: @*/
382: PetscErrorCode MFNDestroy(MFN *mfn)
383: {

387:   if (!*mfn) return(0);
389:   if (--((PetscObject)(*mfn))->refct > 0) { *mfn = 0; return(0); }
390:   MFNReset(*mfn);
391:   if ((*mfn)->ops->destroy) { (*(*mfn)->ops->destroy)(*mfn); }
392:   FNDestroy(&(*mfn)->fn);
393:   MatDestroy(&(*mfn)->AT);
394:   MFNMonitorCancel(*mfn);
395:   PetscHeaderDestroy(mfn);
396:   return(0);
397: }

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

402:    Collective on mfn

404:    Input Parameters:
405: +  mfn - matrix function context obtained from MFNCreate()
406: -  bv  - the basis vectors object

408:    Note:
409:    Use MFNGetBV() to retrieve the basis vectors context (for example,
410:    to free it at the end of the computations).

412:    Level: advanced

414: .seealso: MFNGetBV()
415: @*/
416: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
417: {

424:   PetscObjectReference((PetscObject)bv);
425:   BVDestroy(&mfn->V);
426:   mfn->V = bv;
427:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
428:   return(0);
429: }

431: /*@
432:    MFNGetBV - Obtain the basis vectors object associated to the matrix
433:    function solver.

435:    Not Collective

437:    Input Parameters:
438: .  mfn - matrix function context obtained from MFNCreate()

440:    Output Parameter:
441: .  bv - basis vectors context

443:    Level: advanced

445: .seealso: MFNSetBV()
446: @*/
447: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
448: {

454:   if (!mfn->V) {
455:     BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V);
456:     PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0);
457:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
458:     PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options);
459:   }
460:   *bv = mfn->V;
461:   return(0);
462: }

464: /*@
465:    MFNSetFN - Specifies the function to be computed.

467:    Collective on mfn

469:    Input Parameters:
470: +  mfn - matrix function context obtained from MFNCreate()
471: -  fn  - the math function object

473:    Note:
474:    Use MFNGetFN() to retrieve the math function context (for example,
475:    to free it at the end of the computations).

477:    Level: beginner

479: .seealso: MFNGetFN()
480: @*/
481: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
482: {

489:   PetscObjectReference((PetscObject)fn);
490:   FNDestroy(&mfn->fn);
491:   mfn->fn = fn;
492:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
493:   return(0);
494: }

496: /*@
497:    MFNGetFN - Obtain the math function object associated to the MFN object.

499:    Not Collective

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

504:    Output Parameter:
505: .  fn - math function context

507:    Level: beginner

509: .seealso: MFNSetFN()
510: @*/
511: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
512: {

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