Actual source code: lmebasic.c

slepc-3.13.3 2020-06-14
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:    LMEReasonView - 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:    Level: intermediate

128: .seealso: LMESetTolerances(), LMEGetIterationNumber()
129: @*/
130: PetscErrorCode LMEReasonView(LME lme,PetscViewer viewer)
131: {
133:   PetscBool      isAscii;

136:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
137:   if (isAscii) {
138:     PetscViewerASCIIAddTab(viewer,((PetscObject)lme)->tablevel);
139:     if (lme->reason > 0) {
140:       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);
141:     } else {
142:       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);
143:     }
144:     PetscViewerASCIISubtractTab(viewer,((PetscObject)lme)->tablevel);
145:   }
146:   return(0);
147: }

149: /*@
150:    LMEReasonViewFromOptions - Processes command line options to determine if/how
151:    the LME converged reason is to be viewed.

153:    Collective on lme

155:    Input Parameter:
156: .  lme - the linear matrix equation context

158:    Level: developer
159: @*/
160: PetscErrorCode LMEReasonViewFromOptions(LME lme)
161: {
162:   PetscErrorCode    ierr;
163:   PetscViewer       viewer;
164:   PetscBool         flg;
165:   static PetscBool  incall = PETSC_FALSE;
166:   PetscViewerFormat format;

169:   if (incall) return(0);
170:   incall = PETSC_TRUE;
171:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)lme),((PetscObject)lme)->options,((PetscObject)lme)->prefix,"-lme_converged_reason",&viewer,&format,&flg);
172:   if (flg) {
173:     PetscViewerPushFormat(viewer,format);
174:     LMEReasonView(lme,viewer);
175:     PetscViewerPopFormat(viewer);
176:     PetscViewerDestroy(&viewer);
177:   }
178:   incall = PETSC_FALSE;
179:   return(0);
180: }

182: /*@
183:    LMECreate - Creates the default LME context.

185:    Collective

187:    Input Parameter:
188: .  comm - MPI communicator

190:    Output Parameter:
191: .  lme - location to put the LME context

193:    Note:
194:    The default LME type is LMEKRYLOV

196:    Level: beginner

198: .seealso: LMESetUp(), LMESolve(), LMEDestroy(), LME
199: @*/
200: PetscErrorCode LMECreate(MPI_Comm comm,LME *outlme)
201: {
203:   LME            lme;

207:   *outlme = 0;
208:   LMEInitializePackage();
209:   SlepcHeaderCreate(lme,LME_CLASSID,"LME","Linear Matrix Equation","LME",comm,LMEDestroy,LMEView);

211:   lme->A               = NULL;
212:   lme->B               = NULL;
213:   lme->D               = NULL;
214:   lme->E               = NULL;
215:   lme->C               = NULL;
216:   lme->X               = NULL;
217:   lme->problem_type    = LME_LYAPUNOV;
218:   lme->max_it          = PETSC_DEFAULT;
219:   lme->ncv             = PETSC_DEFAULT;
220:   lme->tol             = PETSC_DEFAULT;
221:   lme->errorifnotconverged = PETSC_FALSE;

223:   lme->numbermonitors  = 0;

225:   lme->V               = NULL;
226:   lme->nwork           = 0;
227:   lme->work            = NULL;
228:   lme->data            = NULL;

230:   lme->its             = 0;
231:   lme->errest          = 0;
232:   lme->setupcalled     = 0;
233:   lme->reason          = LME_CONVERGED_ITERATING;

235:   *outlme = lme;
236:   return(0);
237: }

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

242:    Logically Collective on lme

244:    Input Parameters:
245: +  lme  - the linear matrix equation context
246: -  type - a known method

248:    Options Database Key:
249: .  -lme_type <method> - Sets the method; use -help for a list
250:     of available methods

252:    Notes:
253:    See "slepc/include/slepclme.h" for available methods. The default
254:    is LMEKRYLOV

256:    Normally, it is best to use the LMESetFromOptions() command and
257:    then set the LME type from the options database rather than by using
258:    this routine.  Using the options database provides the user with
259:    maximum flexibility in evaluating the different available methods.
260:    The LMESetType() routine is provided for those situations where it
261:    is necessary to set the iterative solver independently of the command
262:    line or options database.

264:    Level: intermediate

266: .seealso: LMEType
267: @*/
268: PetscErrorCode LMESetType(LME lme,LMEType type)
269: {
270:   PetscErrorCode ierr,(*r)(LME);
271:   PetscBool      match;


277:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
278:   if (match) return(0);

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

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

286:   lme->setupcalled = 0;
287:   PetscObjectChangeTypeName((PetscObject)lme,type);
288:   (*r)(lme);
289:   return(0);
290: }

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

295:    Not Collective

297:    Input Parameter:
298: .  lme - the linear matrix equation context

300:    Output Parameter:
301: .  name - name of LME method

303:    Level: intermediate

305: .seealso: LMESetType()
306: @*/
307: PetscErrorCode LMEGetType(LME lme,LMEType *type)
308: {
312:   *type = ((PetscObject)lme)->type_name;
313:   return(0);
314: }

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

319:    Not Collective

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

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

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

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

338:    Level: advanced

340: .seealso: LMERegisterAll()
341: @*/
342: PetscErrorCode LMERegister(const char *name,PetscErrorCode (*function)(LME))
343: {

347:   LMEInitializePackage();
348:   PetscFunctionListAdd(&LMEList,name,function);
349:   return(0);
350: }

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

356:    Collective on lme

358:    Input Parameter:
359: .  lme - linear matrix equation context obtained from LMECreate()

361:    Level: advanced

363: .seealso: LMEDestroy()
364: @*/
365: PetscErrorCode LMEReset(LME lme)
366: {

371:   if (!lme) return(0);
372:   if (lme->ops->reset) { (lme->ops->reset)(lme); }
373:   MatDestroy(&lme->A);
374:   MatDestroy(&lme->B);
375:   MatDestroy(&lme->D);
376:   MatDestroy(&lme->E);
377:   MatDestroy(&lme->C);
378:   MatDestroy(&lme->X);
379:   BVDestroy(&lme->V);
380:   VecDestroyVecs(lme->nwork,&lme->work);
381:   lme->nwork = 0;
382:   lme->setupcalled = 0;
383:   return(0);
384: }

386: /*@
387:    LMEDestroy - Destroys the LME context.

389:    Collective on lme

391:    Input Parameter:
392: .  lme - linear matrix equation context obtained from LMECreate()

394:    Level: beginner

396: .seealso: LMECreate(), LMESetUp(), LMESolve()
397: @*/
398: PetscErrorCode LMEDestroy(LME *lme)
399: {

403:   if (!*lme) return(0);
405:   if (--((PetscObject)(*lme))->refct > 0) { *lme = 0; return(0); }
406:   LMEReset(*lme);
407:   if ((*lme)->ops->destroy) { (*(*lme)->ops->destroy)(*lme); }
408:   LMEMonitorCancel(*lme);
409:   PetscHeaderDestroy(lme);
410:   return(0);
411: }

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

416:    Collective on lme

418:    Input Parameters:
419: +  lme - linear matrix equation context obtained from LMECreate()
420: -  bv  - the basis vectors object

422:    Note:
423:    Use LMEGetBV() to retrieve the basis vectors context (for example,
424:    to free it at the end of the computations).

426:    Level: advanced

428: .seealso: LMEGetBV()
429: @*/
430: PetscErrorCode LMESetBV(LME lme,BV bv)
431: {

438:   PetscObjectReference((PetscObject)bv);
439:   BVDestroy(&lme->V);
440:   lme->V = bv;
441:   PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
442:   return(0);
443: }

445: /*@
446:    LMEGetBV - Obtain the basis vectors object associated to the matrix
447:    function solver.

449:    Not Collective

451:    Input Parameters:
452: .  lme - linear matrix equation context obtained from LMECreate()

454:    Output Parameter:
455: .  bv - basis vectors context

457:    Level: advanced

459: .seealso: LMESetBV()
460: @*/
461: PetscErrorCode LMEGetBV(LME lme,BV *bv)
462: {

468:   if (!lme->V) {
469:     BVCreate(PetscObjectComm((PetscObject)lme),&lme->V);
470:     PetscObjectIncrementTabLevel((PetscObject)lme->V,(PetscObject)lme,0);
471:     PetscLogObjectParent((PetscObject)lme,(PetscObject)lme->V);
472:     PetscObjectSetOptions((PetscObject)lme->V,((PetscObject)lme)->options);
473:   }
474:   *bv = lme->V;
475:   return(0);
476: }