Actual source code: lmebasic.c

slepc-3.15.1 2021-05-28
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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
 54: @*/
 55: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
 56: {
 58:   PetscBool      isascii;
 59:   const char     *eqname[] = {
 60:                    "continuous-time Lyapunov",
 61:                    "continuous-time Sylvester",
 62:                    "generalized Lyapunov",
 63:                    "generalized Sylvester",
 64:                    "Stein",
 65:                    "discrete-time Lyapunov"
 66:   };

 70:   if (!viewer) {
 71:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer);
 72:   }

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

100: /*@C
101:    LMEViewFromOptions - View from options

103:    Collective on LME

105:    Input Parameters:
106: +  lme  - the linear matrix equation context
107: .  obj  - optional object
108: -  name - command line option

110:    Level: intermediate

112: .seealso: LMEView(), LMECreate()
113: @*/
114: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
115: {

120:   PetscObjectViewFromOptions((PetscObject)lme,obj,name);
121:   return(0);
122: }
123: /*@C
124:    LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.

126:    Collective on lme

128:    Input Parameters:
129: +  lme - the linear matrix equation context
130: -  viewer - the viewer to display the reason

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

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

141:    Level: intermediate

143: .seealso: LMESetTolerances(), LMEGetIterationNumber(), LMEConvergedReasonViewFromOptions()
144: @*/
145: PetscErrorCode LMEConvergedReasonView(LME lme,PetscViewer viewer)
146: {
147:   PetscErrorCode    ierr;
148:   PetscBool         isAscii;
149:   PetscViewerFormat format;

152:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));
153:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
154:   if (isAscii) {
155:     PetscViewerGetFormat(viewer,&format);
156:     PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
157:     if (lme->reason > 0 && format != PETSC_VIEWER_FAILED) {
158:       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);
159:     } else if (lme->reason <= 0) {
160:       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);
161:     }
162:     PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
163:   }
164:   return(0);
165: }

167: /*@
168:    LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
169:    the LME converged reason is to be viewed.

171:    Collective on lme

173:    Input Parameter:
174: .  lme - the linear matrix equation context

176:    Level: developer

178: .seealso: LMEConvergedReasonView()
179: @*/
180: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
181: {
182:   PetscErrorCode    ierr;
183:   PetscViewer       viewer;
184:   PetscBool         flg;
185:   static PetscBool  incall = PETSC_FALSE;
186:   PetscViewerFormat format;

189:   if (incall) return(0);
190:   incall = PETSC_TRUE;
191:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
192:   if (flg) {
193:     PetscViewerPushFormat(viewer,format);
194:     LMEConvergedReasonView(lme,viewer);
195:     PetscViewerPopFormat(viewer);
196:     PetscViewerDestroy(&viewer);
197:   }
198:   incall = PETSC_FALSE;
199:   return(0);
200: }

202: /*@
203:    LMECreate - Creates the default LME context.

205:    Collective

207:    Input Parameter:
208: .  comm - MPI communicator

210:    Output Parameter:
211: .  lme - location to put the LME context

213:    Note:
214:    The default LME type is LMEKRYLOV

216:    Level: beginner

218: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
219: @*/
220: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
221: {
223:   LME            lme;

227:   *outlme = 0;
228:   LMEInitializePackage();
229:   SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);

231:   lme->A               = NULL;
232:   lme->B               = NULL;
233:   lme->D               = NULL;
234:   lme->E               = NULL;
235:   lme->C               = NULL;
236:   lme->X               = NULL;
237:   lme->problem_type    = LME_LYAPUNOV;
238:   lme->max_it          = PETSC_DEFAULT;
239:   lme->ncv             = PETSC_DEFAULT;
240:   lme->tol             = PETSC_DEFAULT;
241:   lme->errorifnotconverged = PETSC_FALSE;

243:   lme->numbermonitors  = 0;

245:   lme->V               = NULL;
246:   lme->nwork           = 0;
247:   lme->work            = NULL;
248:   lme->data            = NULL;

250:   lme->its             = 0;
251:   lme->errest          = 0;
252:   lme->setupcalled     = 0;
253:   lme->reason          = LME_CONVERGED_ITERATING;

255:   *outlme = lme;
256:   return(0);
257: }

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

262:    Logically Collective on lme

264:    Input Parameters:
265: +  lme  - the linear matrix equation context
266: -  type - a known method

268:    Options Database Key:
269: .  -lme_type <method> - Sets the method; use -help for a list
270:     of available methods

272:    Notes:
273:    See "slepc/include/slepclme.h" for available methods. The default
274:    is LMEKRYLOV

276:    Normally, it is best to use the LMESetFromOptions() command and
277:    then set the LME type from the options database rather than by using
278:    this routine.  Using the options database provides the user with
279:    maximum flexibility in evaluating the different available methods.
280:    The LMESetType() routine is provided for those situations where it
281:    is necessary to set the iterative solver independently of the command
282:    line or options database.

284:    Level: intermediate

286: .seealso: LMEType
287: @*/
288: PetscErrorCode LMESetType(LME lme,LMEType type)
289: {
290:   PetscErrorCode ierr,(*r)(LME);
291:   PetscBool      match;


297:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
298:   if (match) return(0);

300:   PetscFunctionListFind(LMEList,type,&r);
301:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown LME type given: %s",type);

303:   if (lme->ops->destroy) { (*lme->ops->destroy)(lme); }
304:   PetscMemzero(lme->ops,sizeof(struct _LMEOps));

306:   lme->setupcalled = 0;
307:   PetscObjectChangeTypeName((PetscObject)lme,type);
308:   (*r)(lme);
309:   return(0);
310: }

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

