Actual source code: mfnbasic.c
slepc-3.21.0 2024-03-30
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: /*@C
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: /*@C
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: /*@C
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(PetscOptionsGetViewer(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(PetscOptionsRestoreViewer(&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: *outmfn = NULL;
206: PetscCall(MFNInitializePackage());
207: PetscCall(SlepcHeaderCreate(mfn,MFN_CLASSID,"MFN","Matrix Function","MFN",comm,MFNDestroy,MFNView));
209: mfn->A = NULL;
210: mfn->fn = NULL;
211: mfn->max_it = PETSC_DEFAULT;
212: mfn->ncv = PETSC_DEFAULT;
213: mfn->tol = PETSC_DEFAULT;
214: mfn->errorifnotconverged = PETSC_FALSE;
216: mfn->numbermonitors = 0;
218: mfn->V = NULL;
219: mfn->nwork = 0;
220: mfn->work = NULL;
221: mfn->data = NULL;
223: mfn->its = 0;
224: mfn->nv = 0;
225: mfn->errest = 0;
226: mfn->setupcalled = 0;
227: mfn->reason = MFN_CONVERGED_ITERATING;
229: *outmfn = mfn;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*@C
234: MFNSetType - Selects the particular solver to be used in the MFN object.
236: Logically Collective
238: Input Parameters:
239: + mfn - the matrix function context
240: - type - a known method
242: Options Database Key:
243: . -mfn_type <method> - Sets the method; use -help for a list
244: of available methods
246: Notes:
247: See "slepc/include/slepcmfn.h" for available methods. The default
248: is MFNKRYLOV
250: Normally, it is best to use the MFNSetFromOptions() command and
251: then set the MFN type from the options database rather than by using
252: this routine. Using the options database provides the user with
253: maximum flexibility in evaluating the different available methods.
254: The MFNSetType() routine is provided for those situations where it
255: is necessary to set the iterative solver independently of the command
256: line or options database.
258: Level: intermediate
260: .seealso: MFNType
261: @*/
262: PetscErrorCode MFNSetType(MFN mfn,MFNType type)
263: {
264: PetscErrorCode (*r)(MFN);
265: PetscBool match;
267: PetscFunctionBegin;
269: PetscAssertPointer(type,2);
271: PetscCall(PetscObjectTypeCompare((PetscObject)mfn,type,&match));
272: if (match) PetscFunctionReturn(PETSC_SUCCESS);
274: PetscCall(PetscFunctionListFind(MFNList,type,&r));
275: PetscCheck(r,PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MFN type given: %s",type);
277: PetscTryTypeMethod(mfn,destroy);
278: PetscCall(PetscMemzero(mfn->ops,sizeof(struct _MFNOps)));
280: mfn->setupcalled = 0;
281: PetscCall(PetscObjectChangeTypeName((PetscObject)mfn,type));
282: PetscCall((*r)(mfn));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@C
287: MFNGetType - Gets the MFN type as a string from the MFN object.
289: Not Collective
291: Input Parameter:
292: . mfn - the matrix function context
294: Output Parameter:
295: . type - name of MFN method
297: Level: intermediate
299: .seealso: MFNSetType()
300: @*/
301: PetscErrorCode MFNGetType(MFN mfn,MFNType *type)
302: {
303: PetscFunctionBegin;
305: PetscAssertPointer(type,2);
306: *type = ((PetscObject)mfn)->type_name;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@C
311: MFNRegister - Adds a method to the matrix function solver package.
313: Not Collective
315: Input Parameters:
316: + name - name of a new user-defined solver
317: - function - routine to create the solver context
319: Notes:
320: MFNRegister() may be called multiple times to add several user-defined solvers.
322: Example Usage:
323: .vb
324: MFNRegister("my_solver",MySolverCreate);
325: .ve
327: Then, your solver can be chosen with the procedural interface via
328: $ MFNSetType(mfn,"my_solver")
329: or at runtime via the option
330: $ -mfn_type my_solver
332: Level: advanced
334: .seealso: MFNRegisterAll()
335: @*/
336: PetscErrorCode MFNRegister(const char *name,PetscErrorCode (*function)(MFN))
337: {
338: PetscFunctionBegin;
339: PetscCall(MFNInitializePackage());
340: PetscCall(PetscFunctionListAdd(&MFNList,name,function));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /*@C
345: MFNMonitorRegister - Adds MFN monitor routine.
347: Not Collective
349: Input Parameters:
350: + name - name of a new monitor routine
351: . vtype - a PetscViewerType for the output
352: . format - a PetscViewerFormat for the output
353: . monitor - monitor routine
354: . create - creation routine, or NULL
355: - destroy - destruction routine, or NULL
357: Notes:
358: MFNMonitorRegister() may be called multiple times to add several user-defined monitors.
360: Example Usage:
361: .vb
362: MFNMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
363: .ve
365: Then, your monitor can be chosen with the procedural interface via
366: $ MFNMonitorSetFromOptions(mfn,"-mfn_monitor_my_monitor","my_monitor",NULL)
367: or at runtime via the option
368: $ -mfn_monitor_my_monitor
370: Level: advanced
372: .seealso: MFNMonitorRegisterAll()
373: @*/
374: PetscErrorCode MFNMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(MFN,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
375: {
376: char key[PETSC_MAX_PATH_LEN];
378: PetscFunctionBegin;
379: PetscCall(MFNInitializePackage());
380: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
381: PetscCall(PetscFunctionListAdd(&MFNMonitorList,key,monitor));
382: if (create) PetscCall(PetscFunctionListAdd(&MFNMonitorCreateList,key,create));
383: if (destroy) PetscCall(PetscFunctionListAdd(&MFNMonitorDestroyList,key,destroy));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@
388: MFNReset - Resets the MFN context to the initial state (prior to setup)
389: and destroys any allocated Vecs and Mats.
391: Collective
393: Input Parameter:
394: . mfn - matrix function context obtained from MFNCreate()
396: Level: advanced
398: .seealso: MFNDestroy()
399: @*/
400: PetscErrorCode MFNReset(MFN mfn)
401: {
402: PetscFunctionBegin;
404: if (!mfn) PetscFunctionReturn(PETSC_SUCCESS);
405: PetscTryTypeMethod(mfn,reset);
406: PetscCall(MatDestroy(&mfn->A));
407: PetscCall(BVDestroy(&mfn->V));
408: PetscCall(VecDestroyVecs(mfn->nwork,&mfn->work));
409: mfn->nwork = 0;
410: mfn->setupcalled = 0;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*@C
415: MFNDestroy - Destroys the MFN context.
417: Collective
419: Input Parameter:
420: . mfn - matrix function context obtained from MFNCreate()
422: Level: beginner
424: .seealso: MFNCreate(), MFNSetUp(), MFNSolve()
425: @*/
426: PetscErrorCode MFNDestroy(MFN *mfn)
427: {
428: PetscFunctionBegin;
429: if (!*mfn) PetscFunctionReturn(PETSC_SUCCESS);
431: if (--((PetscObject)*mfn)->refct > 0) { *mfn = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
432: PetscCall(MFNReset(*mfn));
433: PetscTryTypeMethod(*mfn,destroy);
434: PetscCall(FNDestroy(&(*mfn)->fn));
435: PetscCall(MatDestroy(&(*mfn)->AT));
436: PetscCall(MFNMonitorCancel(*mfn));
437: PetscCall(PetscHeaderDestroy(mfn));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: /*@
442: MFNSetBV - Associates a basis vectors object to the matrix function solver.
444: Collective
446: Input Parameters:
447: + mfn - matrix function context obtained from MFNCreate()
448: - bv - the basis vectors object
450: Note:
451: Use MFNGetBV() to retrieve the basis vectors context (for example,
452: to free it at the end of the computations).
454: Level: advanced
456: .seealso: MFNGetBV()
457: @*/
458: PetscErrorCode MFNSetBV(MFN mfn,BV bv)
459: {
460: PetscFunctionBegin;
463: PetscCheckSameComm(mfn,1,bv,2);
464: PetscCall(PetscObjectReference((PetscObject)bv));
465: PetscCall(BVDestroy(&mfn->V));
466: mfn->V = bv;
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: MFNGetBV - Obtain the basis vectors object associated to the matrix
472: function solver.
474: Not Collective
476: Input Parameters:
477: . mfn - matrix function context obtained from MFNCreate()
479: Output Parameter:
480: . bv - basis vectors context
482: Level: advanced
484: .seealso: MFNSetBV()
485: @*/
486: PetscErrorCode MFNGetBV(MFN mfn,BV *bv)
487: {
488: PetscFunctionBegin;
490: PetscAssertPointer(bv,2);
491: if (!mfn->V) {
492: PetscCall(BVCreate(PetscObjectComm((PetscObject)mfn),&mfn->V));
493: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->V,(PetscObject)mfn,0));
494: PetscCall(PetscObjectSetOptions((PetscObject)mfn->V,((PetscObject)mfn)->options));
495: }
496: *bv = mfn->V;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: MFNSetFN - Specifies the function to be computed.
503: Collective
505: Input Parameters:
506: + mfn - matrix function context obtained from MFNCreate()
507: - fn - the math function object
509: Note:
510: Use MFNGetFN() to retrieve the math function context (for example,
511: to free it at the end of the computations).
513: Level: beginner
515: .seealso: MFNGetFN()
516: @*/
517: PetscErrorCode MFNSetFN(MFN mfn,FN fn)
518: {
519: PetscFunctionBegin;
522: PetscCheckSameComm(mfn,1,fn,2);
523: PetscCall(PetscObjectReference((PetscObject)fn));
524: PetscCall(FNDestroy(&mfn->fn));
525: mfn->fn = fn;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: MFNGetFN - Obtain the math function object associated to the MFN object.
532: Not Collective
534: Input Parameters:
535: . mfn - matrix function context obtained from MFNCreate()
537: Output Parameter:
538: . fn - math function context
540: Level: beginner
542: .seealso: MFNSetFN()
543: @*/
544: PetscErrorCode MFNGetFN(MFN mfn,FN *fn)
545: {
546: PetscFunctionBegin;
548: PetscAssertPointer(fn,2);
549: if (!mfn->fn) {
550: PetscCall(FNCreate(PetscObjectComm((PetscObject)mfn),&mfn->fn));
551: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mfn->fn,(PetscObject)mfn,0));
552: PetscCall(PetscObjectSetOptions((PetscObject)mfn->fn,((PetscObject)mfn)->options));
553: }
554: *fn = mfn->fn;
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }