Actual source code: lmebasic.c

slepc-3.10.2 2019-02-11
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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:   PetscInt       tabs;
 53:   const char     *eqname[] = {
 54:                    "continuous-time Lyapunov",
 55:                    "continuous-time Sylvester",
 56:                    "generalized Lyapunov",
 57:                    "generalized Sylvester",
 58:                    "Stein",
 59:                    "discrete-time Lyapunov"
 60:   };

 64:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)lme));

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

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

 98:    Collective on LME

100:    Parameter:
101: +  lme - the linear matrix equation context
102: -  viewer - the viewer to display the reason

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

107:    Level: intermediate

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

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

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

134:    Collective on LME

136:    Input Parameters:
137: .  lme - the linear matrix equation context

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

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

163: /*@
164:    LMECreate - Creates the default LME context.

166:    Collective on MPI_Comm

168:    Input Parameter:
169: .  comm - MPI communicator

171:    Output Parameter:
172: .  lme - location to put the LME context

174:    Note:
175:    The default LME type is LMEKRYLOV

177:    Level: beginner

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

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

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

204:   lme->numbermonitors  = 0;

206:   lme->V               = NULL;
207:   lme->nwork           = 0;
208:   lme->work            = NULL;
209:   lme->data            = NULL;

211:   lme->its             = 0;
212:   lme->errest          = 0;
213:   lme->setupcalled     = 0;
214:   lme->reason          = LME_CONVERGED_ITERATING;

216:   *outlme = lme;
217:   return(0);
218: }

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

223:    Logically Collective on LME

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

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

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

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

245:    Level: intermediate

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


258:   PetscObjectTypeCompare((PetscObject)lme,type,&match);
259:   if (match) return(0);

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

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

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

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

276:    Not Collective

278:    Input Parameter:
279: .  lme - the linear matrix equation context

281:    Output Parameter:
282: .  name - name of LME method

284:    Level: intermediate

286: .seealso: LMESetType()
287: @*/
288: PetscErrorCode LMEGetType(LME lme,LMEType *type)
289: {
293:   *type = ((PetscObject)lme)->type_name;
294:   return(0);
295: }

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

300:    Not Collective

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

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

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

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

319:    Level: advanced

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

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

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

337:    Collective on LME

339:    Input Parameter:
340: .  lme - linear matrix equation context obtained from LMECreate()

342:    Level: advanced

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

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

367: /*@
368:    LMEDestroy - Destroys the LME context.

370:    Collective on LME

372:    Input Parameter:
373: .  lme - linear matrix equation context obtained from LMECreate()

375:    Level: beginner

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

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

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

397:    Collective on LME

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

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

407:    Level: advanced

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

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

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

430:    Not Collective

432:    Input Parameters:
433: .  lme - linear matrix equation context obtained from LMECreate()

435:    Output Parameter:
436: .  bv - basis vectors context

438:    Level: advanced

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

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