Actual source code: lmebasic.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 LME routines
12: */
14: #include <slepc/private/lmeimpl.h>
16: /* Logging support */
17: PetscClassId LME_CLASSID = 0;
18: PetscLogEvent LME_SetUp = 0,LME_Solve = 0,LME_ComputeError = 0;
20: /* List of registered LME routines */
21: PetscFunctionList LMEList = NULL;
22: PetscBool LMERegisterAllCalled = PETSC_FALSE;
24: /* List of registered LME monitors */
25: PetscFunctionList LMEMonitorList = NULL;
26: PetscFunctionList LMEMonitorCreateList = NULL;
27: PetscFunctionList LMEMonitorDestroyList = NULL;
28: PetscBool LMEMonitorRegisterAllCalled = PETSC_FALSE;
30: /*@
31: LMEView - Prints the LME data structure.
33: Collective
35: Input Parameters:
36: + lme - the linear matrix equation solver context
37: - viewer - optional visualization context
39: Options Database Key:
40: . -lme_view - Calls LMEView() at end of LMESolve()
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: `LMECreate()`
56: @*/
57: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
58: {
59: PetscBool isascii;
61: PetscFunctionBegin;
63: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer));
65: PetscCheckSameComm(lme,1,viewer,2);
67: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
68: if (isascii) {
69: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer));
70: PetscCall(PetscViewerASCIIPushTab(viewer));
71: PetscTryTypeMethod(lme,view,viewer);
72: PetscCall(PetscViewerASCIIPopTab(viewer));
73: PetscCall(PetscViewerASCIIPrintf(viewer," equation type: %s\n",LMEProblemTypes[lme->problem_type]));
74: PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",lme->ncv));
75: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",lme->max_it));
76: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol));
77: } else PetscTryTypeMethod(lme,view,viewer);
78: PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
79: if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
80: PetscCall(BVView(lme->V,viewer));
81: PetscCall(PetscViewerPopFormat(viewer));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: LMEViewFromOptions - View from options
88: Collective
90: Input Parameters:
91: + lme - the linear matrix equation context
92: . obj - optional object
93: - name - command line option
95: Level: intermediate
97: .seealso: `LMEView()`, `LMECreate()`
98: @*/
99: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
100: {
101: PetscFunctionBegin;
103: PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
106: /*@
107: LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.
109: Collective
111: Input Parameters:
112: + lme - the linear matrix equation context
113: - viewer - the viewer to display the reason
115: Options Database Keys:
116: . -lme_converged_reason - print reason for convergence, and number of iterations
118: Note:
119: To change the format of the output call PetscViewerPushFormat(viewer,format) before
120: this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
121: display a reason if it fails. The latter can be set in the command line with
122: -lme_converged_reason ::failed
124: Level: intermediate
126: .seealso: `LMESetTolerances()`, `LMEGetIterationNumber()`, `LMEConvergedReasonViewFromOptions()`
127: @*/
128: PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
129: {
130: PetscBool isAscii;
131: PetscViewerFormat format;
133: PetscFunctionBegin;
134: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
135: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
136: if (isAscii) {
137: PetscCall(PetscViewerGetFormat(viewer,&format));
138: PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel));
139: if (lme->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve converged due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its));
140: else if (lme->reason <= 0) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its));
141: PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel));
142: }
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@
147: LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
148: the LME converged reason is to be viewed.
150: Collective
152: Input Parameter:
153: . lme - the linear matrix equation context
155: Level: developer
157: .seealso: `LMEConvergedReasonView()`
158: @*/
159: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
160: {
161: PetscViewer viewer;
162: PetscBool flg;
163: static PetscBool incall = PETSC_FALSE;
164: PetscViewerFormat format;
166: PetscFunctionBegin;
167: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
168: incall = PETSC_TRUE;
169: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
170: if (flg) {
171: PetscCall(PetscViewerPushFormat(viewer,format));
172: PetscCall(LMEConvergedReasonView(lme,viewer));
173: PetscCall(PetscViewerPopFormat(viewer));
174: PetscCall(PetscViewerDestroy(&viewer));
175: }
176: incall = PETSC_FALSE;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: LMECreate - Creates the default LME context.
183: Collective
185: Input Parameter:
186: . comm - MPI communicator
188: Output Parameter:
189: . outlme - location to put the LME context
191: Note:
192: The default LME type is LMEKRYLOV
194: Level: beginner
196: .seealso: `LMESetUp()`, `LMESolve()`, `LMEDestroy()`, `LME`
197: @*/
198: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
199: {
200: LME lme;
202: PetscFunctionBegin;
203: PetscAssertPointer(outlme,2);
204: PetscCall(LMEInitializePackage());
205: PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));
207: lme->A = NULL;
208: lme->B = NULL;
209: lme->D = NULL;
210: lme->E = NULL;
211: lme->C = NULL;
212: lme->X = NULL;
213: lme->problem_type = LME_LYAPUNOV;
214: lme->max_it = PETSC_DETERMINE;
215: lme->ncv = PETSC_DETERMINE;
216: lme->tol = PETSC_DETERMINE;
217: lme->errorifnotconverged = PETSC_FALSE;
219: lme->numbermonitors = 0;
221: lme->V = NULL;
222: lme->nwork = 0;
223: lme->work = NULL;
224: lme->data = NULL;
226: lme->its = 0;
227: lme->errest = 0;
228: lme->setupcalled = 0;
229: lme->reason = LME_CONVERGED_ITERATING;
231: *outlme = lme;
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@
236: LMESetType - Selects the particular solver to be used in the LME object.
238: Logically Collective
240: Input Parameters:
241: + lme - the linear matrix equation context
242: - type - a known method
244: Options Database Key:
245: . -lme_type <method> - Sets the method; use -help for a list
246: of available methods
248: Notes:
249: See "slepc/include/slepclme.h" for available methods. The default
250: is LMEKRYLOV
252: Normally, it is best to use the LMESetFromOptions() command and
253: then set the LME type from the options database rather than by using
254: this routine. Using the options database provides the user with
255: maximum flexibility in evaluating the different available methods.
256: The LMESetType() routine is provided for those situations where it
257: is necessary to set the iterative solver independently of the command
258: line or options database.
260: Level: intermediate
262: .seealso: `LMEType`
263: @*/
264: PetscErrorCode LMESetType(LME lme,LMEType type)
265: {
266: PetscErrorCode (*r)(LME);
267: PetscBool match;
269: PetscFunctionBegin;
271: PetscAssertPointer(type,2);
273: PetscCall(PetscObjectTypeCompare((PetscObject)lme,type,&match));
274: if (match) PetscFunctionReturn(PETSC_SUCCESS);
276: PetscCall(PetscFunctionListFind(LMEList,type,&r));
277: PetscCheck(r,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
279: PetscTryTypeMethod(lme,destroy);
280: PetscCall(PetscMemzero(lme->ops,sizeof(struct _LMEOps)));
282: lme->setupcalled = 0;
283: PetscCall(PetscObjectChangeTypeName((PetscObject)lme,type));
284: PetscCall((*r)(lme));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: LMEGetType - Gets the LME type as a string from the LME object.
291: Not Collective
293: Input Parameter:
294: . lme - the linear matrix equation context
296: Output Parameter:
297: . type - name of LME method
299: Level: intermediate
301: .seealso: `LMESetType()`
302: @*/
303: PetscErrorCode LMEGetType(LME lme,LMEType *type)
304: {
305: PetscFunctionBegin;
307: PetscAssertPointer(type,2);
308: *type = ((PetscObject)lme)->type_name;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: LMERegister - Adds a method to the linear matrix equation solver package.
315: Not Collective
317: Input Parameters:
318: + name - name of a new user-defined solver
319: - function - routine to create the solver context
321: Notes:
322: LMERegister() may be called multiple times to add several user-defined solvers.
324: Example Usage:
325: .vb
326: LMERegister("my_solver",MySolverCreate);
327: .ve
329: Then, your solver can be chosen with the procedural interface via
330: $ LMESetType(lme,"my_solver")
331: or at runtime via the option
332: $ -lme_type my_solver
334: Level: advanced
336: .seealso: `LMERegisterAll()`
337: @*/
338: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
339: {
340: PetscFunctionBegin;
341: PetscCall(LMEInitializePackage());
342: PetscCall(PetscFunctionListAdd(&LMEList,name,function));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@C
347: LMEMonitorRegister - Registers an LME monitor routine that may be accessed with LMEMonitorSetFromOptions().
349: Not Collective
351: Input Parameters:
352: + name - name of a new monitor routine
353: . vtype - a PetscViewerType for the output
354: . format - a PetscViewerFormat for the output
355: . monitor - monitor routine, see LMEMonitorRegisterFn
356: . create - creation routine, or NULL
357: - destroy - destruction routine, or NULL
359: Notes:
360: LMEMonitorRegister() may be called multiple times to add several user-defined monitors.
362: The calling sequence for the given function matches the calling sequence of LMEMonitorFn
363: functions passed to LMEMonitorSet() with the additional requirement that its final argument
364: be a PetscViewerAndFormat.
366: Example Usage:
367: .vb
368: LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
369: .ve
371: Then, your monitor can be chosen with the procedural interface via
372: $ LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
373: or at runtime via the option
374: $ -lme_monitor_my_monitor
376: Level: advanced
378: .seealso: `LMEMonitorSet()`, `LMEMonitorRegisterAll()`
379: @*/
380: PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,LMEMonitorRegisterFn *monitor,LMEMonitorRegisterCreateFn *create,LMEMonitorRegisterDestroyFn *destroy)
381: {
382: char key[PETSC_MAX_PATH_LEN];
384: PetscFunctionBegin;
385: PetscCall(LMEInitializePackage());
386: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
387: PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
388: if (create) PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
389: if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*@
394: LMEReset - Resets the LME context to the initial state (prior to setup)
395: and destroys any allocated Vecs and Mats.
397: Collective
399: Input Parameter:
400: . lme - linear matrix equation context obtained from LMECreate()
402: Level: advanced
404: .seealso: `LMEDestroy()`
405: @*/
406: PetscErrorCode LMEReset(LME lme)
407: {
408: PetscFunctionBegin;
410: if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
411: PetscTryTypeMethod(lme,reset);
412: PetscCall(MatDestroy(&lme->A));
413: PetscCall(MatDestroy(&lme->B));
414: PetscCall(MatDestroy(&lme->D));
415: PetscCall(MatDestroy(&lme->E));
416: PetscCall(MatDestroy(&lme->C));
417: PetscCall(MatDestroy(&lme->X));
418: PetscCall(BVDestroy(&lme->V));
419: PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
420: lme->nwork = 0;
421: lme->setupcalled = 0;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@
426: LMEDestroy - Destroys the LME context.
428: Collective
430: Input Parameter:
431: . lme - linear matrix equation context obtained from LMECreate()
433: Level: beginner
435: .seealso: `LMECreate()`, `LMESetUp()`, `LMESolve()`
436: @*/
437: PetscErrorCode LMEDestroy(LME *lme)
438: {
439: PetscFunctionBegin;
440: if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
442: if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
443: PetscCall(LMEReset(*lme));
444: PetscTryTypeMethod(*lme,destroy);
445: PetscCall(LMEMonitorCancel(*lme));
446: PetscCall(PetscHeaderDestroy(lme));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
453: Collective
455: Input Parameters:
456: + lme - linear matrix equation context obtained from LMECreate()
457: - bv - the basis vectors object
459: Note:
460: Use LMEGetBV() to retrieve the basis vectors context (for example,
461: to free it at the end of the computations).
463: Level: advanced
465: .seealso: `LMEGetBV()`
466: @*/
467: PetscErrorCode LMESetBV(LME lme,BV bv)
468: {
469: PetscFunctionBegin;
472: PetscCheckSameComm(lme,1,bv,2);
473: PetscCall(PetscObjectReference((PetscObject)bv));
474: PetscCall(BVDestroy(&lme->V));
475: lme->V = bv;
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: LMEGetBV - Obtain the basis vectors object associated to the matrix
481: function solver.
483: Not Collective
485: Input Parameters:
486: . lme - linear matrix equation context obtained from LMECreate()
488: Output Parameter:
489: . bv - basis vectors context
491: Level: advanced
493: .seealso: `LMESetBV()`
494: @*/
495: PetscErrorCode LMEGetBV(LME lme,BV *bv)
496: {
497: PetscFunctionBegin;
499: PetscAssertPointer(bv,2);
500: if (!lme->V) {
501: PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
502: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
503: PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
504: }
505: *bv = lme->V;
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }