Actual source code: lmebasic.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: const char *eqname[] = {
53: "continuous-time Lyapunov",
54: "continuous-time Sylvester",
55: "generalized Lyapunov",
56: "generalized Sylvester",
57: "Stein",
58: "discrete-time Lyapunov"
59: };
63: if (!viewer) {
64: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer);
65: }
69: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
70: if (isascii) {
71: PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer);
72: if (lme->ops->view) {
73: PetscViewerASCIIPushTab(viewer);
74: (*lme->ops->view)(lme,viewer);
75: PetscViewerASCIIPopTab(viewer);
76: }
77: PetscViewerASCIIPrintf(viewer," equation type: %s\n",eqname[lme->problem_type]);
78: PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %D\n",lme->ncv);
79: PetscViewerASCIIPrintf(viewer," maximum number of iterations: %D\n",lme->max_it);
80: PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol);
81: } else {
82: if (lme->ops->view) {
83: (*lme->ops->view)(lme,viewer);
84: }
85: }
86: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
87: if (!lme->V) { LMEGetBV(lme,&lme->V); }
88: BVView(lme->V,viewer);
89: PetscViewerPopFormat(viewer);
90: return(0);
91: }
93: /*@C
94: LMEReasonView - Displays the reason an LME solve converged or diverged.
96: Collective on LME
98: Parameter:
99: + lme - the linear matrix equation context
100: - viewer - the viewer to display the reason
102: Options Database Keys:
103: . -lme_converged_reason - print reason for convergence, and number of iterations
105: Level: intermediate
107: .seealso: LMESetTolerances(), LMEGetIterationNumber()
108: @*/
109: PetscErrorCode LMEReasonView(LME lme,PetscViewer viewer)
110: {
112: PetscBool isAscii;
115: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
116: if (isAscii) {
117: PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
118: if (lme->reason > 0) {
119: 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);
120: } else {
121: 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);
122: }
123: PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
124: }
125: return(0);
126: }
128: /*@
129: LMEReasonViewFromOptions - Processes command line options to determine if/how
130: the LME converged reason is to be viewed.
132: Collective on LME
134: Input Parameters:
135: . lme - the linear matrix equation context
137: Level: developer
138: @*/
139: PetscErrorCode LMEReasonViewFromOptions(LME lme)
140: {
141: PetscErrorCode ierr;
142: PetscViewer viewer;
143: PetscBool flg;
144: static PetscBool incall = PETSC_FALSE;
145: PetscViewerFormat format;
148: if (incall) return(0);
149: incall = PETSC_TRUE;
150: PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
151: if (flg) {
152: PetscViewerPushFormat(viewer,format);
153: LMEReasonView(lme,viewer);
154: PetscViewerPopFormat(viewer);
155: PetscViewerDestroy(&viewer);
156: }
157: incall = PETSC_FALSE;
158: return(0);
159: }
161: /*@
162: LMECreate - Creates the default LME context.
164: Collective on MPI_Comm
166: Input Parameter:
167: . comm - MPI communicator
169: Output Parameter:
170: . lme - location to put the LME context
172: Note:
173: The default LME type is LMEKRYLOV
175: Level: beginner
177: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
178: @*/
179: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
180: {
182: LME lme;
186: *outlme = 0;
187: LMEInitializePackage();
188: SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);
190: lme->A = NULL;
191: lme->B = NULL;
192: lme->D = NULL;
193: lme->E = NULL;
194: lme->C = NULL;
195: lme->X = NULL;
196: lme->problem_type = LME_LYAPUNOV;
197: lme->max_it = 0;
198: lme->ncv = 0;
199: lme->tol = PETSC_DEFAULT;
200: lme->errorifnotconverged = PETSC_FALSE;
202: lme->numbermonitors = 0;
204: lme->V = NULL;
205: lme->nwork = 0;
206: lme->work = NULL;
207: lme->data = NULL;
209: lme->its = 0;
210: lme->errest = 0;
211: lme->setupcalled = 0;
212: lme->reason = LME_CONVERGED_ITERATING;
214: *outlme = lme;
215: return(0);
216: }
218: /*@C
219: LMESetType - Selects the particular solver to be used in the LME object.
221: Logically Collective on LME
223: Input Parameters:
224: + lme - the linear matrix equation context
225: - type - a known method
227: Options Database Key:
228: . -lme_type <method> - Sets the method; use -help for a list
229: of available methods
231: Notes:
232: See "slepc/include/slepclme.h" for available methods. The default
233: is LMEKRYLOV
235: Normally, it is best to use the LMESetFromOptions() command and
236: then set the LME type from the options database rather than by using
237: this routine. Using the options database provides the user with
238: maximum flexibility in evaluating the different available methods.
239: The LMESetType() routine is provided for those situations where it
240: is necessary to set the iterative solver independently of the command
241: line or options database.
243: Level: intermediate
245: .seealso: LMEType
246: @*/
247: PetscErrorCode LMESetType(LME lme,LMEType type)
248: {
249: PetscErrorCode ierr,(*r)(LME);
250: PetscBool match;
256: PetscObjectTypeCompare((PetscObject)lme,type,&match);
257: if (match) return(0);
259: PetscFunctionListFind(LMEList,type,&r);
260: if (!r) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
262: if (lme->ops->destroy) { (*lme->ops->destroy)(lme); }
263: PetscMemzero(lme->ops,sizeof(struct _LMEOps));
265: lme->setupcalled = 0;
266: PetscObjectChangeTypeName((PetscObject)lme,type);
267: (*r)(lme);
268: return(0);
269: }
271: /*@C
272: LMEGetType - Gets the LME type as a string from the LME object.
274: Not Collective
276: Input Parameter:
277: . lme - the linear matrix equation context
279: Output Parameter:
280: . name - name of LME method
282: Level: intermediate
284: .seealso: LMESetType()
285: @*/
286: PetscErrorCode LMEGetType(LME lme,LMEType *type)
287: {
291: *type = ((PetscObject)lme)->type_name;
292: return(0);
293: }
295: /*@C
296: LMERegister - Adds a method to the linear matrix equation solver package.
298: Not Collective
300: Input Parameters:
301: + name - name of a new user-defined solver
302: - function - routine to create the solver context
304: Notes:
305: LMERegister() may be called multiple times to add several user-defined solvers.
307: Sample usage:
308: .vb
309: LMERegister("my_solver",MySolverCreate);
310: .ve
312: Then, your solver can be chosen with the procedural interface via
313: $ LMESetType(lme,"my_solver")
314: or at runtime via the option
315: $ -lme_type my_solver
317: Level: advanced
319: .seealso: LMERegisterAll()
320: @*/
321: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
322: {
326: LMEInitializePackage();
327: PetscFunctionListAdd(&LMEList,name,function);
328: return(0);
329: }
331: /*@
332: LMEReset - Resets the LME context to the initial state (prior to setup)
333: and destroys any allocated Vecs and Mats.
335: Collective on LME
337: Input Parameter:
338: . lme - linear matrix equation context obtained from LMECreate()
340: Level: advanced
342: .seealso: LMEDestroy()
343: @*/
344: PetscErrorCode LMEReset(LME lme)
345: {
350: if (!lme) return(0);
351: if (lme->ops->reset) { (lme->ops->reset)(lme); }
352: MatDestroy(&lme->A);
353: MatDestroy(&lme->B);
354: MatDestroy(&lme->D);
355: MatDestroy(&lme->E);
356: MatDestroy(&lme->C);
357: MatDestroy(&lme->X);
358: BVDestroy(&lme->V);
359: VecDestroyVecs(lme->nwork,&lme->work);
360: lme->nwork = 0;
361: lme->setupcalled = 0;
362: return(0);
363: }
365: /*@
366: LMEDestroy - Destroys the LME context.
368: Collective on LME
370: Input Parameter:
371: . lme - linear matrix equation context obtained from LMECreate()
373: Level: beginner
375: .seealso: LMECreate(), LMESetUp(), LMESolve()
376: @*/
377: PetscErrorCode LMEDestroy(LME *lme)
378: {
382: if (!*lme) return(0);
384: if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
385: LMEReset(*lme);
386: if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
387: LMEMonitorCancel(*lme);
388: PetscHeaderDestroy(lme);
389: return(0);
390: }
392: /*@
393: LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
395: Collective on LME
397: Input Parameters:
398: + lme - linear matrix equation context obtained from LMECreate()
399: - bv - the basis vectors object
401: Note:
402: Use LMEGetBV() to retrieve the basis vectors context (for example,
403: to free it at the end of the computations).
405: Level: advanced
407: .seealso: LMEGetBV()
408: @*/
409: PetscErrorCode LMESetBV(LME lme,BV bv)
410: {
417: PetscObjectReference((PetscObject)bv);
418: BVDestroy(&lme->V);
419: lme->V = bv;
420: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
421: return(0);
422: }
424: /*@
425: LMEGetBV - Obtain the basis vectors object associated to the matrix
426: function solver.
428: Not Collective
430: Input Parameters:
431: . lme - linear matrix equation context obtained from LMECreate()
433: Output Parameter:
434: . bv - basis vectors context
436: Level: advanced
438: .seealso: LMESetBV()
439: @*/
440: PetscErrorCode LMEGetBV(LME lme,BV *bv)
441: {
447: if (!lme->V) {
448: BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
449: PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
450: PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
451: PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
452: }
453: *bv = lme->V;
454: return(0);
455: }