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: Notes:
43: The available visualization contexts include
44: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
45: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
46: first process opens the file; all other processes send their data to the
47: first one to print
49: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
50: to output to a specified file.
52: Level: beginner
54: .seealso: [](ch:mfn), `MFNCreate()`, `MFNViewFromOptions()`
55: @*/
56: PetscErrorCode MFNView(MFN mfn,PetscViewer viewer)
57: {
58: PetscBool isascii;
60: PetscFunctionBegin;
62: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mfn),&viewer));
64: PetscCheckSameComm(mfn,1,viewer,2);
66: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
67: if (isascii) {
68: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)mfn,viewer));
69: PetscCall(PetscViewerASCIIPushTab(viewer));
70: PetscTryTypeMethod(mfn,view,viewer);
71: PetscCall(PetscViewerASCIIPopTab(viewer));
72: PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",mfn->ncv));
73: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",mfn->max_it));
74: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)mfn->tol));
75: } else PetscTryTypeMethod(mfn,view,viewer);
76: PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
77: if (!mfn->V) PetscCall(MFNGetFN(mfn,&mfn->fn));
78: PetscCall(FNView(mfn->fn,viewer));
79: if (!mfn->V) PetscCall(MFNGetBV(mfn,&mfn->V));
80: PetscCall(BVView(mfn->V,viewer));
81: PetscCall(PetscViewerPopFormat(viewer));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: MFNViewFromOptions - View (print) an `MFN` object based on values in the options database.
88: Collective
90: Input Parameters:
91: + mfn - the matrix function solver context
92: . obj - optional object that provides the options prefix used to query the options database
93: - name - command line option
95: Level: intermediate
97: .seealso: [](ch:mfn), `MFNView()`, `MFNCreate()`, `PetscObjectViewFromOptions()`
98: @*/
99: PetscErrorCode MFNViewFromOptions(MFN mfn,PetscObject obj,const char name[])
100: {
101: PetscFunctionBegin;
103: PetscCall(PetscObjectViewFromOptions((PetscObject)mfn,obj,name));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
106: /*@
107: MFNConvergedReasonView - Displays the reason an `MFN` solve converged or diverged.
109: Collective
111: Input Parameters:
112: + mfn - the matrix function solver context
113: - viewer - the viewer to display the reason
115: Options Database Key:
116: . -mfn_converged_reason - print reason for convergence/divergence, and number of iterations
118: Notes:
119: Use `MFNConvergedReasonViewFromOptions()` to display the reason based on values
120: in the options database.
122: To change the format of the output call `PetscViewerPushFormat()` before this
123: call. Use `PETSC_VIEWER_DEFAULT` for the default, or `PETSC_VIEWER_FAILED` to only
124: display a reason if it fails. The latter can be set in the command line with
125: `-mfn_converged_reason ::failed`.
127: Level: intermediate
129: .seealso: [](ch:mfn), `MFNSetTolerances()`, `MFNGetIterationNumber()`, `MFNConvergedReasonViewFromOptions()`
130: @*/
131: PetscErrorCode MFNConvergedReasonView(MFN mfn,PetscViewer viewer)
132: {
133: PetscBool isAscii;
134: PetscViewerFormat format;
136: PetscFunctionBegin;
137: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)mfn));
138: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
139: if (isAscii) {
140: PetscCall(PetscViewerGetFormat(viewer,&format));
141: PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)mfn)->tablevel));
142: 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));
143: 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));
144: PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)mfn)->tablevel));
145: }
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: MFNConvergedReasonViewFromOptions - Processes command line options to determine if/how
151: the `MFN` converged reason is to be viewed.
153: Collective
155: Input Parameter:
156: . mfn - the matrix function solver context
158: Level: intermediate
160: .seealso: [](ch:mfn), `MFNConvergedReasonView()`
161: @*/
162: PetscErrorCode MFNConvergedReasonViewFromOptions(MFN mfn)
163: {
164: PetscViewer viewer;
165: PetscBool flg;
166: static PetscBool incall = PETSC_FALSE;
167: PetscViewerFormat format;
169: PetscFunctionBegin;
170: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
171: incall = PETSC_TRUE;
172: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)mfn),((PetscObject)mfn)->options,((PetscObject)mfn)->prefix,"-mfn_converged_reason",&viewer,&format,&flg));
173: if (flg) {
174: PetscCall(PetscViewerPushFormat(viewer,format));
175: PetscCall(MFNConvergedReasonView(mfn,viewer));
176: PetscCall(PetscViewerPopFormat(viewer));
177: PetscCall(PetscViewerDestroy(&viewer));
178: }
179: incall = PETSC_FALSE;
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: MFNCreate - Creates the `MFN` context.
186: Collective
188: Input Parameter:
189: . comm - MPI communicator
191: Output Parameter:
192: . outmfn - location to put the `MFN` context
194: Note:
195: The default `MFN` type is `MFNKRYLOV`.
197: Level: beginner
199: .seealso: [](ch:mfn), `MFNSetUp()`, `MFNSolve()`, `MFNDestroy()`, `MFN`
200: @*/
201: PetscErrorCode MFNCreate(MPI_Comm comm,MFN *outmfn)
202: {
203: MFN mfn;
205: PetscFunctionBegin;
206: PetscAssertPointer(outmfn,2);
207: PetscCall(MFNInitializePackage());
208: PetscCall(SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView));
210: mfn->A = NULL;
211: mfn->fn = NULL;
212: mfn->max_it = PETSC_DETERMINE;
213: mfn->ncv = PETSC_DETERMINE;
214: mfn->tol = PETSC_DETERMINE;
215: mfn->errorifnotconverged = PETSC_FALSE;
217: mfn->numbermonitors = 0;
219: mfn->V = NULL;
220: mfn->nwork = 0;
221: mfn->work = NULL;
222: mfn->data = NULL;
224: mfn->its = 0;
225: mfn->nv = 0;
226: mfn->errest = 0;
227: mfn->setupcalled = 0;
228: mfn->reason = MFN_CONVERGED_ITERATING;
230: *outmfn = mfn;
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: MFNSetType - Selects the particular solver to be used in the `MFN` object.
237: Logically Collective
239: Input Parameters:
240: + mfn - the matrix function solver context
241: - type - a known method
243: Options Database Key:
244: . -mfn_type type - sets the method; use `-help` for a list of available methods
246: Notes:
247: See `MFNType` for available methods. The default 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: [](ch:mfn), `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 solver context
293: Output Parameter:
294: . type - name of `MFN` method
296: Note:
297: `type` should not be retained for later use as it will be an invalid pointer
298: if the `MFNType` of `mfn` is changed.
300: Level: intermediate
302: .seealso: [](ch:mfn), `MFNSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
303: @*/
304: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
305: {
306: PetscFunctionBegin;
308: PetscAssertPointer(type,2);
309: *type = ((PetscObject)mfn)->type_name;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@C
314: MFNRegister - Adds a method to the matrix function solver package.
316: Not Collective
318: Input Parameters:
319: + name - name of a new user-defined solver
320: - function - routine to create the solver context
322: Note:
323: `MFNRegister()` may be called multiple times to add several user-defined solvers.
325: Example Usage:
326: .vb
327: MFNRegister("my_solver",MySolverCreate);
328: .ve
330: Then, your solver can be chosen with the procedural interface via
331: .vb
332: MFNSetType(mfn,"my_solver")
333: .ve
334: or at runtime via the option `-mfn_type my_solver`.
336: Level: advanced
338: .seealso: [](ch:mfn), `MFNRegisterAll()`
339: @*/
340: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
341: {
342: PetscFunctionBegin;
343: PetscCall(MFNInitializePackage());
344: PetscCall(PetscFunctionListAdd(&MFNList,name,function));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@C
349: MFNMonitorRegister - Registers an `MFN` monitor routine that may be accessed with
350: `MFNMonitorSetFromOptions()`.
352: Not Collective
354: Input Parameters:
355: + name - name of a new monitor routine
356: . vtype - a `PetscViewerType` for the output
357: . format - a `PetscViewerFormat` for the output
358: . monitor - monitor routine, see `MFNMonitorRegisterFn`
359: . create - creation routine, or `NULL`
360: - destroy - destruction routine, or `NULL`
362: Notes:
363: `MFNMonitorRegister()` may be called multiple times to add several user-defined monitors.
365: The calling sequence for the given function matches the calling sequence of `MFNMonitorFn`
366: functions passed to `MFNMonitorSet()` with the additional requirement that its final argument
367: be a `PetscViewerAndFormat`.
369: Example Usage:
370: .vb
371: MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
372: .ve
374: Then, your monitor can be chosen with the procedural interface via
375: .vb
376: MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL);
377: .ve
378: or at runtime via the option `-mfn_monitor_my_monitor`.
380: Level: advanced
382: .seealso: [](ch:mfn), `MFNMonitorSet()`, `MFNMonitorRegisterAll()`, `MFNMonitorSetFromOptions()`
383: @*/
384: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,MFNMonitorRegisterFn *monitor,MFNMonitorRegisterCreateFn *create,MFNMonitorRegisterDestroyFn *destroy)
385: {
386: char key[PETSC_MAX_PATH_LEN];
388: PetscFunctionBegin;
389: PetscCall(MFNInitializePackage());
390: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
391: PetscCall(PetscFunctionListAdd(&MFNMonitorList,key,monitor));
392: if (create) PetscCall(PetscFunctionListAdd(&MFNMonitorCreateList,key,create));
393: if (destroy) PetscCall(PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /*@
398: MFNReset - Resets the `MFN` context to the initial state (prior to setup)
399: and destroys any allocated `Vec`s and `Mat`s.
401: Collective
403: Input Parameter:
404: . mfn - the matrix function solver context
406: Level: advanced
408: .seealso: [](ch:mfn), `MFNDestroy()`
409: @*/
410: PetscErrorCode MFNReset(MFN mfn)
411: {
412: PetscFunctionBegin;
414: if (!mfn) PetscFunctionReturn(PETSC_SUCCESS);
415: PetscTryTypeMethod(mfn,reset);
416: PetscCall(MatDestroy(&mfn->A));
417: PetscCall(BVDestroy(&mfn->V));
418: PetscCall(VecDestroyVecs(mfn->nwork,&mfn->work));
419: mfn->nwork = 0;
420: mfn->setupcalled = 0;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: MFNDestroy - Destroys the `MFN` context.
427: Collective
429: Input Parameter:
430: . mfn - the matrix function solver context
432: Level: beginner
434: .seealso: [](ch:mfn), `MFNCreate()`, `MFNSetUp()`, `MFNSolve()`
435: @*/
436: PetscErrorCode MFNDestroy(MFN *mfn)
437: {
438: PetscFunctionBegin;
439: if (!*mfn) PetscFunctionReturn(PETSC_SUCCESS);
441: if (--((PetscObject)*mfn)->refct > 0) { *mfn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
442: PetscCall(MFNReset(*mfn));
443: PetscTryTypeMethod(*mfn,destroy);
444: PetscCall(FNDestroy(&(*mfn)->fn));
445: PetscCall(MatDestroy(&(*mfn)->AT));
446: PetscCall(MFNMonitorCancel(*mfn));
447: PetscCall(PetscHeaderDestroy(mfn));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@
452: MFNSetBV - Associates a basis vectors object to the matrix function solver.
454: Collective
456: Input Parameters:
457: + mfn - the matrix function solver context
458: - bv - the basis vectors object
460: Note:
461: Use `MFNGetBV()` to retrieve the basis vectors context (for example,
462: to free it at the end of the computations).
464: Level: advanced
466: .seealso: [](ch:mfn), `MFNGetBV()`
467: @*/
468: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
469: {
470: PetscFunctionBegin;
473: PetscCheckSameComm(mfn,1,bv,2);
474: PetscCall(PetscObjectReference((PetscObject)bv));
475: PetscCall(BVDestroy(&mfn->V));
476: mfn->V = bv;
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*@
481: MFNGetBV - Obtain the basis vectors object associated to the matrix
482: function solver.
484: Not Collective
486: Input Parameter:
487: . mfn - the matrix function solver context
489: Output Parameter:
490: . bv - basis vectors context
492: Level: advanced
494: .seealso: [](ch:mfn), `MFNSetBV()`
495: @*/
496: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
497: {
498: PetscFunctionBegin;
500: PetscAssertPointer(bv,2);
501: if (!mfn->V) {
502: PetscCall(BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V));
503: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0));
504: PetscCall(PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options));
505: }
506: *bv = mfn->V;
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: /*@
511: MFNSetFN - Specifies the function to be computed.
513: Collective
515: Input Parameters:
516: + mfn - the matrix function solver context
517: - fn - the math function object
519: Notes:
520: At a later time, use `MFNGetFN()` to retrieve the math function context
521: (for example, to free it at the end of the computations).
523: This function is not called in normal usage. Instead, it is easier to
524: extract the internal `FN` object with `MFNGetFN()` and modify it.
526: Level: beginner
528: .seealso: [](ch:mfn), `MFNGetFN()`
529: @*/
530: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
531: {
532: PetscFunctionBegin;
535: PetscCheckSameComm(mfn,1,fn,2);
536: PetscCall(PetscObjectReference((PetscObject)fn));
537: PetscCall(FNDestroy(&mfn->fn));
538: mfn->fn = fn;
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*@
543: MFNGetFN - Obtain the math function object associated to the `MFN` object.
545: Not Collective
547: Input Parameter:
548: . mfn - the matrix function solver context
550: Output Parameter:
551: . fn - math function context
553: Note:
554: This is the usual way to specify the function that needs to be applied
555: to a given vector in `MFNSolve()`.
557: Level: beginner
559: .seealso: [](ch:mfn), `MFNSetFN()`, `MFNSolve()`
560: @*/
561: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
562: {
563: PetscFunctionBegin;
565: PetscAssertPointer(fn,2);
566: if (!mfn->fn) {
567: PetscCall(FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn));
568: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0));
569: PetscCall(PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options));
570: }
571: *fn = mfn->fn;
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }