Actual source code: mfnbasic.c
slepc-3.11.2 2019-07-30
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: }