Actual source code: lmebasic.c

slepc-3.17.2 2022-08-09
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 = 0;
 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: /*@C
 31:    LMEView - Prints the LME data structure.

 33:    Collective on lme

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

 70:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer);

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

 96: /*@C
 97:    LMEViewFromOptions - View from options

 99:    Collective on LME

101:    Input Parameters:
102: +  lme  - the linear matrix equation context
103: .  obj  - optional object
104: -  name - command line option

106:    Level: intermediate

108: .seealso: LMEView(), LMECreate()
109: @*/
110: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
111: {
113:   PetscObjectViewFromOptions((PetscObject)lme,obj,name);
114:   PetscFunctionReturn(0);
115: }
116: /*@C
117:    LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.

119:    Collective on lme

121:    Input Parameters:
122: +  lme - the linear matrix equation context
123: -  viewer - the viewer to display the reason

125:    Options Database Keys:
126: .  -lme_converged_reason - print reason for convergence, and number of iterations

128:    Note:
129:    To change the format of the output call PetscViewerPushFormat(viewer,format) before
130:    this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
131:    display a reason if it fails. The latter can be set in the command line with
132:    -lme_converged_reason ::failed

134:    Level: intermediate

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

143:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
144:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
145:   if (isAscii) {
146:     PetscViewerGetFormat(viewer,&format);
147:     PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
148:     if (lme->reason > 0 && format != PETSC_VIEWER_FAILED) 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);
149:     else if (lme->reason <= 0) 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);
150:     PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
151:   }
152:   PetscFunctionReturn(0);
153: }

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

159:    Collective on lme

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

164:    Level: developer

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

175:   if (incall) PetscFunctionReturn(0);
176:   incall = PETSC_TRUE;
177:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
178:   if (flg) {
179:     PetscViewerPushFormat(viewer,format);
180:     LMEConvergedReasonView(lme,viewer);
181:     PetscViewerPopFormat(viewer);
182:     PetscViewerDestroy(&viewer);
183:   }
184:   incall = PETSC_FALSE;
185:   PetscFunctionReturn(0);
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;

211:   *outlme = 0;
212:   LMEInitializePackage();
213:   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_DEFAULT;
223:   lme->ncv             = PETSC_DEFAULT;
224:   lme->tol             = PETSC_DEFAULT;
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(0);
241: }

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

246:    Logically Collective on lme

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;


280:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
281:   if (match) PetscFunctionReturn(0);

283:   PetscFunctionListFind(LMEList,type,&r);

286:   if (lme->ops->destroy) (*lme->ops->destroy)(lme);
287:   PetscMemzero(lme->ops,sizeof(struct _LMEOps));

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

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

298:    Not Collective

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

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

306:    Level: intermediate

308: .seealso: LMESetType()
309: @*/
310: PetscErrorCode LMEGetType(LME lme,LMEType *type)
311: {
314:   *type = ((PetscObject)lme)->type_name;
315:   PetscFunctionReturn(0);
316: }

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

321:    Not Collective

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

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

330:    Sample usage:
331: .vb
332:     LMERegister("my_solver",MySolverCreate);
333: .ve

335:    Then, your solver can be chosen with the procedural interface via
336: $     LMESetType(lme,"my_solver")
337:    or at runtime via the option
338: $     -lme_type my_solver

340:    Level: advanced

342: .seealso: LMERegisterAll()
343: @*/
344: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
345: {
346:   LMEInitializePackage();
347:   PetscFunctionListAdd(&LMEList,name,function);
348:   PetscFunctionReturn(0);
349: }

351: /*@C
352:    LMEMonitorRegister - Adds LME monitor routine.

354:    Not Collective

356:    Input Parameters:
357: +  name    - name of a new monitor routine
358: .  vtype   - a PetscViewerType for the output
359: .  format  - a PetscViewerFormat for the output
360: .  monitor - monitor routine
361: .  create  - creation routine, or NULL
362: -  destroy - destruction routine, or NULL

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

367:    Sample usage:
368: .vb
369:    LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
370: .ve

372:    Then, your monitor can be chosen with the procedural interface via
373: $      LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
374:    or at runtime via the option
375: $      -lme_monitor_my_monitor

377:    Level: advanced

379: .seealso: LMEMonitorRegisterAll()
380: @*/
381: PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(LME,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
382: {
383:   char           key[PETSC_MAX_PATH_LEN];

385:   LMEInitializePackage();
386:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
387:   PetscFunctionListAdd(&LMEMonitorList,key,monitor);
388:   if (create)  PetscFunctionListAdd(&LMEMonitorCreateList,key,create);
389:   if (destroy) PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy);
390:   PetscFunctionReturn(0);
391: }

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

397:    Collective on lme

399:    Input Parameter:
400: .  lme - linear matrix equation context obtained from LMECreate()

402:    Level: advanced

404: .seealso: LMEDestroy()
405: @*/
406: PetscErrorCode LMEReset(LME lme)
407: {
409:   if (!lme) PetscFunctionReturn(0);
410:   if (lme->ops->reset) (lme->ops->reset)(lme);
411:   MatDestroy(&lme->A);
412:   MatDestroy(&lme->B);
413:   MatDestroy(&lme->D);
414:   MatDestroy(&lme->E);
415:   MatDestroy(&lme->C);
416:   MatDestroy(&lme->X);
417:   BVDestroy(&lme->V);
418:   VecDestroyVecs(lme->nwork,&lme->work);
419:   lme->nwork = 0;
420:   lme->setupcalled = 0;
421:   PetscFunctionReturn(0);
422: }

424: /*@C
425:    LMEDestroy - Destroys the LME context.

427:    Collective on lme

429:    Input Parameter:
430: .  lme - linear matrix equation context obtained from LMECreate()

432:    Level: beginner

434: .seealso: LMECreate(), LMESetUp(), LMESolve()
435: @*/
436: PetscErrorCode LMEDestroy(LME *lme)
437: {
438:   if (!*lme) PetscFunctionReturn(0);
440:   if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; PetscFunctionReturn(0); }
441:   LMEReset(*lme);
442:   if ((*lme)->ops->destroy) (*(*lme)->ops->destroy)(*lme);
443:   LMEMonitorCancel(*lme);
444:   PetscHeaderDestroy(lme);
445:   PetscFunctionReturn(0);
446: }

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

451:    Collective on lme

453:    Input Parameters:
454: +  lme - linear matrix equation context obtained from LMECreate()
455: -  bv  - the basis vectors object

457:    Note:
458:    Use LMEGetBV() to retrieve the basis vectors context (for example,
459:    to free it at the end of the computations).

461:    Level: advanced

463: .seealso: LMEGetBV()
464: @*/
465: PetscErrorCode LMESetBV(LME lme,BV bv)
466: {
470:   PetscObjectReference((PetscObject)bv);
471:   BVDestroy(&lme->V);
472:   lme->V = bv;
473:   PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
474:   PetscFunctionReturn(0);
475: }

477: /*@
478:    LMEGetBV - Obtain the basis vectors object associated to the matrix
479:    function solver.

481:    Not Collective

483:    Input Parameters:
484: .  lme - linear matrix equation context obtained from LMECreate()

486:    Output Parameter:
487: .  bv - basis vectors context

489:    Level: advanced

491: .seealso: LMESetBV()
492: @*/
493: PetscErrorCode LMEGetBV(LME lme,BV *bv)
494: {
497:   if (!lme->V) {
498:     BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
499:     PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
500:     PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
501:     PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
502:   }
503:   *bv = lme->V;
504:   PetscFunctionReturn(0);
505: }