315:    Not Collective

317:    Input Parameter:
318: .  lme - the linear matrix equation context

320:    Output Parameter:
321: .  name - name of LME method

323:    Level: intermediate

325: .seealso: LMESetType()
326: @*/
327: PetscErrorCode LMEGetType(LME lme,LMEType *type)
328: {
332:   *type = ((PetscObject)lme)->type_name;
333:   return(0);
334: }

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

339:    Not Collective

341:    Input Parameters:
342: +  name - name of a new user-defined solver
343: -  function - routine to create the solver context

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

348:    Sample usage:
349: .vb
350:     LMERegister("my_solver",MySolverCreate);
351: .ve

353:    Then, your solver can be chosen with the procedural interface via
354: $     LMESetType(lme,"my_solver")
355:    or at runtime via the option
356: $     -lme_type my_solver

358:    Level: advanced

360: .seealso: LMERegisterAll()
361: @*/
362: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
363: {

367:   LMEInitializePackage();
368:   PetscFunctionListAdd(&LMEList,name,function);
369:   return(0);
370: }

372: /*@C
373:    LMEMonitorRegister - Adds LME monitor routine.

375:    Not Collective

377:    Input Parameters:
378: +  name    - name of a new monitor routine
379: .  vtype   - a PetscViewerType for the output
380: .  format  - a PetscViewerFormat for the output
381: .  monitor - monitor routine
382: .  create  - creation routine, or NULL
383: -  destroy - destruction routine, or NULL

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

388:    Sample usage:
389: .vb
390:    LMEMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
391: .ve

393:    Then, your monitor can be chosen with the procedural interface via
394: $      LMEMonitorSetFromOptions(lme,"-lme_monitor_my_monitor","my_monitor",NULL)
395:    or at runtime via the option
396: $      -lme_monitor_my_monitor

398:    Level: advanced

400: .seealso: LMEMonitorRegisterAll()
401: @*/
402: PetscErrorCode LMEMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(LME,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
403: {
404:   char           key[PETSC_MAX_PATH_LEN];

408:   LMEInitializePackage();
409:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
410:   PetscFunctionListAdd(&LMEMonitorList,key,monitor);
411:   if (create)  { PetscFunctionListAdd(&LMEMonitorCreateList,key,create); }
412:   if (destroy) { PetscFunctionListAdd(&LMEMonitorDestroyList,key,destroy); }
413:   return(0);
414: }

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

420:    Collective on lme

422:    Input Parameter:
423: .  lme - linear matrix equation context obtained from LMECreate()

425:    Level: advanced

427: .seealso: LMEDestroy()
428: @*/
429: PetscErrorCode LMEReset(LME lme)
430: {

435:   if (!lme) return(0);
436:   if (lme->ops->reset) { (lme->ops->reset)(lme); }
437:   MatDestroy(&lme->A);
438:   MatDestroy(&lme->B);
439:   MatDestroy(&lme->D);
440:   MatDestroy(&lme->E);
441:   MatDestroy(&lme->C);
442:   MatDestroy(&lme->X);
443:   BVDestroy(&lme->V);
444:   VecDestroyVecs(lme->nwork,&lme->work);
445:   lme->nwork = 0;
446:   lme->setupcalled = 0;
447:   return(0);
448: }

450: /*@C
451:    LMEDestroy - Destroys the LME context.

453:    Collective on lme

455:    Input Parameter:
456: .  lme - linear matrix equation context obtained from LMECreate()

458:    Level: beginner

460: .seealso: LMECreate(), LMESetUp(), LMESolve()
461: @*/
462: PetscErrorCode LMEDestroy(LME *lme)
463: {

467:   if (!*lme) return(0);
469:   if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
470:   LMEReset(*lme);
471:   if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
472:   LMEMonitorCancel(*lme);
473:   PetscHeaderDestroy(lme);
474:   return(0);
475: }

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

480:    Collective on lme

482:    Input Parameters:
483: +  lme - linear matrix equation context obtained from LMECreate()
484: -  bv  - the basis vectors object

486:    Note:
487:    Use LMEGetBV() to retrieve the basis vectors context (for example,
488:    to free it at the end of the computations).

490:    Level: advanced

492: .seealso: LMEGetBV()
493: @*/
494: PetscErrorCode LMESetBV(LME lme,BV bv)
495: {

502:   PetscObjectReference((PetscObject)bv);
503:   BVDestroy(&lme->V);
504:   lme->V = bv;
505:   PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
506:   return(0);
507: }

509: /*@
510:    LMEGetBV - Obtain the basis vectors object associated to the matrix
511:    function solver.

513:    Not Collective

515:    Input Parameters:
516: .  lme - linear matrix equation context obtained from LMECreate()

518:    Output Parameter:
519: .  bv - basis vectors context

521:    Level: advanced

523: .seealso: LMESetBV()
524: @*/
525: PetscErrorCode LMEGetBV(LME lme,BV *bv)
526: {

532:   if (!lme->V) {
533:     BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
534:     PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
535:     PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
536:     PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
537:   }
538:   *bv = lme->V;
539:   return(0);
540: }