Actual source code: epsmon.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:    EPS routines related to monitors
 12: */

 14: #include <slepc/private/epsimpl.h>   /*I "slepceps.h" I*/
 15: #include <petscdraw.h>

 17: /*
 18:    Runs the user provided monitor routines, if any.
 19: */
 20: PetscErrorCode EPSMonitor(EPS eps,PetscInt it,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest)
 21: {
 23:   PetscInt       i,n = eps->numbermonitors;

 26:   for (i=0;i<n;i++) {
 27:     (*eps->monitor[i])(eps,it,nconv,eigr,eigi,errest,nest,eps->monitorcontext[i]);
 28:   }
 29:   return(0);
 30: }

 32: /*@C
 33:    EPSMonitorSet - Sets an ADDITIONAL function to be called at every
 34:    iteration to monitor the error estimates for each requested eigenpair.

 36:    Logically Collective on EPS

 38:    Input Parameters:
 39: +  eps     - eigensolver context obtained from EPSCreate()
 40: .  monitor - pointer to function (if this is NULL, it turns off monitoring)
 41: .  mctx    - [optional] context for private data for the
 42:              monitor routine (use NULL if no context is desired)
 43: -  monitordestroy - [optional] routine that frees monitor context (may be NULL)

 45:    Calling Sequence of monitor:
 46: $   monitor(EPS eps,int its,int nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal* errest,int nest,void *mctx)

 48: +  eps    - eigensolver context obtained from EPSCreate()
 49: .  its    - iteration number
 50: .  nconv  - number of converged eigenpairs
 51: .  eigr   - real part of the eigenvalues
 52: .  eigi   - imaginary part of the eigenvalues
 53: .  errest - relative error estimates for each eigenpair
 54: .  nest   - number of error estimates
 55: -  mctx   - optional monitoring context, as set by EPSMonitorSet()

 57:    Options Database Keys:
 58: +    -eps_monitor        - print only the first error estimate
 59: .    -eps_monitor_all    - print error estimates at each iteration
 60: .    -eps_monitor_conv   - print the eigenvalue approximations only when
 61:       convergence has been reached
 62: .    -eps_monitor_lg     - sets line graph monitor for the first unconverged
 63:       approximate eigenvalue
 64: .    -eps_monitor_lg_all - sets line graph monitor for all unconverged
 65:       approximate eigenvalues
 66: -    -eps_monitor_cancel - cancels all monitors that have been hardwired into
 67:       a code by calls to EPSMonitorSet(), but does not cancel those set via
 68:       the options database.

 70:    Notes:
 71:    Several different monitoring routines may be set by calling
 72:    EPSMonitorSet() multiple times; all will be called in the
 73:    order in which they were set.

 75:    Level: intermediate

 77: .seealso: EPSMonitorFirst(), EPSMonitorAll(), EPSMonitorCancel()
 78: @*/
 79: PetscErrorCode EPSMonitorSet(EPS eps,PetscErrorCode (*monitor)(EPS,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
 80: {
 83:   if (eps->numbermonitors >= MAXEPSMONITORS) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Too many EPS monitors set");
 84:   eps->monitor[eps->numbermonitors]           = monitor;
 85:   eps->monitorcontext[eps->numbermonitors]    = (void*)mctx;
 86:   eps->monitordestroy[eps->numbermonitors++]  = monitordestroy;
 87:   return(0);
 88: }

 90: /*@
 91:    EPSMonitorCancel - Clears all monitors for an EPS object.

 93:    Logically Collective on EPS

 95:    Input Parameters:
 96: .  eps - eigensolver context obtained from EPSCreate()

 98:    Options Database Key:
 99: .    -eps_monitor_cancel - Cancels all monitors that have been hardwired
100:       into a code by calls to EPSMonitorSet(),
101:       but does not cancel those set via the options database.

103:    Level: intermediate

105: .seealso: EPSMonitorSet()
106: @*/
107: PetscErrorCode EPSMonitorCancel(EPS eps)
108: {
110:   PetscInt       i;

114:   for (i=0; i<eps->numbermonitors; i++) {
115:     if (eps->monitordestroy[i]) {
116:       (*eps->monitordestroy[i])(&eps->monitorcontext[i]);
117:     }
118:   }
119:   eps->numbermonitors = 0;
120:   return(0);
121: }

123: /*@C
124:    EPSGetMonitorContext - Gets the monitor context, as set by
125:    EPSMonitorSet() for the FIRST monitor only.

127:    Not Collective

129:    Input Parameter:
130: .  eps - eigensolver context obtained from EPSCreate()

132:    Output Parameter:
133: .  ctx - monitor context

135:    Level: intermediate

137: .seealso: EPSMonitorSet()
138: @*/
139: PetscErrorCode EPSGetMonitorContext(EPS eps,void **ctx)
140: {
143:   *ctx = eps->monitorcontext[0];
144:   return(0);
145: }

147: /*@C
148:    EPSMonitorAll - Print the current approximate values and
149:    error estimates at each iteration of the eigensolver.

151:    Collective on EPS

153:    Input Parameters:
154: +  eps    - eigensolver context
155: .  its    - iteration number
156: .  nconv  - number of converged eigenpairs so far
157: .  eigr   - real part of the eigenvalues
158: .  eigi   - imaginary part of the eigenvalues
159: .  errest - error estimates
160: .  nest   - number of error estimates to display
161: -  vf     - viewer and format for monitoring

163:    Level: intermediate

165: .seealso: EPSMonitorSet(), EPSMonitorFirst(), EPSMonitorConverged()
166: @*/
167: PetscErrorCode EPSMonitorAll(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,PetscViewerAndFormat *vf)
168: {
170:   PetscInt       i;
171:   PetscScalar    er,ei;
172:   PetscViewer    viewer;

177:   viewer = vf->viewer;
179:   PetscViewerPushFormat(viewer,vf->format);
180:   PetscViewerASCIIAddTab(viewer,((PetscObject)eps)->tablevel);
181:   if (its==1 && ((PetscObject)eps)->prefix) {
182:     PetscViewerASCIIPrintf(viewer,"  Eigenvalue approximations and residual norms for %s solve.\n",((PetscObject)eps)->prefix);
183:   }
184:   PetscViewerASCIIPrintf(viewer,"%3D EPS nconv=%D Values (Errors)",its,nconv);
185:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
186:   for (i=0;i<nest;i++) {
187:     er = eigr[i]; ei = eigi[i];
188:     STBackTransform(eps->st,1,&er,&ei);
189: #if defined(PETSC_USE_COMPLEX)
190:     PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
191: #else
192:     PetscViewerASCIIPrintf(viewer," %g",(double)er);
193:     if (ei!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
194: #endif
195:     PetscViewerASCIIPrintf(viewer," (%10.8e)",(double)errest[i]);
196:   }
197:   PetscViewerASCIIPrintf(viewer,"\n");
198:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
199:   PetscViewerASCIISubtractTab(viewer,((PetscObject)eps)->tablevel);
200:   PetscViewerPopFormat(viewer);
201:   return(0);
202: }

204: /*@C
205:    EPSMonitorFirst - Print the first approximate value and
206:    error estimate at each iteration of the eigensolver.

208:    Collective on EPS

210:    Input Parameters:
211: +  eps    - eigensolver context
212: .  its    - iteration number
213: .  nconv  - number of converged eigenpairs so far
214: .  eigr   - real part of the eigenvalues
215: .  eigi   - imaginary part of the eigenvalues
216: .  errest - error estimates
217: .  nest   - number of error estimates to display
218: -  vf     - viewer and format for monitoring

220:    Level: intermediate

222: .seealso: EPSMonitorSet(), EPSMonitorAll(), EPSMonitorConverged()
223: @*/
224: PetscErrorCode EPSMonitorFirst(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,PetscViewerAndFormat *vf)
225: {
227:   PetscScalar    er,ei;
228:   PetscViewer    viewer;

233:   viewer = vf->viewer;
235:   if (its==1 && ((PetscObject)eps)->prefix) {
236:     PetscViewerASCIIPrintf(viewer,"  Eigenvalue approximations and residual norms for %s solve.\n",((PetscObject)eps)->prefix);
237:   }
238:   if (nconv<nest) {
239:     PetscViewerPushFormat(viewer,vf->format);
240:     PetscViewerASCIIAddTab(viewer,((PetscObject)eps)->tablevel);
241:     PetscViewerASCIIPrintf(viewer,"%3D EPS nconv=%D first unconverged value (error)",its,nconv);
242:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
243:     er = eigr[nconv]; ei = eigi[nconv];
244:     STBackTransform(eps->st,1,&er,&ei);
245: #if defined(PETSC_USE_COMPLEX)
246:     PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
247: #else
248:     PetscViewerASCIIPrintf(viewer," %g",(double)er);
249:     if (ei!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
250: #endif
251:     PetscViewerASCIIPrintf(viewer," (%10.8e)\n",(double)errest[nconv]);
252:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
253:     PetscViewerASCIISubtractTab(viewer,((PetscObject)eps)->tablevel);
254:     PetscViewerPopFormat(viewer);
255:   }
256:   return(0);
257: }

259: /*@C
260:    EPSMonitorConverged - Print the approximate values and
261:    error estimates as they converge.

263:    Collective on EPS

265:    Input Parameters:
266: +  eps    - eigensolver context
267: .  its    - iteration number
268: .  nconv  - number of converged eigenpairs so far
269: .  eigr   - real part of the eigenvalues
270: .  eigi   - imaginary part of the eigenvalues
271: .  errest - error estimates
272: .  nest   - number of error estimates to display
273: -  ctx    - monitor context

275:    Level: intermediate

277: .seealso: EPSMonitorSet(), EPSMonitorFirst(), EPSMonitorAll()
278: @*/
279: PetscErrorCode EPSMonitorConverged(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,SlepcConvMonitor ctx)
280: {
282:   PetscInt       i;
283:   PetscScalar    er,ei;
284:   PetscViewer    viewer;

289:   viewer = ctx->viewer;
291:   if (its==1 && ((PetscObject)eps)->prefix) {
292:     PetscViewerASCIIPrintf(viewer,"  Convergence history for %s solve.\n",((PetscObject)eps)->prefix);
293:   }
294:   if (its==1) ctx->oldnconv = 0;
295:   if (ctx->oldnconv!=nconv) {
296:     PetscViewerPushFormat(viewer,ctx->format);
297:     PetscViewerASCIIAddTab(viewer,((PetscObject)eps)->tablevel);
298:     for (i=ctx->oldnconv;i<nconv;i++) {
299:       PetscViewerASCIIPrintf(viewer,"%3D EPS converged value (error) #%D",its,i);
300:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
301:       er = eigr[i]; ei = eigi[i];
302:       STBackTransform(eps->st,1,&er,&ei);
303: #if defined(PETSC_USE_COMPLEX)
304:       PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
305: #else
306:       PetscViewerASCIIPrintf(viewer," %g",(double)er);
307:       if (ei!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
308: #endif
309:       PetscViewerASCIIPrintf(viewer," (%10.8e)\n",(double)errest[i]);
310:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
311:     }
312:     PetscViewerASCIISubtractTab(viewer,((PetscObject)eps)->tablevel);
313:     PetscViewerPopFormat(viewer);
314:     ctx->oldnconv = nconv;
315:   }
316:   return(0);
317: }

319: /*@C
320:    EPSMonitorLGCreate - Creates a line graph context for use with
321:    EPS to monitor convergence.

323:    Collective on MPI_Comm

325:    Input Parameters:
326: +  comm - communicator context
327: .  host - the X display to open, or null for the local machine
328: .  label - the title to put in the title bar
329: .  x, y - the screen coordinates of the upper left coordinate of
330:           the window
331: -  m, n - the screen width and height in pixels

333:    Output Parameter:
334: .  lgctx - the drawing context

336:    Options Database Keys:
337: +  -eps_monitor_lg - Sets line graph monitor for the first residual
338: -  -eps_monitor_lg_all - Sets line graph monitor for all residuals

340:    Notes:
341:    Use PetscDrawLGDestroy() to destroy this line graph.

343:    Level: intermediate

345: .seealso: EPSMonitorSet()
346: @*/
347: PetscErrorCode EPSMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx)
348: {
349:   PetscDraw      draw;
350:   PetscDrawLG    lg;

354:   PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
355:   PetscDrawSetFromOptions(draw);
356:   PetscDrawLGCreate(draw,1,&lg);
357:   PetscDrawLGSetFromOptions(lg);
358:   PetscDrawDestroy(&draw);
359:   *lgctx = lg;
360:   return(0);
361: }

363: PetscErrorCode EPSMonitorLG(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
364: {
365:   PetscDrawLG    lg = (PetscDrawLG)ctx;
366:   PetscReal      x,y;

371:   if (its==1) {
372:     PetscDrawLGReset(lg);
373:     PetscDrawLGSetDimension(lg,1);
374:     PetscDrawLGSetLimits(lg,1,1.0,PetscLog10Real(eps->tol)-2,0.0);
375:   }
376:   x = (PetscReal)its;
377:   if (errest[nconv] > 0.0) y = PetscLog10Real(errest[nconv]);
378:   else y = 0.0;
379:   PetscDrawLGAddPoint(lg,&x,&y);
380:   if (its <= 20 || !(its % 5) || eps->reason) {
381:     PetscDrawLGDraw(lg);
382:     PetscDrawLGSave(lg);
383:   }
384:   return(0);
385: }

387: PetscErrorCode EPSMonitorLGAll(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
388: {
389:   PetscDrawLG    lg = (PetscDrawLG)ctx;
390:   PetscInt       i,n = PetscMin(eps->nev,255);
391:   PetscReal      *x,*y;

396:   if (its==1) {
397:     PetscDrawLGReset(lg);
398:     PetscDrawLGSetDimension(lg,n);
399:     PetscDrawLGSetLimits(lg,1,1.0,PetscLog10Real(eps->tol)-2,0.0);
400:   }
401:   PetscMalloc2(n,&x,n,&y);
402:   for (i=0;i<n;i++) {
403:     x[i] = (PetscReal)its;
404:     if (i < nest && errest[i] > 0.0) y[i] = PetscLog10Real(errest[i]);
405:     else y[i] = 0.0;
406:   }
407:   PetscDrawLGAddPoint(lg,x,y);
408:   if (its <= 20 || !(its % 5) || eps->reason) {
409:     PetscDrawLGDraw(lg);
410:     PetscDrawLGSave(lg);
411:   }
412:   PetscFree2(x,y);
413:   return(0);
414: }