Actual source code: lmebasic.c
slepc-3.13.4 2020-09-02
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: LMEReasonView - 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: Level: intermediate
128: .seealso: LMESetTolerances(), LMEGetIterationNumber()
129: @*/
130: PetscErrorCode LMEReasonView(LME lme,PetscViewer viewer)
131: {
133: PetscBool isAscii;
136: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
137: if (isAscii) {
138: PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
139: if (lme->reason > 0) {
140: 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);
141: } else {
142: 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);
143: }
144: PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
145: }
146: return(0);
147: }
149: /*@
150: LMEReasonViewFromOptions - Processes command line options to determine if/how
151: the LME converged reason is to be viewed.
153: Collective on lme
155: Input Parameter:
156: . lme - the linear matrix equation context
158: Level: developer
159: @*/
160: PetscErrorCode LMEReasonViewFromOptions(LME lme)
161: {
162: PetscErrorCode ierr;
163: PetscViewer viewer;
164: PetscBool flg;
165: static PetscBool incall = PETSC_FALSE;
166: PetscViewerFormat format;
169: if (incall) return(0);
170: incall = PETSC_TRUE;
171: PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
172: if (flg) {
173: PetscViewerPushFormat(viewer,format);
174: LMEReasonView(lme,viewer);
175: PetscViewerPopFormat(viewer);
176: PetscViewerDestroy(&viewer);
177: }
178: incall = PETSC_FALSE;
179: return(0);
180: }
182: /*@
183: LMECreate - Creates the default LME context.
185: Collective
187: Input Parameter:
188: . comm - MPI communicator
190: Output Parameter:
191: . lme - location to put the LME context
193: Note:
194: The default LME type is LMEKRYLOV
196: Level: beginner
198: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
199: @*/
200: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
201: {
203: LME lme;
207: *outlme = 0;
208: LMEInitializePackage();
209: SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);
211: lme->A = NULL;
212: lme->B = NULL;
213: lme->D = NULL;
214: lme->E = NULL;
215: lme->C = NULL;
216: lme->X = NULL;
217: lme->problem_type = LME_LYAPUNOV;
218: lme->max_it = PETSC_DEFAULT;
219: lme->ncv = PETSC_DEFAULT;
220: lme->tol = PETSC_DEFAULT;
221: lme->errorifnotconverged = PETSC_FALSE;
223: lme->numbermonitors = 0;
225: lme->V = NULL;
226: lme->nwork = 0;
227: lme->work = NULL;
228: lme->data = NULL;
230: lme->its = 0;
231: lme->errest = 0;
232: lme->setupcalled = 0;
233: lme->reason = LME_CONVERGED_ITERATING;
235: *outlme = lme;
236: return(0);
237: }
239: /*@C
240: LMESetType - Selects the particular solver to be used in the LME object.
242: Logically Collective on lme
244: Input Parameters:
245: + lme - the linear matrix equation context
246: - type - a known method
248: Options Database Key:
249: . -lme_type <method> - Sets the method; use -help for a list
250: of available methods
252: Notes:
253: See "slepc/include/slepclme.h" for available methods. The default
254: is LMEKRYLOV
256: Normally, it is best to use the LMESetFromOptions() command and
257: then set the LME type from the options database rather than by using
258: this routine. Using the options database provides the user with
259: maximum flexibility in evaluating the different available methods.
260: The LMESetType() routine is provided for those situations where it
261: is necessary to set the iterative solver independently of the command
262: line or options database.
264: Level: intermediate
266: .seealso: LMEType
267: @*/
268: PetscErrorCode LMESetType(LME lme,LMEType type)
269: {
270: PetscErrorCode ierr,(*r)(LME);
271: PetscBool match;
277: PetscObjectTypeCompare((PetscObject)lme,type,&match);
278: if (match) return(0);
280: PetscFunctionListFind(LMEList,type,&r);
281: if (!r) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
283: if (lme->ops->destroy) { (*lme->ops->destroy)(lme); }
284: PetscMemzero(lme->ops,sizeof(struct _LMEOps));
286: lme->setupcalled = 0;
287: PetscObjectChangeTypeName((PetscObject)lme,type);
288: (*r)(lme);
289: return(0);
290: }
292: /*@C
293: LMEGetType - Gets the LME type as a string from the LME object.
295: Not Collective
297: Input Parameter:
298: . lme - the linear matrix equation context
300: Output Parameter:
301: . name - name of LME method
303: Level: intermediate
305: .seealso: LMESetType()
306: @*/
307: PetscErrorCode LMEGetType(LME lme,LMEType *type)
308: {
312: *type = ((PetscObject)lme)->type_name;
313: return(0);
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: Notes:
326: LMERegister() may be called multiple times to add several user-defined solvers.
328: Sample usage:
329: .vb
330: LMERegister("my_solver",MySolverCreate);
331: .ve
333: Then, your solver can be chosen with the procedural interface via
334: $ LMESetType(lme,"my_solver")
335: or at runtime via the option
336: $ -lme_type my_solver
338: Level: advanced
340: .seealso: LMERegisterAll()
341: @*/
342: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
343: {
347: LMEInitializePackage();
348: PetscFunctionListAdd(&LMEList,name,function);
349: return(0);
350: }
352: /*@
353: LMEReset - Resets the LME context to the initial state (prior to setup)
354: and destroys any allocated Vecs and Mats.
356: Collective on lme
358: Input Parameter:
359: . lme - linear matrix equation context obtained from LMECreate()
361: Level: advanced
363: .seealso: LMEDestroy()
364: @*/
365: PetscErrorCode LMEReset(LME lme)
366: {
371: if (!lme) return(0);
372: if (lme->ops->reset) { (lme->ops->reset)(lme); }
373: MatDestroy(&lme->A);
374: MatDestroy(&lme->B);
375: MatDestroy(&lme->D);
376: MatDestroy(&lme->E);
377: MatDestroy(&lme->C);
378: MatDestroy(&lme->X);
379: BVDestroy(&lme->V);
380: VecDestroyVecs(lme->nwork,&lme->work);
381: lme->nwork = 0;
382: lme->setupcalled = 0;
383: return(0);
384: }
386: /*@
387: LMEDestroy - Destroys the LME context.
389: Collective on lme
391: Input Parameter:
392: . lme - linear matrix equation context obtained from LMECreate()
394: Level: beginner
396: .seealso: LMECreate(), LMESetUp(), LMESolve()
397: @*/
398: PetscErrorCode LMEDestroy(LME *lme)
399: {
403: if (!*lme) return(0);
405: if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
406: LMEReset(*lme);
407: if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
408: LMEMonitorCancel(*lme);
409: PetscHeaderDestroy(lme);
410: return(0);
411: }
413: /*@
414: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
416: Collective on lme
418: Input Parameters:
419: + lme - linear matrix equation context obtained from LMECreate()
420: - bv - the basis vectors object
422: Note:
423: Use LMEGetBV() to retrieve the basis vectors context (for example,
424: to free it at the end of the computations).
426: Level: advanced
428: .seealso: LMEGetBV()
429: @*/
430: PetscErrorCode LMESetBV(LME lme,BV bv)
431: {
438: PetscObjectReference((PetscObject)bv);
439: BVDestroy(&lme->V);
440: lme->V = bv;
441: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
442: return(0);
443: }
445: /*@
446: LMEGetBV - Obtain the basis vectors object associated to the matrix
447: function solver.
449: Not Collective
451: Input Parameters:
452: . lme - linear matrix equation context obtained from LMECreate()
454: Output Parameter:
455: . bv - basis vectors context
457: Level: advanced
459: .seealso: LMESetBV()
460: @*/
461: PetscErrorCode LMEGetBV(LME lme,BV *bv)
462: {
468: if (!lme->V) {
469: BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
470: PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
471: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
472: PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
473: }
474: *bv = lme->V;
475: return(0);
476: }