Actual source code: lmebasic.c

slepc-3.14.0 2020-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: PetscFunctionList LMEList = 0;
 17: PetscBool         LMERegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      LME_CLASSID = 0;
 19: PetscLogEvent     LME_SetUp = 0,LME_Solve = 0,LME_ComputeError = 0;

 21: /*@C
 22:    LMEView - Prints the LME data structure.

 24:    Collective on lme

 26:    Input Parameters:
 27: +  lme - the linear matrix equation solver context
 28: -  viewer - optional visualization context

 30:    Options Database Key:
 31: .  -lme_view -  Calls LMEView() at end of LMESolve()

 33:    Note:
 34:    The available visualization contexts include
 35: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 36: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 37:          output where only the first processor opens
 38:          the file.  All other processors send their
 39:          data to the first processor to print.

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

 44:    Level: beginner
 45: @*/
 46: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
 47: {
 49:   PetscBool      isascii;
 50:   const char     *eqname[] = {
 51:                    "continuous-time Lyapunov",
 52:                    "continuous-time Sylvester",
 53:                    "generalized Lyapunov",
 54:                    "generalized Sylvester",
 55:                    "Stein",
 56:                    "discrete-time Lyapunov"
 57:   };

 61:   if (!viewer) {
 62:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lme),&viewer);
 63:   }

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

 91: /*@C
 92:    LMEViewFromOptions - View from options

 94:    Collective on LME

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

101:    Level: intermediate

103: .seealso: LMEView(), LMECreate()
104: @*/
105: PetscErrorCode LMEViewFromOptions(LME lme,PetscObject obj,const char name[])
106: {

111:   PetscObjectViewFromOptions((PetscObject)lme,obj,name);
112:   return(0);
113: }
114: /*@C
115:    LMEConvergedReasonView - Displays the reason an LME solve converged or diverged.

117:    Collective on lme

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:   PetscErrorCode    ierr;
139:   PetscBool         isAscii;
140:   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) {
149:       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);
150:     } else if (lme->reason <= 0) {
151:       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);
152:     }
153:     PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
154:   }
155:   return(0);
156: }

158: /*@
159:    LMEConvergedReasonViewFromOptions - Processes command line options to determine if/how
160:    the LME converged reason is to be viewed.

162:    Collective on lme

164:    Input Parameter:
165: .  lme - the linear matrix equation context

167:    Level: developer

169: .seealso: LMEConvergedReasonView()
170: @*/
171: PetscErrorCode LMEConvergedReasonViewFromOptions(LME lme)
172: {
173:   PetscErrorCode    ierr;
174:   PetscViewer       viewer;
175:   PetscBool         flg;
176:   static PetscBool  incall = PETSC_FALSE;
177:   PetscViewerFormat format;

180:   if (incall) return(0);
181:   incall = PETSC_TRUE;
182:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
183:   if (flg) {
184:     PetscViewerPushFormat(viewer,format);
185:     LMEConvergedReasonView(lme,viewer);
186:     PetscViewerPopFormat(viewer);
187:     PetscViewerDestroy(&viewer);
188:   }
189:   incall = PETSC_FALSE;
190:   return(0);
191: }

193: /*@
194:    LMECreate - Creates the default LME context.

196:    Collective

198:    Input Parameter:
199: .  comm - MPI communicator

201:    Output Parameter:
202: .  lme - location to put the LME context

204:    Note:
205:    The default LME type is LMEKRYLOV

207:    Level: beginner

209: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
210: @*/
211: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
212: {
214:   LME            lme;

218:   *outlme = 0;
219:   LMEInitializePackage();
220:   SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);

222:   lme->A               = NULL;
223:   lme->B               = NULL;
224:   lme->D               = NULL;
225:   lme->E               = NULL;
226:   lme->C               = NULL;
227:   lme->X               = NULL;
228:   lme->problem_type    = LME_LYAPUNOV;
229:   lme->max_it          = PETSC_DEFAULT;
230:   lme->ncv             = PETSC_DEFAULT;
231:   lme->tol             = PETSC_DEFAULT;
232:   lme->errorifnotconverged = PETSC_FALSE;

234:   lme->numbermonitors  = 0;

236:   lme->V               = NULL;
237:   lme->nwork           = 0;
238:   lme->work            = NULL;
239:   lme->data            = NULL;

241:   lme->its             = 0;
242:   lme->errest          = 0;
243:   lme->setupcalled     = 0;
244:   lme->reason          = LME_CONVERGED_ITERATING;

246:   *outlme = lme;
247:   return(0);
248: }

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

253:    Logically Collective on lme

255:    Input Parameters:
256: +  lme  - the linear matrix equation context
257: -  type - a known method

259:    Options Database Key:
260: .  -lme_type <method> - Sets the method; use -help for a list
261:     of available methods

263:    Notes:
264:    See "slepc/include/slepclme.h" for available methods. The default
265:    is LMEKRYLOV

267:    Normally, it is best to use the LMESetFromOptions() command and
268:    then set the LME type from the options database rather than by using
269:    this routine.  Using the options database provides the user with
270:    maximum flexibility in evaluating the different available methods.
271:    The LMESetType() routine is provided for those situations where it
272:    is necessary to set the iterative solver independently of the command
273:    line or options database.

275:    Level: intermediate

277: .seealso: LMEType
278: @*/
279: PetscErrorCode LMESetType(LME lme,LMEType type)
280: {
281:   PetscErrorCode ierr,(*r)(LME);
282:   PetscBool      match;


288:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
289:   if (match) return(0);

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

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

297:   lme->setupcalled = 0;
298:   PetscObjectChangeTypeName((PetscObject)lme,type);
299:   (*r)(lme);
300:   return(0);
301: }

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

306:    Not Collective

308:    Input Parameter:
309: .  lme - the linear matrix equation context

311:    Output Parameter:
312: .  name - name of LME method

314:    Level: intermediate

316: .seealso: LMESetType()
317: @*/
318: PetscErrorCode LMEGetType(LME lme,LMEType *type)
319: {
323:   *type = ((PetscObject)lme)->type_name;
324:   return(0);
325: }

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

330:    Not Collective

332:    Input Parameters:
333: +  name - name of a new user-defined solver
334: -  function - routine to create the solver context

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

339:    Sample usage:
340: .vb
341:     LMERegister("my_solver",MySolverCreate);
342: .ve

344:    Then, your solver can be chosen with the procedural interface via
345: $     LMESetType(lme,"my_solver")
346:    or at runtime via the option
347: $     -lme_type my_solver

349:    Level: advanced

351: .seealso: LMERegisterAll()
352: @*/
353: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
354: {

358:   LMEInitializePackage();
359:   PetscFunctionListAdd(&LMEList,name,function);
360:   return(0);
361: }

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

367:    Collective on lme

369:    Input Parameter:
370: .  lme - linear matrix equation context obtained from LMECreate()

372:    Level: advanced

374: .seealso: LMEDestroy()
375: @*/
376: PetscErrorCode LMEReset(LME lme)
377: {

382:   if (!lme) return(0);
383:   if (lme->ops->reset) { (lme->ops->reset)(lme); }
384:   MatDestroy(&lme->A);
385:   MatDestroy(&lme->B);
386:   MatDestroy(&lme->D);
387:   MatDestroy(&lme->E);
388:   MatDestroy(&lme->C);
389:   MatDestroy(&lme->X);
390:   BVDestroy(&lme->V);
391:   VecDestroyVecs(lme->nwork,&lme->work);
392:   lme->nwork = 0;
393:   lme->setupcalled = 0;
394:   return(0);
395: }

397: /*@
398:    LMEDestroy - Destroys the LME context.

400:    Collective on lme

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

405:    Level: beginner

407: .seealso: LMECreate(), LMESetUp(), LMESolve()
408: @*/
409: PetscErrorCode LMEDestroy(LME *lme)
410: {

414:   if (!*lme) return(0);
416:   if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
417:   LMEReset(*lme);
418:   if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
419:   LMEMonitorCancel(*lme);
420:   PetscHeaderDestroy(lme);
421:   return(0);
422: }

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

427:    Collective on lme

429:    Input Parameters:
430: +  lme - linear matrix equation context obtained from LMECreate()
431: -  bv  - the basis vectors object

433:    Note:
434:    Use LMEGetBV() to retrieve the basis vectors context (for example,
435:    to free it at the end of the computations).

437:    Level: advanced

439: .seealso: LMEGetBV()
440: @*/
441: PetscErrorCode LMESetBV(LME lme,BV bv)
442: {

449:   PetscObjectReference((PetscObject)bv);
450:   BVDestroy(&lme->V);
451:   lme->V = bv;
452:   PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
453:   return(0);
454: }

456: /*@
457:    LMEGetBV - Obtain the basis vectors object associated to the matrix
458:    function solver.

460:    Not Collective

462:    Input Parameters:
463: .  lme - linear matrix equation context obtained from LMECreate()

465:    Output Parameter:
466: .  bv - basis vectors context

468:    Level: advanced

470: .seealso: LMESetBV()
471: @*/
472: PetscErrorCode LMEGetBV(LME lme,BV *bv)
473: {

479:   if (!lme->V) {
480:     BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
481:     PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
482:     PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
483:     PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
484:   }
485:   *bv = lme->V;
486:   return(0);
487: }