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 : /*@
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 : /*@
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 64 : PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
108 : {
109 64 : PetscFunctionBegin;
110 64 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
111 64 : PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
112 64 : PetscFunctionReturn(PETSC_SUCCESS);
113 : }
114 : /*@
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 32 : PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
168 : {
169 32 : PetscViewer viewer;
170 32 : PetscBool flg;
171 32 : static PetscBool incall = PETSC_FALSE;
172 32 : PetscViewerFormat format;
173 :
174 32 : PetscFunctionBegin;
175 32 : if (incall) PetscFunctionReturn(PETSC_SUCCESS);
176 32 : incall = PETSC_TRUE;
177 32 : PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
178 32 : if (flg) {
179 1 : PetscCall(PetscViewerPushFormat(viewer,format));
180 1 : PetscCall(LMEConvergedReasonView(lme,viewer));
181 1 : PetscCall(PetscViewerPopFormat(viewer));
182 1 : PetscCall(PetscViewerDestroy(&viewer));
183 : }
184 32 : incall = PETSC_FALSE;
185 32 : 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 8 : PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
207 : {
208 8 : LME lme;
209 :
210 8 : PetscFunctionBegin;
211 8 : PetscAssertPointer(outlme,2);
212 8 : PetscCall(LMEInitializePackage());
213 8 : PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));
214 :
215 8 : lme->A = NULL;
216 8 : lme->B = NULL;
217 8 : lme->D = NULL;
218 8 : lme->E = NULL;
219 8 : lme->C = NULL;
220 8 : lme->X = NULL;
221 8 : lme->problem_type = LME_LYAPUNOV;
222 8 : lme->max_it = PETSC_DETERMINE;
223 8 : lme->ncv = PETSC_DETERMINE;
224 8 : lme->tol = PETSC_DETERMINE;
225 8 : lme->errorifnotconverged = PETSC_FALSE;
226 :
227 8 : lme->numbermonitors = 0;
228 :
229 8 : lme->V = NULL;
230 8 : lme->nwork = 0;
231 8 : lme->work = NULL;
232 8 : lme->data = NULL;
233 :
234 8 : lme->its = 0;
235 8 : lme->errest = 0;
236 8 : lme->setupcalled = 0;
237 8 : lme->reason = LME_CONVERGED_ITERATING;
238 :
239 8 : *outlme = lme;
240 8 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 : /*@
244 : LMESetType - Selects the particular solver to be used in the LME object.
245 :
246 : Logically Collective
247 :
248 : Input Parameters:
249 : + lme - the linear matrix equation context
250 : - type - a known method
251 :
252 : Options Database Key:
253 : . -lme_type <method> - Sets the method; use -help for a list
254 : of available methods
255 :
256 : Notes:
257 : See "slepc/include/slepclme.h" for available methods. The default
258 : is LMEKRYLOV
259 :
260 : Normally, it is best to use the LMESetFromOptions() command and
261 : then set the LME type from the options database rather than by using
262 : this routine. Using the options database provides the user with
263 : maximum flexibility in evaluating the different available methods.
264 : The LMESetType() routine is provided for those situations where it
265 : is necessary to set the iterative solver independently of the command
266 : line or options database.
267 :
268 : Level: intermediate
269 :
270 : .seealso: LMEType
271 : @*/
272 7 : PetscErrorCode LMESetType(LME lme,LMEType type)
273 : {
274 7 : PetscErrorCode (*r)(LME);
275 7 : PetscBool match;
276 :
277 7 : PetscFunctionBegin;
278 7 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
279 7 : PetscAssertPointer(type,2);
280 :
281 7 : PetscCall(PetscObjectTypeCompare((PetscObject)lme,type,&match));
282 7 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
283 :
284 7 : PetscCall(PetscFunctionListFind(LMEList,type,&r));
285 7 : PetscCheck(r,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);
286 :
287 7 : PetscTryTypeMethod(lme,destroy);
288 7 : PetscCall(PetscMemzero(lme->ops,sizeof(struct _LMEOps)));
289 :
290 7 : lme->setupcalled = 0;
291 7 : PetscCall(PetscObjectChangeTypeName((PetscObject)lme,type));
292 7 : PetscCall((*r)(lme));
293 7 : PetscFunctionReturn(PETSC_SUCCESS);
294 : }
295 :
296 : /*@
297 : LMEGetType - Gets the LME type as a string from the LME object.
298 :
299 : Not Collective
300 :
301 : Input Parameter:
302 : . lme - the linear matrix equation context
303 :
304 : Output Parameter:
305 : . type - name of LME method
306 :
307 : Level: intermediate
308 :
309 : .seealso: LMESetType()
310 : @*/
311 2 : PetscErrorCode LMEGetType(LME lme,LMEType *type)
312 : {
313 2 : PetscFunctionBegin;
314 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
315 2 : PetscAssertPointer(type,2);
316 2 : *type = ((PetscObject)lme)->type_name;
317 2 : PetscFunctionReturn(PETSC_SUCCESS);
318 : }
319 :
320 : /*@C
321 : LMERegister - Adds a method to the linear matrix equation solver package.
322 :
323 : Not Collective
324 :
325 : Input Parameters:
326 : + name - name of a new user-defined solver
327 : - function - routine to create the solver context
328 :
329 : Notes:
330 : LMERegister() may be called multiple times to add several user-defined solvers.
331 :
332 : Example Usage:
333 : .vb
334 : LMERegister("my_solver",MySolverCreate);
335 : .ve
336 :
337 : Then, your solver can be chosen with the procedural interface via
338 : $ LMESetType(lme,"my_solver")
339 : or at runtime via the option
340 : $ -lme_type my_solver
341 :
342 : Level: advanced
343 :
344 : .seealso: LMERegisterAll()
345 : @*/
346 9 : PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
347 : {
348 9 : PetscFunctionBegin;
349 9 : PetscCall(LMEInitializePackage());
350 9 : PetscCall(PetscFunctionListAdd(&LMEList,name,function));
351 9 : PetscFunctionReturn(PETSC_SUCCESS);
352 : }
353 :
354 : /*@C
355 : LMEMonitorRegister - Adds LME monitor routine.
356 :
357 : Not Collective
358 :
359 : Input Parameters:
360 : + name - name of a new monitor routine
361 : . vtype - a PetscViewerType for the output
362 : . format - a PetscViewerFormat for the output
363 : . monitor - monitor routine
364 : . create - creation routine, or NULL
365 : - destroy - destruction routine, or NULL
366 :
367 : Notes:
368 : LMEMonitorRegister() may be called multiple times to add several user-defined monitors.
369 :
370 : Example Usage:
371 : .vb
372 : LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
373 : .ve
374 :
375 : Then, your monitor can be chosen with the procedural interface via
376 : $ LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
377 : or at runtime via the option
378 : $ -lme_monitor_my_monitor
379 :
380 : Level: advanced
381 :
382 : .seealso: LMEMonitorRegisterAll()
383 : @*/
384 18 : PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(LME,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
385 : {
386 18 : char key[PETSC_MAX_PATH_LEN];
387 :
388 18 : PetscFunctionBegin;
389 18 : PetscCall(LMEInitializePackage());
390 18 : PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
391 18 : PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
392 18 : if (create) PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
393 18 : if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
394 18 : PetscFunctionReturn(PETSC_SUCCESS);
395 : }
396 :
397 : /*@
398 : LMEReset - Resets the LME context to the initial state (prior to setup)
399 : and destroys any allocated Vecs and Mats.
400 :
401 : Collective
402 :
403 : Input Parameter:
404 : . lme - linear matrix equation context obtained from LMECreate()
405 :
406 : Level: advanced
407 :
408 : .seealso: LMEDestroy()
409 : @*/
410 8 : PetscErrorCode LMEReset(LME lme)
411 : {
412 8 : PetscFunctionBegin;
413 8 : if (lme) PetscValidHeaderSpecific(lme,LME_CLASSID,1);
414 8 : if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
415 8 : PetscTryTypeMethod(lme,reset);
416 8 : PetscCall(MatDestroy(&lme->A));
417 8 : PetscCall(MatDestroy(&lme->B));
418 8 : PetscCall(MatDestroy(&lme->D));
419 8 : PetscCall(MatDestroy(&lme->E));
420 8 : PetscCall(MatDestroy(&lme->C));
421 8 : PetscCall(MatDestroy(&lme->X));
422 8 : PetscCall(BVDestroy(&lme->V));
423 8 : PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
424 8 : lme->nwork = 0;
425 8 : lme->setupcalled = 0;
426 8 : PetscFunctionReturn(PETSC_SUCCESS);
427 : }
428 :
429 : /*@
430 : LMEDestroy - Destroys the LME context.
431 :
432 : Collective
433 :
434 : Input Parameter:
435 : . lme - linear matrix equation context obtained from LMECreate()
436 :
437 : Level: beginner
438 :
439 : .seealso: LMECreate(), LMESetUp(), LMESolve()
440 : @*/
441 8 : PetscErrorCode LMEDestroy(LME *lme)
442 : {
443 8 : PetscFunctionBegin;
444 8 : if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
445 8 : PetscValidHeaderSpecific(*lme,LME_CLASSID,1);
446 8 : if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
447 8 : PetscCall(LMEReset(*lme));
448 8 : PetscTryTypeMethod(*lme,destroy);
449 8 : PetscCall(LMEMonitorCancel(*lme));
450 8 : PetscCall(PetscHeaderDestroy(lme));
451 8 : PetscFunctionReturn(PETSC_SUCCESS);
452 : }
453 :
454 : /*@
455 : LMESetBV - Associates a basis vectors object to the linear matrix equation solver.
456 :
457 : Collective
458 :
459 : Input Parameters:
460 : + lme - linear matrix equation context obtained from LMECreate()
461 : - bv - the basis vectors object
462 :
463 : Note:
464 : Use LMEGetBV() to retrieve the basis vectors context (for example,
465 : to free it at the end of the computations).
466 :
467 : Level: advanced
468 :
469 : .seealso: LMEGetBV()
470 : @*/
471 0 : PetscErrorCode LMESetBV(LME lme,BV bv)
472 : {
473 0 : PetscFunctionBegin;
474 0 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
475 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,2);
476 0 : PetscCheckSameComm(lme,1,bv,2);
477 0 : PetscCall(PetscObjectReference((PetscObject)bv));
478 0 : PetscCall(BVDestroy(&lme->V));
479 0 : lme->V = bv;
480 0 : PetscFunctionReturn(PETSC_SUCCESS);
481 : }
482 :
483 : /*@
484 : LMEGetBV - Obtain the basis vectors object associated to the matrix
485 : function solver.
486 :
487 : Not Collective
488 :
489 : Input Parameters:
490 : . lme - linear matrix equation context obtained from LMECreate()
491 :
492 : Output Parameter:
493 : . bv - basis vectors context
494 :
495 : Level: advanced
496 :
497 : .seealso: LMESetBV()
498 : @*/
499 7 : PetscErrorCode LMEGetBV(LME lme,BV *bv)
500 : {
501 7 : PetscFunctionBegin;
502 7 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
503 7 : PetscAssertPointer(bv,2);
504 7 : if (!lme->V) {
505 7 : PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
506 7 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
507 7 : PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
508 : }
509 7 : *bv = lme->V;
510 7 : PetscFunctionReturn(PETSC_SUCCESS);
511 : }
|
