Actual source code: lmebasic.c
slepc-3.10.1 2018-10-23
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: LMEInitializePackage();
329: PetscFunctionListAdd(&LMEList,name,function);
330: return(0);
331: }
333: /*@
334: LMEReset - Resets the LME context to the initial state (prior to setup)
335: and destroys any allocated Vecs and Mats.
337: Collective on LME
339: Input Parameter:
340: . lme - linear matrix equation context obtained from LMECreate()
342: Level: advanced
344: .seealso: LMEDestroy()
345: @*/
346: PetscErrorCode LMEReset(LME lme)
347: {
352: if (!lme) return(0);
353: if (lme->ops->reset) { (lme->ops->reset)(lme); }
354: MatDestroy(&lme->A);
355: MatDestroy(&lme->B);
356: MatDestroy(&lme->D);
357: MatDestroy(&lme->E);
358: MatDestroy(&lme->C);
359: MatDestroy(&lme->X);
360: BVDestroy(&lme->V);
361: VecDestroyVecs(lme->nwork,&lme->work);
362: lme->nwork = 0;
363: lme->setupcalled = 0;
364: return(0);
365: }
367: /*@
368: LMEDestroy - Destroys the LME context.
370: Collective on LME
372: Input Parameter:
373: . lme - linear matrix equation context obtained from LMECreate()
375: Level: beginner
377: .seealso: LMECreate(), LMESetUp(), LMESolve()
378: @*/
379: PetscErrorCode LMEDestroy(LME *lme)
380: {
384: if (!*lme) return(0);
386: if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
387: LMEReset(*lme);
388: if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
389: LMEMonitorCancel(*lme);
390: PetscHeaderDestroy(lme);
391: return(0);
392: }
394: /*@
395: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
397: Collective on LME
399: Input Parameters:
400: + lme - linear matrix equation context obtained from LMECreate()
401: - bv - the basis vectors object
403: Note:
404: Use LMEGetBV() to retrieve the basis vectors context (for example,
405: to free it at the end of the computations).
407: Level: advanced
409: .seealso: LMEGetBV()
410: @*/
411: PetscErrorCode LMESetBV(LME lme,BV bv)
412: {
419: PetscObjectReference((PetscObject)bv);
420: BVDestroy(&lme->V);
421: lme->V = bv;
422: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
423: return(0);
424: }
426: /*@
427: LMEGetBV - Obtain the basis vectors object associated to the matrix
428: function solver.
430: Not Collective
432: Input Parameters:
433: . lme - linear matrix equation context obtained from LMECreate()
435: Output Parameter:
436: . bv - basis vectors context
438: Level: advanced
440: .seealso: LMESetBV()
441: @*/
442: PetscErrorCode LMEGetBV(LME lme,BV *bv)
443: {
449: if (!lme->V) {
450: BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
451: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
452: }
453: *bv = lme->V;
454: return(0);
455: }