Actual source code: lmebasic.c
slepc-3.14.2 2021-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: PetscFunctionList LMEList = 0;
17: PetscBool LMERegisterAllCalled = PETSC_FALSE;
18: PetscClassId LME_CLASSID = 0;
19: PetscLogEvent LME_SetUp = 0,LME_Solve = 0,LME_ComputeError = 0;
21: /*@C
22: LMEView - Prints the LME data structure.
24: Collective on lme
26: Input Parameters:
27: + lme - the linear matrix equation solver context
28: - viewer - optional visualization context
30: Options Database Key:
31: . -lme_view - Calls LMEView() at end of LMESolve()
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
45: @*/
46: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
47: {
49: PetscBool isascii;
50: const char *eqname[] = {
51: "continuous-time Lyapunov",
52: "continuous-time Sylvester",
53: "generalized Lyapunov",
54: "generalized Sylvester",
55: "Stein",
56: "discrete-time Lyapunov"
57: };
61: if (!viewer) {
62: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer);
63: }
67: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
68: if (isascii) {
69: PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer);
70: if (lme->ops->view) {
71: PetscViewerASCIIPushTab(viewer);
72: (*lme->ops->view)(lme,viewer);
73: PetscViewerASCIIPopTab(viewer);
74: }
75: PetscViewerASCIIPrintf(viewer," equation type: %s\n",eqname[lme->problem_type]);
76: PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %D\n",lme->ncv);
77: PetscViewerASCIIPrintf(viewer," maximum number of iterations: %D\n",lme->max_it);
78: PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol);
79: } else {
80: if (lme->ops->view) {
81: (*lme->ops->view)(lme,viewer);
82: }
83: }
84: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
85: if (!lme->V) { LMEGetBV(lme,&lme->V); }
86: BVView(lme->V,viewer);
87: PetscViewerPopFormat(viewer);
88: return(0);
89: }
91: /*@C
92: LMEViewFromOptions - View from options
94: Collective on LME
96: Input Parameters:
97: + lme - the linear matrix equation context
98: . obj - optional object
99: - name - command line option
101: Level: intermediate
103: .seealso: LMEView(), LMECreate()
104: @*/
105: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
106: {
111: PetscObjectViewFromOptions((PetscObject)lme,obj,name);
112: return(0);
113: }
114: /*@C
115: LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.
117: Collective on lme
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: PetscErrorCode ierr;
139: PetscBool isAscii;
140: PetscViewerFormat format;
143: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
144: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
145: if (isAscii) {
146: PetscViewerGetFormat(viewer,&format);
147: PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
148: if (lme->reason > 0 && format != PETSC_VIEWER_FAILED) {
149: PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve converged due to %s; iterations %D\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its);
150: } else if (lme->reason <= 0) {
151: PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve did not converge due to %s; iterations %D\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its);
152: }
153: PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
154: }
155: return(0);
156: }
158: /*@
159: LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
160: the LME converged reason is to be viewed.
162: Collective on lme
164: Input Parameter:
165: . lme - the linear matrix equation context
167: Level: developer
169: .seealso: LMEConvergedReasonView()
170: @*/
171: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
172: {
173: PetscErrorCode ierr;
174: PetscViewer viewer;
175: PetscBool flg;
176: static PetscBool incall = PETSC_FALSE;
177: PetscViewerFormat format;
180: if (incall) return(0);
181: incall = PETSC_TRUE;
182: PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
183: if (flg) {
184: PetscViewerPushFormat(viewer,format);
185: LMEConvergedReasonView(lme,viewer);
186: PetscViewerPopFormat(viewer);
187: PetscViewerDestroy(&viewer);
188: }
189: incall = PETSC_FALSE;
190: return(0);
191: }
193: /*@
194: LMECreate - Creates the default LME context.
196: Collective
198: Input Parameter:
199: . comm - MPI communicator
201: Output Parameter:
202: . lme - location to put the LME context
204: Note:
205: The default LME type is LMEKRYLOV
207: Level: beginner
209: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
210: @*/
211: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
212: {
214: LME lme;
218: *outlme = 0;
219: LMEInitializePackage();
220: SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);
222: lme->A = NULL;
223: lme->B = NULL;
224: lme->D = NULL;
225: lme->E = NULL;
226: lme->C = NULL;
227: lme->X = NULL;
228: lme->problem_type = LME_LYAPUNOV;
229: lme->max_it = PETSC_DEFAULT;
230: lme->ncv = PETSC_DEFAULT;
231: lme->tol = PETSC_DEFAULT;
232: lme->errorifnotconverged = PETSC_FALSE;
234: lme->numbermonitors = 0;
236: lme->V = NULL;
237: lme->nwork = 0;
238: lme->work = NULL;
239: lme->data = NULL;
241: lme->its = 0;
242: lme->errest = 0;
243: lme->setupcalled = 0;
244: lme->reason = LME_CONVERGED_ITERATING;
246: *outlme = lme;
247: return(0);
248: }
250: /*@C
251: LMESetType - Selects the particular solver to be used in the LME object.
253: Logically Collective on lme
255: Input Parameters:
256: + lme - the linear matrix equation context
257: - type - a known method
259: Options Database Key:
260: . -lme_type <method> - Sets the method; use -help for a list
261: of available methods
263: Notes:
264: See "slepc/include/slepclme.h" for available methods. The default
265: is LMEKRYLOV
267: Normally, it is best to use the LMESetFromOptions() command and
268: then set the LME type from the options database rather than by using
269: this routine. Using the options database provides the user with
270: maximum flexibility in evaluating the different available methods.
271: The LMESetType() routine is provided for those situations where it
272: is necessary to set the iterative solver independently of the command
273: line or options database.
275: Level: intermediate
277: .seealso: LMEType
278: @*/
279: PetscErrorCode LMESetType(LME lme,LMEType type)
280: {
281: PetscErrorCode ierr,(*r)(LME);
282: PetscBool match;
288: PetscObjectTypeCompare((PetscObject)lme,type,&match);
289: if (match) return(0);
291: PetscFunctionListFind(LMEList,type,&r);
292: if (!r) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
294: if (lme->ops->destroy) { (*lme->ops->destroy)(lme); }
295: PetscMemzero(lme->ops,sizeof(struct _LMEOps));
297: lme->setupcalled = 0;
298: PetscObjectChangeTypeName((PetscObject)lme,type);
299: (*r)(lme);
300: return(0);
301: }
303: /*@C
304: LMEGetType - Gets the LME type as a string from the LME object.
306: Not Collective
308: Input Parameter:
309: . lme - the linear matrix equation context
311: Output Parameter:
312: . name - name of LME method
314: Level: intermediate
316: .seealso: LMESetType()
317: @*/
318: PetscErrorCode LMEGetType(LME lme,LMEType *type)
319: {
323: *type = ((PetscObject)lme)->type_name;
324: return(0);
325: }
327: /*@C
328: LMERegister - Adds a method to the linear matrix equation solver package.
330: Not Collective
332: Input Parameters:
333: + name - name of a new user-defined solver
334: - function - routine to create the solver context
336: Notes:
337: LMERegister() may be called multiple times to add several user-defined solvers.
339: Sample usage:
340: .vb
341: LMERegister("my_solver",MySolverCreate);
342: .ve
344: Then, your solver can be chosen with the procedural interface via
345: $ LMESetType(lme,"my_solver")
346: or at runtime via the option
347: $ -lme_type my_solver
349: Level: advanced
351: .seealso: LMERegisterAll()
352: @*/
353: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
354: {
358: LMEInitializePackage();
359: PetscFunctionListAdd(&LMEList,name,function);
360: return(0);
361: }
363: /*@
364: LMEReset - Resets the LME context to the initial state (prior to setup)
365: and destroys any allocated Vecs and Mats.
367: Collective on lme
369: Input Parameter:
370: . lme - linear matrix equation context obtained from LMECreate()
372: Level: advanced
374: .seealso: LMEDestroy()
375: @*/
376: PetscErrorCode LMEReset(LME lme)
377: {
382: if (!lme) return(0);
383: if (lme->ops->reset) { (lme->ops->reset)(lme); }
384: MatDestroy(&lme->A);
385: MatDestroy(&lme->B);
386: MatDestroy(&lme->D);
387: MatDestroy(&lme->E);
388: MatDestroy(&lme->C);
389: MatDestroy(&lme->X);
390: BVDestroy(&lme->V);
391: VecDestroyVecs(lme->nwork,&lme->work);
392: lme->nwork = 0;
393: lme->setupcalled = 0;
394: return(0);
395: }
397: /*@
398: LMEDestroy - Destroys the LME context.
400: Collective on lme
402: Input Parameter:
403: . lme - linear matrix equation context obtained from LMECreate()
405: Level: beginner
407: .seealso: LMECreate(), LMESetUp(), LMESolve()
408: @*/
409: PetscErrorCode LMEDestroy(LME *lme)
410: {
414: if (!*lme) return(0);
416: if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
417: LMEReset(*lme);
418: if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
419: LMEMonitorCancel(*lme);
420: PetscHeaderDestroy(lme);
421: return(0);
422: }
424: /*@
425: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
427: Collective on lme
429: Input Parameters:
430: + lme - linear matrix equation context obtained from LMECreate()
431: - bv - the basis vectors object
433: Note:
434: Use LMEGetBV() to retrieve the basis vectors context (for example,
435: to free it at the end of the computations).
437: Level: advanced
439: .seealso: LMEGetBV()
440: @*/
441: PetscErrorCode LMESetBV(LME lme,BV bv)
442: {
449: PetscObjectReference((PetscObject)bv);
450: BVDestroy(&lme->V);
451: lme->V = bv;
452: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
453: return(0);
454: }
456: /*@
457: LMEGetBV - Obtain the basis vectors object associated to the matrix
458: function solver.
460: Not Collective
462: Input Parameters:
463: . lme - linear matrix equation context obtained from LMECreate()
465: Output Parameter:
466: . bv - basis vectors context
468: Level: advanced
470: .seealso: LMESetBV()
471: @*/
472: PetscErrorCode LMEGetBV(LME lme,BV *bv)
473: {
479: if (!lme->V) {
480: BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
481: PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
482: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
483: PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
484: }
485: *bv = lme->V;
486: return(0);
487: }