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: Notes:
43: The available visualization contexts include
44: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
45: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
46: first process opens the file; all other processes send their data to the
47: first one to print
49: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
50: to output to a specified file.
52: Level: beginner
54: .seealso: [](ch:lme), `LMECreate()`, `LMEViewFromOptions()`
55: @*/
56: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
57: {
58: PetscBool isascii;
60: PetscFunctionBegin;
62: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer));
64: PetscCheckSameComm(lme,1,viewer,2);
66: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
67: if (isascii) {
68: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer));
69: PetscCall(PetscViewerASCIIPushTab(viewer));
70: PetscTryTypeMethod(lme,view,viewer);
71: PetscCall(PetscViewerASCIIPopTab(viewer));
72: PetscCall(PetscViewerASCIIPrintf(viewer," equation type: %s\n",LMEProblemTypes[lme->problem_type]));
73: PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",lme->ncv));
74: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",lme->max_it));
75: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol));
76: } else PetscTryTypeMethod(lme,view,viewer);
77: PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
78: if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
79: PetscCall(BVView(lme->V,viewer));
80: PetscCall(PetscViewerPopFormat(viewer));
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: /*@
85: LMEViewFromOptions - View (print) an `LME` object based on values in the options database.
87: Collective
89: Input Parameters:
90: + lme - the linear matrix equation solver context
91: . obj - optional object that provides the options prefix used to query the options database
92: - name - command line option
94: Level: intermediate
96: .seealso: [](ch:lme), `LMEView()`, `LMECreate()`, `PetscObjectViewFromOptions()`
97: @*/
98: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
99: {
100: PetscFunctionBegin;
102: PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
105: /*@
106: LMEConvergedReasonView - Displays the reason an `LME` solve converged or diverged.
108: Collective
110: Input Parameters:
111: + lme - the linear matrix equation solver context
112: - viewer - the viewer to display the reason
114: Options Database Key:
115: . -lme_converged_reason - print reason for convergence/divergence, and number of iterations
117: Notes:
118: Use `LMEConvergedReasonViewFromOptions()` to display the reason based on values
119: in the options database.
121: To change the format of the output call `PetscViewerPushFormat()` before this
122: call. Use `PETSC_VIEWER_DEFAULT` for the default, or `PETSC_VIEWER_FAILED` to only
123: display a reason if it fails. The latter can be set in the command line with
124: `-lme_converged_reason ::failed`.
126: Level: intermediate
128: .seealso: [](ch:lme), `LMESetTolerances()`, `LMEGetIterationNumber()`, `LMEConvergedReasonViewFromOptions()`
129: @*/
130: PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
131: {
132: PetscBool isAscii;
133: PetscViewerFormat format;
135: PetscFunctionBegin;
136: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
137: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
138: if (isAscii) {
139: PetscCall(PetscViewerGetFormat(viewer,&format));
140: PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel));
141: 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));
142: 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));
143: PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel));
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*@
149: LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
150: the `LME` converged reason is to be viewed.
152: Collective
154: Input Parameter:
155: . lme - the linear matrix equation solver context
157: Level: intermediate
159: .seealso: [](ch:lme), `LMEConvergedReasonView()`
160: @*/
161: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
162: {
163: PetscViewer viewer;
164: PetscBool flg;
165: static PetscBool incall = PETSC_FALSE;
166: PetscViewerFormat format;
168: PetscFunctionBegin;
169: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
170: incall = PETSC_TRUE;
171: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
172: if (flg) {
173: PetscCall(PetscViewerPushFormat(viewer,format));
174: PetscCall(LMEConvergedReasonView(lme,viewer));
175: PetscCall(PetscViewerPopFormat(viewer));
176: PetscCall(PetscViewerDestroy(&viewer));
177: }
178: incall = PETSC_FALSE;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*@
183: LMECreate - Creates the `LME` context.
185: Collective
187: Input Parameter:
188: . comm - MPI communicator
190: Output Parameter:
191: . outlme - location to put the `LME` context
193: Note:
194: The default `LME` type is `LMEKRYLOV`.
196: Level: beginner
198: .seealso: [](ch:lme), `LMESetUp()`, `LMESolve()`, `LMEDestroy()`, `LME`
199: @*/
200: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
201: {
202: LME lme;
204: PetscFunctionBegin;
205: PetscAssertPointer(outlme,2);
206: PetscCall(LMEInitializePackage());
207: PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));
209: lme->A = NULL;
210: lme->B = NULL;
211: lme->D = NULL;
212: lme->E = NULL;
213: lme->C = NULL;
214: lme->X = NULL;
215: lme->problem_type = LME_LYAPUNOV;
216: lme->max_it = PETSC_DETERMINE;
217: lme->ncv = PETSC_DETERMINE;
218: lme->tol = PETSC_DETERMINE;
219: lme->errorifnotconverged = PETSC_FALSE;
221: lme->numbermonitors = 0;
223: lme->V = NULL;
224: lme->nwork = 0;
225: lme->work = NULL;
226: lme->data = NULL;
228: lme->its = 0;
229: lme->errest = 0;
230: lme->setupcalled = 0;
231: lme->reason = LME_CONVERGED_ITERATING;
233: *outlme = lme;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*@
238: LMESetType - Selects the particular solver to be used in the `LME` object.
240: Logically Collective
242: Input Parameters:
243: + lme - the linear matrix equation solver context
244: - type - a known method
246: Options Database Key:
247: . -lme_type type - sets the method; use -help for a list of available methods
249: Notes:
250: See `LMEType` for available methods. The default 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: [](ch:lme), `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 solver context
296: Output Parameter:
297: . type - name of `LME` method
299: Note:
300: `type` should not be retained for later use as it will be an invalid pointer
301: if the `LMEType` of `lme` is changed.
303: Level: intermediate
305: .seealso: [](ch:lme), `LMESetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
306: @*/
307: PetscErrorCode LMEGetType(LME lme,LMEType *type)
308: {
309: PetscFunctionBegin;
311: PetscAssertPointer(type,2);
312: *type = ((PetscObject)lme)->type_name;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@C
317: LMERegister - Adds a method to the linear matrix equation solver package.
319: Not Collective
321: Input Parameters:
322: + name - name of a new user-defined solver
323: - function - routine to create the solver context
325: Note:
326: `LMERegister()` may be called multiple times to add several user-defined solvers.
328: Example Usage:
329: .vb
330: LMERegister("my_solver",MySolverCreate);
331: .ve
333: Then, your solver can be chosen with the procedural interface via
334: .vb
335: LMESetType(lme,"my_solver")
336: .ve
337: or at runtime via the option `-lme_type my_solver`.
339: Level: advanced
341: .seealso: [](ch:lme), `LMERegisterAll()`
342: @*/
343: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
344: {
345: PetscFunctionBegin;
346: PetscCall(LMEInitializePackage());
347: PetscCall(PetscFunctionListAdd(&LMEList,name,function));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@C
352: LMEMonitorRegister - Registers an `LME` monitor routine that may be accessed with
353: `LMEMonitorSetFromOptions()`.
355: Not Collective
357: Input Parameters:
358: + name - name of a new monitor routine
359: . vtype - a `PetscViewerType` for the output
360: . format - a `PetscViewerFormat` for the output
361: . monitor - monitor routine, see `LMEMonitorRegisterFn`
362: . create - creation routine, or `NULL`
363: - destroy - destruction routine, or `NULL`
365: Notes:
366: `LMEMonitorRegister()` may be called multiple times to add several user-defined monitors.
368: The calling sequence for the given function matches the calling sequence of `LMEMonitorFn`
369: functions passed to `LMEMonitorSet()` with the additional requirement that its final argument
370: be a `PetscViewerAndFormat`.
372: Example Usage:
373: .vb
374: LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
375: .ve
377: Then, your monitor can be chosen with the procedural interface via
378: .vb
379: LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL);
380: .ve
381: or at runtime via the option `-lme_monitor_my_monitor`.
383: Level: advanced
385: .seealso: [](ch:lme), `LMEMonitorSet()`, `LMEMonitorRegisterAll()`, `LMEMonitorSetFromOptions()`
386: @*/
387: PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,LMEMonitorRegisterFn *monitor,LMEMonitorRegisterCreateFn *create,LMEMonitorRegisterDestroyFn *destroy)
388: {
389: char key[PETSC_MAX_PATH_LEN];
391: PetscFunctionBegin;
392: PetscCall(LMEInitializePackage());
393: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
394: PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
395: if (create) PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
396: if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*@
401: LMEReset - Resets the `LME` context to the initial state (prior to setup)
402: and destroys any allocated `Vec`s and `Mat`s.
404: Collective
406: Input Parameter:
407: . lme - the linear matrix equation solver context
409: Level: advanced
411: .seealso: [](ch:lme), `LMEDestroy()`
412: @*/
413: PetscErrorCode LMEReset(LME lme)
414: {
415: PetscFunctionBegin;
417: if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
418: PetscTryTypeMethod(lme,reset);
419: PetscCall(MatDestroy(&lme->A));
420: PetscCall(MatDestroy(&lme->B));
421: PetscCall(MatDestroy(&lme->D));
422: PetscCall(MatDestroy(&lme->E));
423: PetscCall(MatDestroy(&lme->C));
424: PetscCall(MatDestroy(&lme->X));
425: PetscCall(BVDestroy(&lme->V));
426: PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
427: lme->nwork = 0;
428: lme->setupcalled = 0;
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: /*@
433: LMEDestroy - Destroys the `LME` context.
435: Collective
437: Input Parameter:
438: . lme - the linear matrix equation solver context
440: Level: beginner
442: .seealso: [](ch:lme), `LMECreate()`, `LMESetUp()`, `LMESolve()`
443: @*/
444: PetscErrorCode LMEDestroy(LME *lme)
445: {
446: PetscFunctionBegin;
447: if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
449: if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
450: PetscCall(LMEReset(*lme));
451: PetscTryTypeMethod(*lme,destroy);
452: PetscCall(LMEMonitorCancel(*lme));
453: PetscCall(PetscHeaderDestroy(lme));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
460: Collective
462: Input Parameters:
463: + lme - the linear matrix equation solver context
464: - bv - the basis vectors object
466: Note:
467: Use `LMEGetBV()` to retrieve the basis vectors context (for example,
468: to free it at the end of the computations).
470: Level: advanced
472: .seealso: [](ch:lme), `LMEGetBV()`
473: @*/
474: PetscErrorCode LMESetBV(LME lme,BV bv)
475: {
476: PetscFunctionBegin;
479: PetscCheckSameComm(lme,1,bv,2);
480: PetscCall(PetscObjectReference((PetscObject)bv));
481: PetscCall(BVDestroy(&lme->V));
482: lme->V = bv;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /*@
487: LMEGetBV - Obtain the basis vectors object associated to the matrix
488: function solver.
490: Not Collective
492: Input Parameter:
493: . lme - the linear matrix equation solver context
495: Output Parameter:
496: . bv - basis vectors context
498: Level: advanced
500: .seealso: [](ch:lme), `LMESetBV()`
501: @*/
502: PetscErrorCode LMEGetBV(LME lme,BV *bv)
503: {
504: PetscFunctionBegin;
506: PetscAssertPointer(bv,2);
507: if (!lme->V) {
508: PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
509: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
510: PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
511: }
512: *bv = lme->V;
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }