Actual source code: lmebasic.c

slepc-main 2024-11-15
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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: /* Logging support */
 17: PetscClassId      LME_CLASSID = 0;
 18: PetscLogEvent     LME_SetUp = 0,LME_Solve = 0,LME_ComputeError = 0;

 20: /* List of registered LME routines */
 21: PetscFunctionList LMEList = NULL;
 22: PetscBool         LMERegisterAllCalled = PETSC_FALSE;

 24: /* List of registered LME monitors */
 25: PetscFunctionList LMEMonitorList              = NULL;
 26: PetscFunctionList LMEMonitorCreateList        = NULL;
 27: PetscFunctionList LMEMonitorDestroyList       = NULL;
 28: PetscBool         LMEMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@
 31:    LMEView - Prints the LME data structure.

 33:    Collective

 35:    Input Parameters:
 36: +  lme - the linear matrix equation solver context
 37: -  viewer - optional visualization context

 39:    Options Database Key:
 40: .  -lme_view -  Calls LMEView() at end of LMESolve()

 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.

 50:    The user can open an alternative visualization context with
 51:    PetscViewerASCIIOpen() - output to a specified file.

 53:    Level: beginner

 55: .seealso: LMECreate()
 56: @*/
 57: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
 58: {
 59:   PetscBool      isascii;
 60:   const char     *eqname[] = {
 61:                    "continuous-time Lyapunov",
 62:                    "continuous-time Sylvester",
 63:                    "generalized Lyapunov",
 64:                    "generalized Sylvester",
 65:                    "Stein",
 66:                    "discrete-time Lyapunov"
 67:   };

 69:   PetscFunctionBegin;
 71:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer));
 73:   PetscCheckSameComm(lme,1,viewer,2);

 75:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
 76:   if (isascii) {
 77:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)lme,viewer));
 78:     PetscCall(PetscViewerASCIIPushTab(viewer));
 79:     PetscTryTypeMethod(lme,view,viewer);
 80:     PetscCall(PetscViewerASCIIPopTab(viewer));
 81:     PetscCall(PetscViewerASCIIPrintf(viewer,"  equation type: %s\n",eqname[lme->problem_type]));
 82:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",lme->ncv));
 83:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",lme->max_it));
 84:     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)lme->tol));
 85:   } else PetscTryTypeMethod(lme,view,viewer);
 86:   PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
 87:   if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
 88:   PetscCall(BVView(lme->V,viewer));
 89:   PetscCall(PetscViewerPopFormat(viewer));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: /*@
 94:    LMEViewFromOptions - View from options

 96:    Collective

 98:    Input Parameters:
 99: +  lme  - the linear matrix equation context
100: .  obj  - optional object
101: -  name - command line option

103:    Level: intermediate

105: .seealso: LMEView(), LMECreate()
106: @*/
107: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
108: {
109:   PetscFunctionBegin;
111:   PetscCall(PetscObjectViewFromOptions((PetscObject)lme,obj,name));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }
114: /*@
115:    LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.

117:    Collective

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:    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

132:    Level: intermediate

134: .seealso: LMESetTolerances(), LMEGetIterationNumber(), LMEConvergedReasonViewFromOptions()
135: @*/
136: PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
137: {
138:   PetscBool         isAscii;
139:   PetscViewerFormat format;

141:   PetscFunctionBegin;
142:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
143:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
144:   if (isAscii) {
145:     PetscCall(PetscViewerGetFormat(viewer,&format));
146:     PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel));
147:     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:     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:     PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel));
150:   }
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: /*@
155:    LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
156:    the LME converged reason is to be viewed.

158:    Collective

160:    Input Parameter:
161: .  lme - the linear matrix equation context

163:    Level: developer

165: .seealso: LMEConvergedReasonView()
166: @*/
167: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
168: {
169:   PetscViewer       viewer;
170:   PetscBool         flg;
171:   static PetscBool  incall = PETSC_FALSE;
172:   PetscViewerFormat format;

174:   PetscFunctionBegin;
175:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
176:   incall = PETSC_TRUE;
177:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg));
178:   if (flg) {
179:     PetscCall(PetscViewerPushFormat(viewer,format));
180:     PetscCall(LMEConvergedReasonView(lme,viewer));
181:     PetscCall(PetscViewerPopFormat(viewer));
182:     PetscCall(PetscViewerDestroy(&viewer));
183:   }
184:   incall = PETSC_FALSE;
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: /*@
189:    LMECreate - Creates the default LME context.

191:    Collective

193:    Input Parameter:
194: .  comm - MPI communicator

196:    Output Parameter:
197: .  outlme - location to put the LME context

199:    Note:
200:    The default LME type is LMEKRYLOV

202:    Level: beginner

204: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
205: @*/
206: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
207: {
208:   LME            lme;

210:   PetscFunctionBegin;
211:   PetscAssertPointer(outlme,2);
212:   PetscCall(LMEInitializePackage());
213:   PetscCall(SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView));

215:   lme->A               = NULL;
216:   lme->B               = NULL;
217:   lme->D               = NULL;
218:   lme->E               = NULL;
219:   lme->C               = NULL;
220:   lme->X               = NULL;
221:   lme->problem_type    = LME_LYAPUNOV;
222:   lme->max_it          = PETSC_DETERMINE;
223:   lme->ncv             = PETSC_DETERMINE;
224:   lme->tol             = PETSC_DETERMINE;
225:   lme->errorifnotconverged = PETSC_FALSE;

227:   lme->numbermonitors  = 0;

229:   lme->V               = NULL;
230:   lme->nwork           = 0;
231:   lme->work            = NULL;
232:   lme->data            = NULL;

234:   lme->its             = 0;
235:   lme->errest          = 0;
236:   lme->setupcalled     = 0;
237:   lme->reason          = LME_CONVERGED_ITERATING;

239:   *outlme = lme;
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*@
244:    LMESetType - Selects the particular solver to be used in the LME object.

246:    Logically Collective

248:    Input Parameters:
249: +  lme  - the linear matrix equation context
250: -  type - a known method

252:    Options Database Key:
253: .  -lme_type <method> - Sets the method; use -help for a list
254:     of available methods

256:    Notes:
257:    See "slepc/include/slepclme.h" for available methods. The default
258:    is LMEKRYLOV

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.

268:    Level: intermediate

