Actual source code: lmebasic.c
slepc-3.9.2 2018-07-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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> /*I "slepclme.h" I*/
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
46: .seealso: PetscViewerASCIIOpen()
47: @*/
48: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
49: {
51: PetscBool isascii;
52: PetscInt tabs;
53: const char *eqname[] = {
54: "continuous-time Lyapunov",
55: "continuous-time Sylvester",
56: "generalized Lyapunov",
57: "generalized Sylvester",
58: "Stein",
59: "discrete-time Lyapunov"
60: };
64: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
68: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
69: if (isascii) {
70: PetscViewerASCIIGetTab(viewer,&tabs);
71: PetscViewerASCIISetTab(viewer,((PetscObject)lme)->tablevel);
72: PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer);
73: if (lme->ops->view) {
74: PetscViewerASCIIPushTab(viewer);
75: (*lme->ops->view)(lme,viewer);
76: PetscViewerASCIIPopTab(viewer);
77: }
78: PetscViewerASCIIPrintf(viewer," equation type: %s\n",eqname[lme->problem_type]);
79: PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %D\n",lme->ncv);
80: PetscViewerASCIIPrintf(viewer," maximum number of iterations: %D\n",lme->max_it);
81: PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol);
82: PetscViewerASCIISetTab(viewer,tabs);
83: } else {
84: if (lme->ops->view) {
85: (*lme->ops->view)(lme,viewer);
86: }
87: }
88: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
89: if (!lme->V) { LMEGetBV(lme,&lme->V); }
90: BVView(lme->V,viewer);
91: PetscViewerPopFormat(viewer);
92: return(0);
93: }
95: /*@C
96: LMEReasonView - Displays the reason an LME solve converged or diverged.
98: Collective on LME
100: Parameter:
101: + lme - the linear matrix equation context
102: - viewer - the viewer to display the reason
104: Options Database Keys:
105: . -lme_converged_reason - print reason for convergence, and number of iterations
107: Level: intermediate
109: .seealso: LMESetTolerances(), LMEGetIterationNumber()
110: @*/
111: PetscErrorCode LMEReasonView(LME lme,PetscViewer viewer)
112: {
114: PetscBool isAscii;
117: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
118: if (isAscii) {
119: PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
120: if (lme->reason > 0) {
121: 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);
122: } else {
123: 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);
124: }
125: PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
126: }
127: return(0);
128: }
130: /*@
131: LMEReasonViewFromOptions - Processes command line options to determine if/how
132: the LME converged reason is to be viewed.
134: Collective on LME
136: Input Parameters:
137: . lme - the linear matrix equation context
139: Level: developer
140: @*/
141: PetscErrorCode LMEReasonViewFromOptions(LME lme)
142: {
143: PetscErrorCode ierr;
144: PetscViewer viewer;
145: PetscBool flg;
146: static PetscBool incall = PETSC_FALSE;
147: PetscViewerFormat format;
150: if (incall) return(0);
151: incall = PETSC_TRUE;
152: PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
153: if (flg) {
154: PetscViewerPushFormat(viewer,format);
155: LMEReasonView(lme,viewer);
156: PetscViewerPopFormat(viewer);
157: PetscViewerDestroy(&viewer);
158: }
159: incall = PETSC_FALSE;
160: return(0);
161: }
163: /*@
164: LMECreate - Creates the default LME context.
166: Collective on MPI_Comm
168: Input Parameter:
169: . comm - MPI communicator
171: Output Parameter:
172: . lme - location to put the LME context
174: Note:
175: The default LME type is LMEKRYLOV
177: Level: beginner
179: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
180: @*/
181: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
182: {
184: LME lme;
188: *outlme = 0;
189: LMEInitializePackage();
190: SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);
192: lme->A = NULL;
193: lme->B = NULL;
194: lme->D = NULL;
195: lme->E = NULL;
196: lme->C = NULL;
197: lme->X = NULL;
198: lme->problem_type = LME_LYAPUNOV;
199: lme->max_it = 0;
200: lme->ncv = 0;
201: lme->tol = PETSC_DEFAULT;
202: lme->errorifnotconverged = PETSC_FALSE;
204: lme->numbermonitors = 0;
206: lme->V = NULL;
207: lme->nwork = 0;
208: lme->work = NULL;
209: lme->data = NULL;
211: lme->its = 0;
212: lme->errest = 0;
213: lme->setupcalled = 0;
214: lme->reason = LME_CONVERGED_ITERATING;
216: *outlme = lme;
217: return(0);
218: }
220: /*@C
221: LMESetType - Selects the particular solver to be used in the LME object.
223: Logically Collective on LME
225: Input Parameters:
226: + lme - the linear matrix equation context
227: - type - a known method
229: Options Database Key:
230: . -lme_type <method> - Sets the method; use -help for a list
231: of available methods
233: Notes:
234: See "slepc/include/slepclme.h" for available methods. The default
235: is LMEKRYLOV
237: Normally, it is best to use the LMESetFromOptions() command and
238: then set the LME type from the options database rather than by using
239: this routine. Using the options database provides the user with
240: maximum flexibility in evaluating the different available methods.
241: The LMESetType() routine is provided for those situations where it
242: is necessary to set the iterative solver independently of the command
243: line or options database.
245: Level: intermediate
247: .seealso: LMEType
248: @*/
249: PetscErrorCode LMESetType(LME lme,LMEType type)
250: {
251: PetscErrorCode ierr,(*r)(LME);
252: PetscBool match;
258: PetscObjectTypeCompare((PetscObject)lme,type,&match);
259: if (match) return(0);
261: PetscFunctionListFind(LMEList,type,&r);
262: if (!r) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
264: if (lme->ops->destroy) { (*lme->ops->destroy)(lme); }
265: PetscMemzero(lme->ops,sizeof(struct _LMEOps));
267: lme->setupcalled = 0;
268: PetscObjectChangeTypeName((PetscObject)lme,type);
269: (*r)(lme);
270: return(0);
271: }
273: /*@C
274: LMEGetType - Gets the LME type as a string from the LME object.
276: Not Collective
278: Input Parameter:
279: . lme - the linear matrix equation context
281: Output Parameter:
282: . name - name of LME method
284: Level: intermediate
286: .seealso: LMESetType()
287: @*/
288: PetscErrorCode LMEGetType(LME lme,LMEType *type)
289: {
293: *type = ((PetscObject)lme)->type_name;
294: return(0);
295: }
297: /*@C
298: LMERegister - Adds a method to the linear matrix equation solver package.
300: Not Collective
302: Input Parameters:
303: + name - name of a new user-defined solver
304: - function - routine to create the solver context
306: Notes:
307: LMERegister() may be called multiple times to add several user-defined solvers.
309: Sample usage:
310: .vb
311: LMERegister("my_solver",MySolverCreate);
312: .ve
314: Then, your solver can be chosen with the procedural interface via
315: $ LMESetType(lme,"my_solver")
316: or at runtime via the option
317: $ -lme_type my_solver
319: Level: advanced
321: .seealso: LMERegisterAll()
322: @*/
323: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
324: {
328: PetscFunctionListAdd(&LMEList,name,function);
329: return(0);
330: }
332: /*@
333: LMEReset - Resets the LME context to the initial state (prior to setup)
334: and destroys any allocated Vecs and Mats.
336: Collective on LME
338: Input Parameter:
339: . lme - linear matrix equation context obtained from LMECreate()
341: Level: advanced
343: .seealso: LMEDestroy()
344: @*/
345: PetscErrorCode LMEReset(LME lme)
346: {
351: if (!lme) return(0);
352: if (lme->ops->reset) { (lme->ops->reset)(lme); }
353: MatDestroy(&lme->A);
354: MatDestroy(&lme->B);
355: MatDestroy(&lme->D);
356: MatDestroy(&lme->E);
357: MatDestroy(&lme->C);
358: MatDestroy(&lme->X);
359: BVDestroy(&lme->V);
360: VecDestroyVecs(lme->nwork,&lme->work);
361: lme->nwork = 0;
362: lme->setupcalled = 0;
363: return(0);
364: }
366: /*@
367: LMEDestroy - Destroys the LME context.
369: Collective on LME
371: Input Parameter:
372: . lme - linear matrix equation context obtained from LMECreate()
374: Level: beginner
376: .seealso: LMECreate(), LMESetUp(), LMESolve()
377: @*/
378: PetscErrorCode LMEDestroy(LME *lme)
379: {
383: if (!*lme) return(0);
385: if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
386: LMEReset(*lme);
387: if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
388: LMEMonitorCancel(*lme);
389: PetscHeaderDestroy(lme);
390: return(0);
391: }
393: /*@
394: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
396: Collective on LME
398: Input Parameters:
399: + lme - linear matrix equation context obtained from LMECreate()
400: - bv - the basis vectors object
402: Note:
403: Use LMEGetBV() to retrieve the basis vectors context (for example,
404: to free it at the end of the computations).
406: Level: advanced
408: .seealso: LMEGetBV()
409: @*/
410: PetscErrorCode LMESetBV(LME lme,BV bv)
411: {
418: PetscObjectReference((PetscObject)bv);
419: BVDestroy(&lme->V);
420: lme->V = bv;
421: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
422: return(0);
423: }
425: /*@
426: LMEGetBV - Obtain the basis vectors object associated to the matrix
427: function solver.
429: Not Collective
431: Input Parameters:
432: . lme - linear matrix equation context obtained from LMECreate()
434: Output Parameter:
435: . bv - basis vectors context
437: Level: advanced
439: .seealso: LMESetBV()
440: @*/
441: PetscErrorCode LMEGetBV(LME lme,BV *bv)
442: {
448: if (!lme->V) {
449: BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
450: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
451: }
452: *bv = lme->V;
453: return(0);
454: }