Actual source code: mfnbasic.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepcmfn.h" I*/

 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

 46: .seealso: PetscViewerASCIIOpen()
 47: @*/
 48: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
 49: {
 51:   PetscBool      isascii;

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

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

 86: /*@C
 87:    MFNReasonView - Displays the reason an MFN solve converged or diverged.

 89:    Collective on MFN

 91:    Parameter:
 92: +  mfn - the matrix function context
 93: -  viewer - the viewer to display the reason

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

 98:    Level: intermediate

100: .seealso: MFNSetTolerances(), MFNGetIterationNumber()
101: @*/
102: PetscErrorCode MFNReasonView(MFN mfn,PetscViewer viewer)
103: {
105:   PetscBool      isAscii;

108:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
109:   if (isAscii) {
110:     PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel);
111:     if (mfn->reason > 0) {
112:       PetscViewerASCIIPrintf(viewer,"%s Matrix function solve converged due to %s; iterations %D\n",((PetscObject)mfn)->prefix?((PetscObject)mfn)->prefix:"",MFNConvergedReasons[mfn->reason],mfn->its);
113:     } else {
114:       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);
115:     }
116:     PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel);
117:   }
118:   return(0);
119: }

121: /*@
122:    MFNReasonViewFromOptions - Processes command line options to determine if/how
123:    the MFN converged reason is to be viewed.

125:    Collective on MFN

127:    Input Parameters:
128: .  mfn - the matrix function context

130:    Level: developer
131: @*/
132: PetscErrorCode MFNReasonViewFromOptions(MFN mfn)
133: {
134:   PetscErrorCode    ierr;
135:   PetscViewer       viewer;
136:   PetscBool         flg;
137:   static PetscBool  incall = PETSC_FALSE;
138:   PetscViewerFormat format;

141:   if (incall) return(0);
142:   incall = PETSC_TRUE;
143:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg);
144:   if (flg) {
145:     PetscViewerPushFormat(viewer,format);
146:     MFNReasonView(mfn,viewer);
147:     PetscViewerPopFormat(viewer);
148:     PetscViewerDestroy(&viewer);
149:   }
150:   incall = PETSC_FALSE;
151:   return(0);
152: }

154: /*@
155:    MFNCreate - Creates the default MFN context.

157:    Collective on MPI_Comm

159:    Input Parameter:
160: .  comm - MPI communicator

162:    Output Parameter:
163: .  mfn - location to put the MFN context

165:    Note:
166:    The default MFN type is MFNKRYLOV

168:    Level: beginner

170: .seealso: MFNSetUp(), MFNSolve(), MFNDestroy(), MFN
171: @*/
172: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
173: {
175:   MFN            mfn;

179:   *outmfn = 0;
180:   MFNInitializePackage();
181:   SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView);

183:   mfn->A               = NULL;
184:   mfn->fn              = NULL;
185:   mfn->max_it          = 0;
186:   mfn->ncv             = 0;
187:   mfn->tol             = PETSC_DEFAULT;
188:   mfn->errorifnotconverged = PETSC_FALSE;

190:   mfn->numbermonitors  = 0;

192:   mfn->V               = NULL;
193:   mfn->nwork           = 0;
194:   mfn->work            = NULL;
195:   mfn->data            = NULL;

197:   mfn->its             = 0;
198:   mfn->nv              = 0;
199:   mfn->errest          = 0;
200:   mfn->setupcalled     = 0;
201:   mfn->reason          = MFN_CONVERGED_ITERATING;

203:   *outmfn = mfn;
204:   return(0);
205: }

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

210:    Logically Collective on MFN

212:    Input Parameters:
213: +  mfn  - the matrix function context
214: -  type - a known method

216:    Options Database Key:
217: .  -mfn_type <method> - Sets the method; use -help for a list
218:     of available methods

220:    Notes:
221:    See "slepc/include/slepcmfn.h" for available methods. The default
222:    is MFNKRYLOV

224:    Normally, it is best to use the MFNSetFromOptions() command and
225:    then set the MFN type from the options database rather than by using
226:    this routine.  Using the options database provides the user with
227:    maximum flexibility in evaluating the different available methods.
228:    The MFNSetType() routine is provided for those situations where it
229:    is necessary to set the iterative solver independently of the command
230:    line or options database.

232:    Level: intermediate

234: .seealso: MFNType
235: @*/
236: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
237: {
238:   PetscErrorCode ierr,(*r)(MFN);
239:   PetscBool      match;


245:   PetscObjectTypeCompare((PetscObject)mfn,type,&match);
246:   if (match) return(0);

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

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

254:   mfn->setupcalled = 0;
255:   PetscObjectChangeTypeName((PetscObject)mfn,type);
256:   (*r)(mfn);
257:   return(0);
258: }

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

