Actual source code: mfnbasic.c
slepc-main 2024-11-09
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: }