Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/lmeimpl.h> /*I "slepclme.h" I*/
15 :
16 : /* Logging support */
17 : PetscClassId LME_CLASSID = 0;
18 : PetscLogEvent LME_SetUp = 0,LME_Solve = 0,LME_ComputeError = 0;
19 :
20 : /* List of registered LME routines */
21 : PetscFunctionList LMEList = NULL;
22 : PetscBool LMERegisterAllCalled = PETSC_FALSE;
23 :
24 : /* List of registered LME monitors */
25 : PetscFunctionList LMEMonitorList = NULL;
26 : PetscFunctionList LMEMonitorCreateList = NULL;
27 : PetscFunctionList LMEMonitorDestroyList = NULL;
28 : PetscBool LMEMonitorRegisterAllCalled = PETSC_FALSE;
29 :
30 : /*@C
31 : LMEView - Prints the LME data structure.
32 :
33 : Collective
34 :
35 : Input Parameters:
36 : + lme - the linear matrix equation solver context
37 : - viewer - optional visualization context
38 :
39 : Options Database Key:
40 : . -lme_view - Calls LMEView() at end of LMESolve()
41 :
42 : Note:
43 : The available visualization contexts include
44 : + PETSC_VIEWER_STDOUT_SELF - standard output (default)
45 : - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
46 : output where only the first processor opens
47 : the file. All other processors send their
48 : data to the first processor to print.
49 :
50 : The user can open an alternative visualization context with
51 : PetscViewerASCIIOpen() - output to a specified file.
52 :
53 : Level: beginner
54 :
55 : .seealso: LMECreate()
56 : @*/
57 2 : PetscErrorCode LMEView(LME lme,PetscViewer viewer)
58 : {
59 2 : PetscBool isascii;
60 2 : const char *eqname[] = {
61 : "continuous-time Lyapunov",
62 : "continuous-time Sylvester",
63 : "generalized Lyapunov",
64 : "generalized Sylvester",
65 : "Stein",
66 : "discrete-time Lyapunov"
67 : };
68 :
69 2 : PetscFunctionBegin;
70 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
71 2 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer));
72 2 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
73 2 : PetscCheckSameComm(lme,1,viewer,2);
74 :
75 2 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
76 2 : if (isascii) {
77 2 : PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer));
78 2 : PetscCall(PetscViewerASCIIPushTab(viewer));
79 2 : PetscTryTypeMethod(lme,view,viewer);
80 2 : PetscCall(PetscViewerASCIIPopTab(viewer));
81 2 : PetscCall(PetscViewerASCIIPrintf(viewer," equation type: %s\n",eqname[lme->problem_type]));
82 2 : PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",lme->ncv));
83 2 : PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",lme->max_it));
84 2 : PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)lme->tol));
85 0 : } else PetscTryTypeMethod(lme,view,viewer);
86 2 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
87 2 : if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
88 2 : PetscCall(BVView(lme->V,viewer));
89 2 : PetscCall(PetscViewerPopFormat(viewer));
90 2 : PetscFunctionReturn(PETSC_SUCCESS);
91 : }
92 :
93 : /*@C
94 : LMEViewFromOptions - View from options
95 :
96 : Collective
97 :
98 : Input Parameters:
99 : + lme - the linear matrix equation context
100 : . obj - optional object
101 : - name - command line option
102 :
103 : Level: intermediate
104 :
105 : .seealso: LMEView(), LMECreate()
106 : @*/
107 66 : PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
108 : {
109 66 : PetscFunctionBegin;
110 66 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
111 66 : PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
112 66 : PetscFunctionReturn(PETSC_SUCCESS);
113 : }
114 : /*@C
115 : LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.
116 :
117 : Collective
118 :
119 : Input Parameters:
120 : + lme - the linear matrix equation context
121 : - viewer - the viewer to display the reason
122 :
123 : Options Database Keys:
124 : . -lme_converged_reason - print reason for convergence, and number of iterations
125 :
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
131 :
132 : Level: intermediate
133 :
134 : .seealso: LMESetTolerances(), LMEGetIterationNumber(), LMEConvergedReasonViewFromOptions()
135 : @*/
136 1 : PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
137 : {
138 1 : PetscBool isAscii;
139 1 : PetscViewerFormat format;
140 :
141 1 : PetscFunctionBegin;
142 1 : if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
143 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
144 1 : if (isAscii) {
145 1 : PetscCall(PetscViewerGetFormat(viewer,&format));
146 1 : PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel));
147 1 : if (lme->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve converged due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its));
148 0 : else if (lme->reason <= 0) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Linear matrix equation solve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)lme)->prefix?((PetscObject)lme)->prefix:"",LMEConvergedReasons[lme->reason],lme->its));
149 1 : PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel));
150 : }
151 1 : PetscFunctionReturn(PETSC_SUCCESS);
152 : }
153 :
154 : /*@
155 : LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
156 : the LME converged reason is to be viewed.
157 :
158 : Collective
159 :
160 : Input Parameter:
161 : . lme - the linear matrix equation context
162 :
163 : Level: developer
164 :
165 : .seealso: LMEConvergedReasonView()
166 : @*/
167 33 : PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
168 : {
169 33 : PetscViewer viewer;
170 33 : PetscBool flg;
171 33 : static PetscBool incall = PETSC_FALSE;
172 33 : PetscViewerFormat format;
173 :
174 33 : PetscFunctionBegin;
175 33 : if (incall) PetscFunctionReturn(PETSC_SUCCESS);
176 33 : incall = PETSC_TRUE;
177 33 : PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
178 33 : if (flg) {
179 1 : PetscCall(PetscViewerPushFormat(viewer,format));
180 1 : PetscCall(LMEConvergedReasonView(lme,viewer));
181 1 : PetscCall(PetscViewerPopFormat(viewer));
182 1 : PetscCall(PetscOptionsRestoreViewer(&viewer));
183 : }
184 33 : incall = PETSC_FALSE;
185 33 : PetscFunctionReturn(PETSC_SUCCESS);
186 : }
187 :
188 : /*@
189 : LMECreate - Creates the default LME context.
190 :
191 : Collective
192 :
193 : Input Parameter:
194 : . comm - MPI communicator
195 :
196 : Output Parameter:
197 : . outlme - location to put the LME context
198 :
199 : Note:
200 : The default LME type is LMEKRYLOV
201 :
202 : Level: beginner
203 :
204 : .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
205 : @*/
206 9 : PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
207 : {
208 9 : LME lme;
209 :
210 9 : PetscFunctionBegin;
211 9 : PetscAssertPointer(outlme,2);
212 9 : *outlme = NULL;
213 9 : PetscCall(LMEInitializePackage());
214 9 : PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));
215 :
216 9 : lme->A = NULL;
217 9 : lme->B = NULL;
218 9 : lme->D = NULL;
219 9 : lme->E = NULL;
220 9 : lme->C = NULL;
221 9 : lme->X = NULL;
222 9 : lme->problem_type = LME_LYAPUNOV;
223 9 : lme->max_it = PETSC_DEFAULT;
224 9 : lme->ncv = PETSC_DEFAULT;
225 9 : lme->tol = PETSC_DEFAULT;
226 9 : lme->errorifnotconverged = PETSC_FALSE;
227 :
228 9 : lme->numbermonitors = 0;
229 :
230 9 : lme->V = NULL;
231 9 : lme->nwork = 0;
232 9 : lme->work = NULL;
233 9 : lme->data = NULL;
234 :
235 9 : lme->its = 0;
236 9 : lme->errest = 0;
237 9 : lme->setupcalled = 0;
238 9 : lme->reason = LME_CONVERGED_ITERATING;
239 :
240 9 : *outlme = lme;
241 9 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 : /*@C
245 : LMESetType - Selects the particular solver to be used in the LME object.
246 :
247 : Logically Collective
248 :
249 : Input Parameters:
250 : + lme - the linear matrix equation context
251 : - type - a known method
252 :
253 : Options Database Key:
254 : . -lme_type <method> - Sets the method; use -help for a list
255 : of available methods
256 :
257 : Notes:
258 : See "slepc/include/slepclme.h" for available methods. The default
259 : is LMEKRYLOV
260 :
261 : Normally, it is best to use the LMESetFromOptions() command and
262 : then set the LME type from the options database rather than by using
263 : this routine. Using the options database provides the user with
264 : maximum flexibility in evaluating the different available methods.
265 : The LMESetType() routine is provided for those situations where it
266 : is necessary to set the iterative solver independently of the command
267 : line or options database.
268 :
269 : Level: intermediate
270 :
271 : .seealso: LMEType
272 : @*/
273 8 : PetscErrorCode LMESetType(LME lme,LMEType type)
274 : {
275 8 : PetscErrorCode (*r)(LME);
276 8 : PetscBool match;
277 :
278 8 : PetscFunctionBegin;
279 8 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
280 8 : PetscAssertPointer(type,2);
281 :
282 8 : PetscCall(PetscObjectTypeCompare((PetscObject)lme,type,&match));
283 8 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
284 :
285 8 : PetscCall(PetscFunctionListFind(LMEList,type,&r));
286 8 : PetscCheck(r,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
287 :
288 8 : PetscTryTypeMethod(lme,destroy);
289 8 : PetscCall(PetscMemzero(lme->ops,sizeof(struct _LMEOps)));
290 :
291 8 : lme->setupcalled = 0;
292 8 : PetscCall(PetscObjectChangeTypeName((PetscObject)lme,type));
293 8 : PetscCall((*r)(lme));
294 8 : PetscFunctionReturn(PETSC_SUCCESS);
295 : }
296 :
297 : /*@C
298 : LMEGetType - Gets the LME type as a string from the LME object.
299 :
300 : Not Collective
301 :
302 : Input Parameter:
303 : . lme - the linear matrix equation context
304 :
305 : Output Parameter:
306 : . type - name of LME method
307 :
308 : Level: intermediate
309 :
310 : .seealso: LMESetType()
311 : @*/
312 3 : PetscErrorCode LMEGetType(LME lme,LMEType *type)
313 : {
314 3 : PetscFunctionBegin;
315 3 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
316 3 : PetscAssertPointer(type,2);
317 3 : *type = ((PetscObject)lme)->type_name;
318 3 : PetscFunctionReturn(PETSC_SUCCESS);
319 : }
320 :
321 : /*@C
322 : LMERegister - Adds a method to the linear matrix equation solver package.
323 :
324 : Not Collective
325 :
326 : Input Parameters:
327 : + name - name of a new user-defined solver
328 : - function - routine to create the solver context
329 :
330 : Notes:
331 : LMERegister() may be called multiple times to add several user-defined solvers.
332 :
333 : Example Usage:
334 : .vb
335 : LMERegister("my_solver",MySolverCreate);
336 : .ve
337 :
338 : Then, your solver can be chosen with the procedural interface via
339 : $ LMESetType(lme,"my_solver")
340 : or at runtime via the option
341 : $ -lme_type my_solver
342 :
343 : Level: advanced
344 :
345 : .seealso: LMERegisterAll()
346 : @*/
347 10 : PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
348 : {
349 10 : PetscFunctionBegin;
350 10 : PetscCall(LMEInitializePackage());
351 10 : PetscCall(PetscFunctionListAdd(&LMEList,name,function));
352 10 : PetscFunctionReturn(PETSC_SUCCESS);
353 : }
354 :
355 : /*@C
356 : LMEMonitorRegister - Adds LME monitor routine.
357 :
358 : Not Collective
359 :
360 : Input Parameters:
361 : + name - name of a new monitor routine
362 : . vtype - a PetscViewerType for the output
363 : . format - a PetscViewerFormat for the output
364 : . monitor - monitor routine
365 : . create - creation routine, or NULL
366 : - destroy - destruction routine, or NULL
367 :
368 : Notes:
369 : LMEMonitorRegister() may be called multiple times to add several user-defined monitors.
370 :
371 : Example Usage:
372 : .vb
373 : LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
374 : .ve
375 :
376 : Then, your monitor can be chosen with the procedural interface via
377 : $ LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
378 : or at runtime via the option
379 : $ -lme_monitor_my_monitor
380 :
381 : Level: advanced
382 :
383 : .seealso: LMEMonitorRegisterAll()
384 : @*/
385 20 : PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(LME,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
386 : {
387 20 : char key[PETSC_MAX_PATH_LEN];
388 :
389 20 : PetscFunctionBegin;
390 20 : PetscCall(LMEInitializePackage());
391 20 : PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
392 20 : PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
393 20 : if (create) PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
394 20 : if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
395 20 : PetscFunctionReturn(PETSC_SUCCESS);
396 : }
397 :
398 : /*@
399 : LMEReset - Resets the LME context to the initial state (prior to setup)
400 : and destroys any allocated Vecs and Mats.
401 :
402 : Collective
403 :
404 : Input Parameter:
405 : . lme - linear matrix equation context obtained from LMECreate()
406 :
407 : Level: advanced
408 :
409 : .seealso: LMEDestroy()
410 : @*/
411 9 : PetscErrorCode LMEReset(LME lme)
412 : {
413 9 : PetscFunctionBegin;
414 9 : if (lme) PetscValidHeaderSpecific(lme,LME_CLASSID,1);
415 0 : if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
416 9 : PetscTryTypeMethod(lme,reset);
417 9 : PetscCall(MatDestroy(&lme->A));
418 9 : PetscCall(MatDestroy(&lme->B));
419 9 : PetscCall(MatDestroy(&lme->D));
420 9 : PetscCall(MatDestroy(&lme->E));
421 9 : PetscCall(MatDestroy(&lme->C));
422 9 : PetscCall(MatDestroy(&lme->X));
423 9 : PetscCall(BVDestroy(&lme->V));
424 9 : PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
425 9 : lme->nwork = 0;
426 9 : lme->setupcalled = 0;
427 9 : PetscFunctionReturn(PETSC_SUCCESS);
428 : }
429 :
430 : /*@C
431 : LMEDestroy - Destroys the LME context.
432 :
433 : Collective
434 :
435 : Input Parameter:
436 : . lme - linear matrix equation context obtained from LMECreate()
437 :
438 : Level: beginner
439 :
440 : .seealso: LMECreate(), LMESetUp(), LMESolve()
441 : @*/
442 9 : PetscErrorCode LMEDestroy(LME *lme)
443 : {
444 9 : PetscFunctionBegin;
445 9 : if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
446 9 : PetscValidHeaderSpecific(*lme,LME_CLASSID,1);
447 9 : if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
448 9 : PetscCall(LMEReset(*lme));
449 9 : PetscTryTypeMethod(*lme,destroy);
450 9 : PetscCall(LMEMonitorCancel(*lme));
451 9 : PetscCall(PetscHeaderDestroy(lme));
452 9 : PetscFunctionReturn(PETSC_SUCCESS);
453 : }
454 :
455 : /*@
456 : LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
457 :
458 : Collective
459 :
460 : Input Parameters:
461 : + lme - linear matrix equation context obtained from LMECreate()
462 : - bv - the basis vectors object
463 :
464 : Note:
465 : Use LMEGetBV() to retrieve the basis vectors context (for example,
466 : to free it at the end of the computations).
467 :
468 : Level: advanced
469 :
470 : .seealso: LMEGetBV()
471 : @*/
472 0 : PetscErrorCode LMESetBV(LME lme,BV bv)
473 : {
474 0 : PetscFunctionBegin;
475 0 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
476 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,2);
477 0 : PetscCheckSameComm(lme,1,bv,2);
478 0 : PetscCall(PetscObjectReference((PetscObject)bv));
479 0 : PetscCall(BVDestroy(&lme->V));
480 0 : lme->V = bv;
481 0 : PetscFunctionReturn(PETSC_SUCCESS);
482 : }
483 :
484 : /*@
485 : LMEGetBV - Obtain the basis vectors object associated to the matrix
486 : function solver.
487 :
488 : Not Collective
489 :
490 : Input Parameters:
491 : . lme - linear matrix equation context obtained from LMECreate()
492 :
493 : Output Parameter:
494 : . bv - basis vectors context
495 :
496 : Level: advanced
497 :
498 : .seealso: LMESetBV()
499 : @*/
500 8 : PetscErrorCode LMEGetBV(LME lme,BV *bv)
501 : {
502 8 : PetscFunctionBegin;
503 8 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
504 8 : PetscAssertPointer(bv,2);
505 8 : if (!lme->V) {
506 8 : PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
507 8 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
508 8 : PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
509 : }
510 8 : *bv = lme->V;
511 8 : PetscFunctionReturn(PETSC_SUCCESS);
512 : }
|