263:    Not Collective

265:    Input Parameter:
266: .  mfn - the matrix function context

268:    Output Parameter:
269: .  name - name of MFN method

271:    Level: intermediate

273: .seealso: MFNSetType()
274: @*/
275: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
276: {
280:   *type = ((PetscObject)mfn)->type_name;
281:   return(0);
282: }

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

287:    Not Collective

289:    Input Parameters:
290: +  name - name of a new user-defined solver
291: -  function - routine to create the solver context

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

296:    Sample usage:
297: .vb
298:     MFNRegister("my_solver",MySolverCreate);
299: .ve

301:    Then, your solver can be chosen with the procedural interface via
302: $     MFNSetType(mfn,"my_solver")
303:    or at runtime via the option
304: $     -mfn_type my_solver

306:    Level: advanced

308: .seealso: MFNRegisterAll()
309: @*/
310: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
311: {

315:   MFNInitializePackage();
316:   PetscFunctionListAdd(&MFNList,name,function);
317:   return(0);
318: }

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

324:    Collective on MFN

326:    Input Parameter:
327: .  mfn - matrix function context obtained from MFNCreate()

329:    Level: advanced

331: .seealso: MFNDestroy()
332: @*/
333: PetscErrorCode MFNReset(MFN mfn)
334: {

339:   if (!mfn) return(0);
340:   if (mfn->ops->reset) { (mfn->ops->reset)(mfn); }
341:   MatDestroy(&mfn->A);
342:   BVDestroy(&mfn->V);
343:   VecDestroyVecs(mfn->nwork,&mfn->work);
344:   mfn->nwork = 0;
345:   mfn->setupcalled = 0;
346:   return(0);
347: }

349: /*@
350:    MFNDestroy - Destroys the MFN context.

352:    Collective on MFN

354:    Input Parameter:
355: .  mfn - matrix function context obtained from MFNCreate()

357:    Level: beginner

359: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
360: @*/
361: PetscErrorCode MFNDestroy(MFN *mfn)
362: {

366:   if (!*mfn) return(0);
368:   if (--((PetscObject)(*mfn))->refct > 0) { *mfn = 0; return(0); }
369:   MFNReset(*mfn);
370:   if ((*mfn)->ops->destroy) { (*(*mfn)->ops->destroy)(*mfn); }
371:   FNDestroy(&(*mfn)->fn);
372:   MFNMonitorCancel(*mfn);
373:   PetscHeaderDestroy(mfn);
374:   return(0);
375: }

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

380:    Collective on MFN

382:    Input Parameters:
383: +  mfn - matrix function context obtained from MFNCreate()
384: -  bv  - the basis vectors object

386:    Note:
387:    Use MFNGetBV() to retrieve the basis vectors context (for example,
388:    to free it at the end of the computations).

390:    Level: advanced

392: .seealso: MFNGetBV()
393: @*/
394: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
395: {

402:   PetscObjectReference((PetscObject)bv);
403:   BVDestroy(&mfn->V);
404:   mfn->V = bv;
405:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
406:   return(0);
407: }

409: /*@
410:    MFNGetBV - Obtain the basis vectors object associated to the matrix
411:    function solver.

413:    Not Collective

415:    Input Parameters:
416: .  mfn - matrix function context obtained from MFNCreate()

418:    Output Parameter:
419: .  bv - basis vectors context

421:    Level: advanced

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

432:   if (!mfn->V) {
433:     BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V);
434:     PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0);
435:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->V);
436:     PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options);
437:   }
438:   *bv = mfn->V;
439:   return(0);
440: }

442: /*@
443:    MFNSetFN - Specifies the function to be computed.

445:    Collective on MFN

447:    Input Parameters:
448: +  mfn - matrix function context obtained from MFNCreate()
449: -  fn  - the math function object

451:    Note:
452:    Use MFNGetFN() to retrieve the math function context (for example,
453:    to free it at the end of the computations).

455:    Level: beginner

457: .seealso: MFNGetFN()
458: @*/
459: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
460: {

467:   PetscObjectReference((PetscObject)fn);
468:   FNDestroy(&mfn->fn);
469:   mfn->fn = fn;
470:   PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
471:   return(0);
472: }

474: /*@
475:    MFNGetFN - Obtain the math function object associated to the MFN object.

477:    Not Collective

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

482:    Output Parameter:
483: .  fn - math function context

485:    Level: beginner

487: .seealso: MFNSetFN()
488: @*/
489: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
490: {

496:   if (!mfn->fn) {
497:     FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn);
498:     PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0);
499:     PetscLogObjectParent((PetscObject)mfn,(PetscObject)mfn->fn);
500:     PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options);
501:   }
502:   *fn = mfn->fn;
503:   return(0);
504: }