270: .seealso: LMEType
271: @*/
272: PetscErrorCode LMESetType(LME lme,LMEType type)
273: {
274:   PetscErrorCode (*r)(LME);
275:   PetscBool      match;

277:   PetscFunctionBegin;
279:   PetscAssertPointer(type,2);

281:   PetscCall(PetscObjectTypeCompare((PetscObject)lme,type,&match));
282:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

284:   PetscCall(PetscFunctionListFind(LMEList,type,&r));
285:   PetscCheck(r,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);

287:   PetscTryTypeMethod(lme,destroy);
288:   PetscCall(PetscMemzero(lme->ops,sizeof(struct _LMEOps)));

290:   lme->setupcalled = 0;
291:   PetscCall(PetscObjectChangeTypeName((PetscObject)lme,type));
292:   PetscCall((*r)(lme));
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:    LMEGetType - Gets the LME type as a string from the LME object.

299:    Not Collective

301:    Input Parameter:
302: .  lme - the linear matrix equation context

304:    Output Parameter:
305: .  type - name of LME method

307:    Level: intermediate

309: .seealso: LMESetType()
310: @*/
311: PetscErrorCode LMEGetType(LME lme,LMEType *type)
312: {
313:   PetscFunctionBegin;
315:   PetscAssertPointer(type,2);
316:   *type = ((PetscObject)lme)->type_name;
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: /*@C
321:    LMERegister - Adds a method to the linear matrix equation solver package.

323:    Not Collective

325:    Input Parameters:
326: +  name - name of a new user-defined solver
327: -  function - routine to create the solver context

329:    Notes:
330:    LMERegister() may be called multiple times to add several user-defined solvers.

332:    Example Usage:
333: .vb
334:     LMERegister("my_solver",MySolverCreate);
335: .ve

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

342:    Level: advanced

344: .seealso: LMERegisterAll()
345: @*/
346: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
347: {
348:   PetscFunctionBegin;
349:   PetscCall(LMEInitializePackage());
350:   PetscCall(PetscFunctionListAdd(&LMEList,name,function));
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

354: /*@C
355:    LMEMonitorRegister - Adds LME monitor routine.

357:    Not Collective

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

367:    Notes:
368:    LMEMonitorRegister() may be called multiple times to add several user-defined monitors.

370:    Example Usage:
371: .vb
372:    LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
373: .ve

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

380:    Level: advanced

382: .seealso: LMEMonitorRegisterAll()
383: @*/
384: 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:   char           key[PETSC_MAX_PATH_LEN];

388:   PetscFunctionBegin;
389:   PetscCall(LMEInitializePackage());
390:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
391:   PetscCall(PetscFunctionListAdd(&LMEMonitorList,key,monitor));
392:   if (create)  PetscCall(PetscFunctionListAdd(&LMEMonitorCreateList,key,create));
393:   if (destroy) PetscCall(PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy));
394:   PetscFunctionReturn(PETSC_SUCCESS);
395: }

397: /*@
398:    LMEReset - Resets the LME context to the initial state (prior to setup)
399:    and destroys any allocated Vecs and Mats.

401:    Collective

403:    Input Parameter:
404: .  lme - linear matrix equation context obtained from LMECreate()

406:    Level: advanced

408: .seealso: LMEDestroy()
409: @*/
410: PetscErrorCode LMEReset(LME lme)
411: {
412:   PetscFunctionBegin;
414:   if (!lme) PetscFunctionReturn(PETSC_SUCCESS);
415:   PetscTryTypeMethod(lme,reset);
416:   PetscCall(MatDestroy(&lme->A));
417:   PetscCall(MatDestroy(&lme->B));
418:   PetscCall(MatDestroy(&lme->D));
419:   PetscCall(MatDestroy(&lme->E));
420:   PetscCall(MatDestroy(&lme->C));
421:   PetscCall(MatDestroy(&lme->X));
422:   PetscCall(BVDestroy(&lme->V));
423:   PetscCall(VecDestroyVecs(lme->nwork,&lme->work));
424:   lme->nwork = 0;
425:   lme->setupcalled = 0;
426:   PetscFunctionReturn(PETSC_SUCCESS);
427: }

429: /*@
430:    LMEDestroy - Destroys the LME context.

432:    Collective

434:    Input Parameter:
435: .  lme - linear matrix equation context obtained from LMECreate()

437:    Level: beginner

439: .seealso: LMECreate(), LMESetUp(), LMESolve()
440: @*/
441: PetscErrorCode LMEDestroy(LME *lme)
442: {
443:   PetscFunctionBegin;
444:   if (!*lme) PetscFunctionReturn(PETSC_SUCCESS);
446:   if (--((PetscObject)*lme)->refct > 0) { *lme = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
447:   PetscCall(LMEReset(*lme));
448:   PetscTryTypeMethod(*lme,destroy);
449:   PetscCall(LMEMonitorCancel(*lme));
450:   PetscCall(PetscHeaderDestroy(lme));
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: /*@
455:    LMESetBV - Associates a basis vectors object to the linear matrix equation solver.

457:    Collective

459:    Input Parameters:
460: +  lme - linear matrix equation context obtained from LMECreate()
461: -  bv  - the basis vectors object

463:    Note:
464:    Use LMEGetBV() to retrieve the basis vectors context (for example,
465:    to free it at the end of the computations).

467:    Level: advanced

469: .seealso: LMEGetBV()
470: @*/
471: PetscErrorCode LMESetBV(LME lme,BV bv)
472: {
473:   PetscFunctionBegin;
476:   PetscCheckSameComm(lme,1,bv,2);
477:   PetscCall(PetscObjectReference((PetscObject)bv));
478:   PetscCall(BVDestroy(&lme->V));
479:   lme->V = bv;
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

483: /*@
484:    LMEGetBV - Obtain the basis vectors object associated to the matrix
485:    function solver.

487:    Not Collective

489:    Input Parameters:
490: .  lme - linear matrix equation context obtained from LMECreate()

492:    Output Parameter:
493: .  bv - basis vectors context

495:    Level: advanced

497: .seealso: LMESetBV()
498: @*/
499: PetscErrorCode LMEGetBV(LME lme,BV *bv)
500: {
501:   PetscFunctionBegin;
503:   PetscAssertPointer(bv,2);
504:   if (!lme->V) {
505:     PetscCall(BVCreate(PetscObjectComm((PetscObject)lme),&lme->V));
506:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0));
507:     PetscCall(PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options));
508:   }
509:   *bv = lme->V;
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }