Actual source code: lmebasic.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepclme.h" I*/

 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

 46: .seealso: PetscViewerASCIIOpen()
 47: @*/
 48: PetscErrorCode LMEView(LME lme,PetscViewer viewer)
 49: {
 51:   PetscBool      isascii;
 52:   const char     *eqname[] = {
 53:                    "continuous-time Lyapunov",
 54:                    "continuous-time Sylvester",
 55:                    "generalized Lyapunov",
 56:                    "generalized Sylvester",
 57:                    "Stein",
 58:                    "discrete-time Lyapunov"
 59:   };

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

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

 93: /*@C
 94:    LMEReasonView - Displays the reason an LME solve converged or diverged.

 96:    Collective on LME

 98:    Parameter:
 99: +  lme - the linear matrix equation context
100: -  viewer - the viewer to display the reason

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

105:    Level: intermediate

107: .seealso: LMESetTolerances(), LMEGetIterationNumber()
108: @*/
109: PetscErrorCode LMEReasonView(LME lme,PetscViewer viewer)
110: {
112:   PetscBool      isAscii;

115:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
116:   if (isAscii) {
117:     PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
118:     if (lme->reason > 0) {
119:       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);
120:     } else {
121:       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);
122:     }
123:     PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
124:   }
125:   return(0);
126: }

128: /*@
129:    LMEReasonViewFromOptions - Processes command line options to determine if/how
130:    the LME converged reason is to be viewed.

132:    Collective on LME

134:    Input Parameters:
135: .  lme - the linear matrix equation context

137:    Level: developer
138: @*/
139: PetscErrorCode LMEReasonViewFromOptions(LME lme)
140: {
141:   PetscErrorCode    ierr;
142:   PetscViewer       viewer;
143:   PetscBool         flg;
144:   static PetscBool  incall = PETSC_FALSE;
145:   PetscViewerFormat format;

148:   if (incall) return(0);
149:   incall = PETSC_TRUE;
150:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
151:   if (flg) {
152:     PetscViewerPushFormat(viewer,format);
153:     LMEReasonView(lme,viewer);
154:     PetscViewerPopFormat(viewer);
155:     PetscViewerDestroy(&viewer);
156:   }
157:   incall = PETSC_FALSE;
158:   return(0);
159: }

161: /*@
162:    LMECreate - Creates the default LME context.

164:    Collective on MPI_Comm

166:    Input Parameter:
167: .  comm - MPI communicator

169:    Output Parameter:
170: .  lme - location to put the LME context

172:    Note:
173:    The default LME type is LMEKRYLOV

175:    Level: beginner

177: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
178: @*/
179: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
180: {
182:   LME            lme;

186:   *outlme = 0;
187:   LMEInitializePackage();
188:   SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);

190:   lme->A               = NULL;
191:   lme->B               = NULL;
192:   lme->D               = NULL;
193:   lme->E               = NULL;
194:   lme->C               = NULL;
195:   lme->X               = NULL;
196:   lme->problem_type    = LME_LYAPUNOV;
197:   lme->max_it          = 0;
198:   lme->ncv             = 0;
199:   lme->tol             = PETSC_DEFAULT;
200:   lme->errorifnotconverged = PETSC_FALSE;

202:   lme->numbermonitors  = 0;

204:   lme->V               = NULL;
205:   lme->nwork           = 0;
206:   lme->work            = NULL;
207:   lme->data            = NULL;

209:   lme->its             = 0;
210:   lme->errest          = 0;
211:   lme->setupcalled     = 0;
212:   lme->reason          = LME_CONVERGED_ITERATING;

214:   *outlme = lme;
215:   return(0);
216: }

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

221:    Logically Collective on LME

223:    Input Parameters:
224: +  lme  - the linear matrix equation context
225: -  type - a known method

227:    Options Database Key:
228: .  -lme_type <method> - Sets the method; use -help for a list
229:     of available methods

231:    Notes:
232:    See "slepc/include/slepclme.h" for available methods. The default
233:    is LMEKRYLOV

235:    Normally, it is best to use the LMESetFromOptions() command and
236:    then set the LME type from the options database rather than by using
237:    this routine.  Using the options database provides the user with
238:    maximum flexibility in evaluating the different available methods.
239:    The LMESetType() routine is provided for those situations where it
240:    is necessary to set the iterative solver independently of the command
241:    line or options database.

243:    Level: intermediate

245: .seealso: LMEType
246: @*/
247: PetscErrorCode LMESetType(LME lme,LMEType type)
248: {
249:   PetscErrorCode ierr,(*r)(LME);
250:   PetscBool      match;


256:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
257:   if (match) return(0);

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

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

265:   lme->setupcalled = 0;
266:   PetscObjectChangeTypeName((PetscObject)lme,type);
267:   (*r)(lme);
268:   return(0);
269: }

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

274:    Not Collective

276:    Input Parameter:
277: .  lme - the linear matrix equation context

279:    Output Parameter:
280: .  name - name of LME method

282:    Level: intermediate

284: .seealso: LMESetType()
285: @*/
286: PetscErrorCode LMEGetType(LME lme,LMEType *type)
287: {
291:   *type = ((PetscObject)lme)->type_name;
292:   return(0);
293: }

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

298:    Not Collective

300:    Input Parameters:
301: +  name - name of a new user-defined solver
302: -  function - routine to create the solver context

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

307:    Sample usage:
308: .vb
309:     LMERegister("my_solver",MySolverCreate);
310: .ve

312:    Then, your solver can be chosen with the procedural interface via
313: $     LMESetType(lme,"my_solver")
314:    or at runtime via the option
315: $     -lme_type my_solver

317:    Level: advanced

319: .seealso: LMERegisterAll()
320: @*/
321: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
322: {

326:   LMEInitializePackage();
327:   PetscFunctionListAdd(&LMEList,name,function);
328:   return(0);
329: }

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

335:    Collective on LME

337:    Input Parameter:
338: .  lme - linear matrix equation context obtained from LMECreate()

340:    Level: advanced

342: .seealso: LMEDestroy()
343: @*/
344: PetscErrorCode LMEReset(LME lme)
345: {

350:   if (!lme) return(0);
351:   if (lme->ops->reset) { (lme->ops->reset)(lme); }
352:   MatDestroy(&lme->A);
353:   MatDestroy(&lme->B);
354:   MatDestroy(&lme->D);
355:   MatDestroy(&lme->E);
356:   MatDestroy(&lme->C);
357:   MatDestroy(&lme->X);
358:   BVDestroy(&lme->V);
359:   VecDestroyVecs(lme->nwork,&lme->work);
360:   lme->nwork = 0;
361:   lme->setupcalled = 0;
362:   return(0);
363: }

365: /*@
366:    LMEDestroy - Destroys the LME context.

368:    Collective on LME

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

373:    Level: beginner

375: .seealso: LMECreate(), LMESetUp(), LMESolve()
376: @*/
377: PetscErrorCode LMEDestroy(LME *lme)
378: {

382:   if (!*lme) return(0);
384:   if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
385:   LMEReset(*lme);
386:   if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
387:   LMEMonitorCancel(*lme);
388:   PetscHeaderDestroy(lme);
389:   return(0);
390: }

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

395:    Collective on LME

397:    Input Parameters:
398: +  lme - linear matrix equation context obtained from LMECreate()
399: -  bv  - the basis vectors object

401:    Note:
402:    Use LMEGetBV() to retrieve the basis vectors context (for example,
403:    to free it at the end of the computations).

405:    Level: advanced

407: .seealso: LMEGetBV()
408: @*/
409: PetscErrorCode LMESetBV(LME lme,BV bv)
410: {

417:   PetscObjectReference((PetscObject)bv);
418:   BVDestroy(&lme->V);
419:   lme->V = bv;
420:   PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
421:   return(0);
422: }

424: /*@
425:    LMEGetBV - Obtain the basis vectors object associated to the matrix
426:    function solver.

428:    Not Collective

430:    Input Parameters:
431: .  lme - linear matrix equation context obtained from LMECreate()

433:    Output Parameter:
434: .  bv - basis vectors context

436:    Level: advanced

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

447:   if (!lme->V) {
448:     BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
449:     PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
450:     PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
451:     PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
452:   }
453:   *bv = lme->V;
454:   return(0);
455: }