Actual source code: lmebasic.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 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: /*@C
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;
60: const char *eqname[] = {
61: "continuous-time Lyapunov",
62: "continuous-time Sylvester",
63: "generalized Lyapunov",
64: "generalized Sylvester",
65: "Stein",
66: "discrete-time Lyapunov"
67: };
69: PetscFunctionBegin;
71: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer));
73: PetscCheckSameComm(lme,1,viewer,2);
75: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
76: if (isascii) {
77: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer));
78: PetscCall(PetscViewerASCIIPushTab(viewer));
79: PetscTryTypeMethod(lme,view,viewer);
80: PetscCall(PetscViewerASCIIPopTab(viewer));
81: PetscCall(PetscViewerASCIIPrintf(viewer," equation type: %s\n",eqname[lme->problem_type]));
82: PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",lme->ncv));
83: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",lme->max_it));
84: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol));
85: } else PetscTryTypeMethod(lme,view,viewer);
86: PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
87: if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
88: PetscCall(BVView(lme->V,viewer));
89: PetscCall(PetscViewerPopFormat(viewer));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*@C
94: LMEViewFromOptions - View from options
96: Collective
98: Input Parameters:
99: + lme - the linear matrix equation context
100: . obj - optional object
101: - name - command line option
103: Level: intermediate
105: .seealso: LMEView(), LMECreate()
106: @*/
107: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
108: {
109: PetscFunctionBegin;
111: PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
114: /*@C
115: LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.
117: Collective
119: Input Parameters:
120: + lme - the linear matrix equation context
121: - viewer - the viewer to display the reason
123: Options Database Keys:
124: . -lme_converged_reason - print reason for convergence, and number of iterations
126: Note:
127: To change the format of the output call PetscViewerPushFormat(viewer,format) before
128: this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
129: display a reason if it fails. The latter can be set in the command line with
130: -lme_converged_reason ::failed
132: Level: intermediate
134: .seealso: LMESetTolerances(), LMEGetIterationNumber(), LMEConvergedReasonViewFromOptions()
135: @*/
136: PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
137: {
138: PetscBool isAscii;
139: PetscViewerFormat format;
141: PetscFunctionBegin;
142: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
143: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
144: if (isAscii) {
145: PetscCall(PetscViewerGetFormat(viewer,&format));
146: PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel));
147: 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));
148: 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));
149: PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel));
150: }
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*@
155: LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
156: the LME converged reason is to be viewed.
158: Collective
160: Input Parameter:
161: . lme - the linear matrix equation context
163: Level: developer
165: .seealso: LMEConvergedReasonView()
166: @*/
167: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
168: {
169: PetscViewer viewer;
170: PetscBool flg;
171: static PetscBool incall = PETSC_FALSE;
172: PetscViewerFormat format;
174: PetscFunctionBegin;
175: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
176: incall = PETSC_TRUE;
177: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
178: if (flg) {
179: PetscCall(PetscViewerPushFormat(viewer,format));
180: PetscCall(LMEConvergedReasonView(lme,viewer));
181: PetscCall(PetscViewerPopFormat(viewer));
182: PetscCall(PetscOptionsRestoreViewer(&viewer));
183: }
184: incall = PETSC_FALSE;
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@
189: LMECreate - Creates the default LME context.
191: Collective
193: Input Parameter:
194: . comm - MPI communicator
196: Output Parameter:
197: . outlme - location to put the LME context
199: Note:
200: The default LME type is LMEKRYLOV
202: Level: beginner
204: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
205: @*/
206: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
207: {
208: LME lme;
210: PetscFunctionBegin;
211: PetscAssertPointer(outlme,2);
212: *outlme = NULL;
213: PetscCall(LMEInitializePackage());
214: PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));
216: lme->A = NULL;
217: lme->B = NULL;
218: lme->D = NULL;
219: lme->E = NULL;
220: lme->C = NULL;
221: lme->X = NULL;
222: lme->problem_type = LME_LYAPUNOV;
223: lme->max_it = PETSC_DEFAULT;
224: lme->ncv = PETSC_DEFAULT;
225: lme->tol = PETSC_DEFAULT;
226: lme->errorifnotconverged = PETSC_FALSE;
228: lme->numbermonitors = 0;
230: lme->V = NULL;
231: lme->nwork = 0;
232: lme->work = NULL;
233: lme->data = NULL;
235: lme->its = 0;
236: lme->errest = 0;
237: lme->setupcalled = 0;
238: lme->reason = LME_CONVERGED_ITERATING;
240: *outlme = lme;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*@C
245: LMESetType - Selects the particular solver to be used in the LME object.
247: Logically Collective
249: Input Parameters:
250: + lme - the linear matrix equation context
251: - type - a known method
253: Options Database Key:
254: . -lme_type <method> - Sets the method; use -help for a list
255: of available methods
257: Notes:
258: See "slepc/include/slepclme.h" for available methods. The default
259: is LMEKRYLOV
261: Normally, it is best to use the LMESetFromOptions() command and
262: then set the LME type from the options database rather than by using
263: this routine. Using the options database provides the user with
264: maximum flexibility in evaluating the different available methods.
265: The LMESetType() routine is provided for those situations where it
266: is necessary to set the iterative solver independently of the command
267: line or options database.
269: Level: intermediate
271: .seealso: LMEType
272: @*/
273: PetscErrorCode LMESetType(LME lme,LMEType type)
274: {
275: PetscErrorCode (*r)(LME);
276: PetscBool match;
278: PetscFunctionBegin;
280: PetscAssertPointer(type,2);
282: PetscCall(PetscObjectTypeCompare((PetscObject)lme,type,&match));
283: if (match) PetscFunctionReturn(PETSC_SUCCESS);
285: PetscCall(PetscFunctionListFind(LMEList,type,&r));
286: PetscCheck(r,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
288: PetscTryTypeMethod(lme,destroy);
289: PetscCall(PetscMemzero(lme->ops,sizeof(struct _LMEOps)));
291: lme->setupcalled = 0;
292: PetscCall(PetscObjectChangeTypeName((PetscObject)lme,type));
293: PetscCall((*r)(lme));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@C
298: LMEGetType - Gets the LME type as a string from the LME object.
300: Not Collective
302: Input Parameter:
303: . lme - the linear matrix equation context
305: Output Parameter:
306: . type - name of LME method
308: Level: intermediate
310: .seealso: LMESetType()
311: @*/
312: PetscErrorCode LMEGetType(LME lme,LMEType *type)
313: {
314: PetscFunctionBegin;
316: PetscAssertPointer(type,2);
317: *type = ((PetscObject)lme)->type_name;
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@C
322: LMERegister - Adds a method to the linear matrix equation solver package.
324: Not Collective
326: Input Parameters:
327: + name - name of a new user-defined solver
328: - function - routine to create the solver context
330: Notes:
331: LMERegister() may be called multiple times to add several user-defined solvers.
333: Example Usage:
334: .vb
335: LMERegister("my_solver",MySolverCreate);
336: .ve
338: Then, your solver can be chosen with the procedural interface via
339: $ LMESetType(lme,"my_solver")
340: or at runtime via the option
341: $ -lme_type my_solver
343: Level: advanced
345: .seealso: LMERegisterAll()
346: @*/
347: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
348: {
349: PetscFunctionBegin;
350: PetscCall(LMEInitializePackage());
351: PetscCall(PetscFunctionListAdd(&LMEList,name,function));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@C
356: LMEMonitorRegister - Adds LME monitor routine.
358: Not Collective
360: Input Parameters:
361: + name - name of a new monitor routine
362: . vtype - a PetscViewerType for the output
363: . format - a PetscViewerFormat for the output
364: . monitor - monitor routine
365: . create - creation routine, or NULL
366: - destroy - destruction routine, or NULL
368: Notes:
369: LMEMonitorRegister() may be called multiple times to add several user-defined monitors.
371: Example Usage:
372: .vb
373: LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
374: .ve
376: Then, your monitor can be chosen with the procedural interface via
377: $ LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
378: or at runtime via the option
379: $ -lme_monitor_my_monitor
381: Level: advanced
383: .seealso: LMEMonitorRegisterAll()
384: @*/
385: PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(LME,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
386: {
387: char key[PETSC_MAX_PATH_LEN];
389: PetscFunctionBegin;
390: PetscCall(LMEInitializePackage());
391: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
392: PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
393: if (create) PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
394: if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: LMEReset - Resets the LME context to the initial state (prior to setup)
400: and destroys any allocated Vecs and Mats.
402: Collective
404: Input Parameter:
405: . lme - linear matrix equation context obtained from LMECreate()
407: Level: advanced
409: .seealso: LMEDestroy()
410: @*/
411: PetscErrorCode LMEReset(LME lme)
412: {
413: PetscFunctionBegin;
415: if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
416: PetscTryTypeMethod(lme,reset);
417: PetscCall(MatDestroy(&lme->A));
418: PetscCall(MatDestroy(&lme->B));
419: PetscCall(MatDestroy(&lme->D));
420: PetscCall(MatDestroy(&lme->E));
421: PetscCall(MatDestroy(&lme->C));
422: PetscCall(MatDestroy(&lme->X));
423: PetscCall(BVDestroy(&lme->V));
424: PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
425: lme->nwork = 0;
426: lme->setupcalled = 0;
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: /*@C
431: LMEDestroy - Destroys the LME context.
433: Collective
435: Input Parameter:
436: . lme - linear matrix equation context obtained from LMECreate()
438: Level: beginner
440: .seealso: LMECreate(), LMESetUp(), LMESolve()
441: @*/
442: PetscErrorCode LMEDestroy(LME *lme)
443: {
444: PetscFunctionBegin;
445: if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
447: if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
448: PetscCall(LMEReset(*lme));
449: PetscTryTypeMethod(*lme,destroy);
450: PetscCall(LMEMonitorCancel(*lme));
451: PetscCall(PetscHeaderDestroy(lme));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: /*@
456: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
458: Collective
460: Input Parameters:
461: + lme - linear matrix equation context obtained from LMECreate()
462: - bv - the basis vectors object
464: Note:
465: Use LMEGetBV() to retrieve the basis vectors context (for example,
466: to free it at the end of the computations).
468: Level: advanced
470: .seealso: LMEGetBV()
471: @*/
472: PetscErrorCode LMESetBV(LME lme,BV bv)
473: {
474: PetscFunctionBegin;
477: PetscCheckSameComm(lme,1,bv,2);
478: PetscCall(PetscObjectReference((PetscObject)bv));
479: PetscCall(BVDestroy(&lme->V));
480: lme->V = bv;
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*@
485: LMEGetBV - Obtain the basis vectors object associated to the matrix
486: function solver.
488: Not Collective
490: Input Parameters:
491: . lme - linear matrix equation context obtained from LMECreate()
493: Output Parameter:
494: . bv - basis vectors context
496: Level: advanced
498: .seealso: LMESetBV()
499: @*/
500: PetscErrorCode LMEGetBV(LME lme,BV *bv)
501: {
502: PetscFunctionBegin;
504: PetscAssertPointer(bv,2);
505: if (!lme->V) {
506: PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
507: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
508: PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
509: }
510: *bv = lme->V;
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }