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

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

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

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

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

 36:    Logically Collective on PEP

 38:    Input Parameters:
 39: +  pep     - eigensolver context obtained from PEPCreate()
 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(PEP pep,int its,int nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal* errest,int nest,void *mctx)

 48: +  pep    - polynomial eigensolver context obtained from PEPCreate()
 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 PEPMonitorSet()

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

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

 75:    Level: intermediate

 77: .seealso: PEPMonitorFirst(), PEPMonitorAll(), PEPMonitorCancel()
 78: @*/
 79: PetscErrorCode PEPMonitorSet(PEP pep,PetscErrorCode (*monitor)(PEP,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
 80: {
 83:   if (pep->numbermonitors >= MAXPEPMONITORS) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Too many PEP monitors set");
 84:   pep->monitor[pep->numbermonitors]           = monitor;
 85:   pep->monitorcontext[pep->numbermonitors]    = (void*)mctx;
 86:   pep->monitordestroy[pep->numbermonitors++]  = monitordestroy;
 87:   return(0);
 88: }

 90: /*@
 91:    PEPMonitorCancel - Clears all monitors for a PEP object.

 93:    Logically Collective on PEP

 95:    Input Parameters:
 96: .  pep - eigensolver context obtained from PEPCreate()

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

103:    Level: intermediate

105: .seealso: PEPMonitorSet()
106: @*/
107: PetscErrorCode PEPMonitorCancel(PEP pep)
108: {
110:   PetscInt       i;

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

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

127:    Not Collective

129:    Input Parameter:
130: .  pep - eigensolver context obtained from PEPCreate()

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

135:    Level: intermediate

137: .seealso: PEPMonitorSet(), PEPDefaultMonitor()
138: @*/
139: PetscErrorCode PEPGetMonitorContext(PEP pep,void **ctx)
140: {
143:   *ctx = pep->monitorcontext[0];
144:   return(0);
145: }

147: /*
148:    Helper function to compute eigenvalue that must be viewed in monitor
149:  */
150: static PetscErrorCode PEPMonitorGetTrueEig(PEP pep,PetscScalar *er,PetscScalar *ei)
151: {
153:   PetscBool      flg,sinv;

156:   STBackTransform(pep->st,1,er,ei);
157:   if (pep->sfactor!=1.0) {
158:     STGetTransform(pep->st,&flg);
159:     PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
160:     if (flg && sinv) {
161:       *er /= pep->sfactor; *ei /= pep->sfactor;
162:     } else {
163:       *er *= pep->sfactor; *ei *= pep->sfactor;
164:     }
165:   }
166:   return(0);
167: }

169: /*@C
170:    PEPMonitorAll - Print the current approximate values and
171:    error estimates at each iteration of the polynomial eigensolver.

173:    Collective on PEP

175:    Input Parameters:
176: +  pep    - polynomial eigensolver context
177: .  its    - iteration number
178: .  nconv  - number of converged eigenpairs so far
179: .  eigr   - real part of the eigenvalues
180: .  eigi   - imaginary part of the eigenvalues
181: .  errest - error estimates
182: .  nest   - number of error estimates to display
183: -  vf     - viewer and format for monitoring

185:    Level: intermediate

187: .seealso: PEPMonitorSet(), PEPMonitorFirst(), PEPMonitorConverged()
188: @*/
189: PetscErrorCode PEPMonitorAll(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,PetscViewerAndFormat *vf)
190: {
192:   PetscInt       i;
193:   PetscScalar    er,ei;
194:   PetscViewer    viewer;

199:   viewer = vf->viewer;
201:   PetscViewerPushFormat(viewer,vf->format);
202:   PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
203:   if (its==1 && ((PetscObject)pep)->prefix) {
204:     PetscViewerASCIIPrintf(viewer,"  Eigenvalue approximations and residual norms for %s solve.\n",((PetscObject)pep)->prefix);
205:   }
206:   PetscViewerASCIIPrintf(viewer,"%3D PEP nconv=%D Values (Errors)",its,nconv);
207:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
208:   for (i=0;i<nest;i++) {
209:     er = eigr[i]; ei = eigi[i];
210:     PEPMonitorGetTrueEig(pep,&er,&ei);
211: #if defined(PETSC_USE_COMPLEX)
212:     PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
213: #else
214:     PetscViewerASCIIPrintf(viewer," %g",(double)er);
215:     if (eigi[i]!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
216: #endif
217:     PetscViewerASCIIPrintf(viewer," (%10.8e)",(double)errest[i]);
218:   }
219:   PetscViewerASCIIPrintf(viewer,"\n");
220:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
221:   PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
222:   PetscViewerPopFormat(viewer);
223:   return(0);
224: }

226: /*@C
227:    PEPMonitorFirst - Print the first unconverged approximate value and
228:    error estimate at each iteration of the polynomial eigensolver.

230:    Collective on PEP

232:    Input Parameters:
233: +  pep    - polynomial eigensolver context
234: .  its    - iteration number
235: .  nconv  - number of converged eigenpairs so far
236: .  eigr   - real part of the eigenvalues
237: .  eigi   - imaginary part of the eigenvalues
238: .  errest - error estimates
239: .  nest   - number of error estimates to display
240: -  vf     - viewer and format for monitoring

242:    Level: intermediate

244: .seealso: PEPMonitorSet(), PEPMonitorAll(), PEPMonitorConverged()
245: @*/
246: PetscErrorCode PEPMonitorFirst(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,PetscViewerAndFormat *vf)
247: {
249:   PetscScalar    er,ei;
250:   PetscViewer    viewer;

255:   viewer = vf->viewer;
257:   if (its==1 && ((PetscObject)pep)->prefix) {
258:     PetscViewerASCIIPrintf(viewer,"  Eigenvalue approximations and residual norms for %s solve.\n",((PetscObject)pep)->prefix);
259:   }
260:   if (nconv<nest) {
261:     PetscViewerPushFormat(viewer,vf->format);
262:     PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
263:     PetscViewerASCIIPrintf(viewer,"%3D PEP nconv=%D first unconverged value (error)",its,nconv);
264:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
265:     er = eigr[nconv]; ei = eigi[nconv];
266:     PEPMonitorGetTrueEig(pep,&er,&ei);
267: #if defined(PETSC_USE_COMPLEX)
268:     PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
269: #else
270:     PetscViewerASCIIPrintf(viewer," %g",(double)er);
271:     if (eigi[nconv]!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
272: #endif
273:     PetscViewerASCIIPrintf(viewer," (%10.8e)\n",(double)errest[nconv]);
274:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
275:     PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
276:     PetscViewerPopFormat(viewer);
277:   }
278:   return(0);
279: }

281: /*@C
282:    PEPMonitorConverged - Print the approximate values and
283:    error estimates as they converge.

285:    Collective on PEP

287:    Input Parameters:
288: +  pep    - polynomial eigensolver context
289: .  its    - iteration number
290: .  nconv  - number of converged eigenpairs so far
291: .  eigr   - real part of the eigenvalues
292: .  eigi   - imaginary part of the eigenvalues
293: .  errest - error estimates
294: .  nest   - number of error estimates to display
295: -  ctx    - monitor context

297:    Level: intermediate

299: .seealso: PEPMonitorSet(), PEPMonitorFirst(), PEPMonitorAll()
300: @*/
301: PetscErrorCode PEPMonitorConverged(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,SlepcConvMonitor ctx)
302: {
304:   PetscInt       i;
305:   PetscScalar    er,ei;
306:   PetscViewer    viewer;

311:   viewer = ctx->viewer;
313:   if (its==1 && ((PetscObject)pep)->prefix) {
314:     PetscViewerASCIIPrintf(viewer,"  Convergence history for %s solve.\n",((PetscObject)pep)->prefix);
315:   }
316:   if (its==1) ctx->oldnconv = 0;
317:   if (ctx->oldnconv!=nconv) {
318:     PetscViewerPushFormat(viewer,ctx->format);
319:     PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
320:     for (i=ctx->oldnconv;i<nconv;i++) {
321:       PetscViewerASCIIPrintf(viewer,"%3D PEP converged value (error) #%D",its,i);
322:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
323:       er = eigr[i]; ei = eigi[i];
324:       PEPMonitorGetTrueEig(pep,&er,&ei);
325: #if defined(PETSC_USE_COMPLEX)
326:       PetscViewerASCIIPrintf(viewer," %g%+gi",(double)PetscRealPart(er),(double)PetscImaginaryPart(er));
327: #else
328:       PetscViewerASCIIPrintf(viewer," %g",(double)er);
329:       if (eigi[i]!=0.0) { PetscViewerASCIIPrintf(viewer,"%+gi",(double)ei); }
330: #endif
331:       PetscViewerASCIIPrintf(viewer," (%10.8e)\n",(double)errest[i]);
332:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
333:     }
334:     PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
335:     PetscViewerPopFormat(viewer);
336:     ctx->oldnconv = nconv;
337:   }
338:   return(0);
339: }

341: /*@C
342:    PEPMonitorLGCreate - Creates a line graph context for use with
343:    PEP to monitor convergence.

345:    Collective on MPI_Comm

347:    Input Parameters:
348: +  comm - communicator context
349: .  host - the X display to open, or null for the local machine
350: .  label - the title to put in the title bar
351: .  x, y - the screen coordinates of the upper left coordinate of
352:           the window
353: -  m, n - the screen width and height in pixels

355:    Output Parameter:
356: .  lgctx - the drawing context

358:    Options Database Keys:
359: +  -pep_monitor_lg - Sets line graph monitor for the first residual
360: -  -pep_monitor_lg_all - Sets line graph monitor for all residuals

362:    Notes:
363:    Use PetscDrawLGDestroy() to destroy this line graph.

365:    Level: intermediate

367: .seealso: PEPMonitorSet()
368: @*/
369: PetscErrorCode PEPMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx)
370: {
371:   PetscDraw      draw;
372:   PetscDrawLG    lg;

376:   PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
377:   PetscDrawSetFromOptions(draw);
378:   PetscDrawLGCreate(draw,1,&lg);
379:   PetscDrawLGSetFromOptions(lg);
380:   PetscDrawDestroy(&draw);
381:   *lgctx = lg;
382:   return(0);
383: }

385: PetscErrorCode PEPMonitorLG(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
386: {
387:   PetscDrawLG    lg = (PetscDrawLG)ctx;
388:   PetscReal      x,y;

393:   if (its==1) {
394:     PetscDrawLGReset(lg);
395:     PetscDrawLGSetDimension(lg,1);
396:     PetscDrawLGSetLimits(lg,1,1.0,PetscLog10Real(pep->tol)-2,0.0);
397:   }
398:   x = (PetscReal)its;
399:   if (errest[nconv] > 0.0) y = PetscLog10Real(errest[nconv]);
400:   else y = 0.0;
401:   PetscDrawLGAddPoint(lg,&x,&y);
402:   if (its <= 20 || !(its % 5) || pep->reason) {
403:     PetscDrawLGDraw(lg);
404:     PetscDrawLGSave(lg);
405:   }
406:   return(0);
407: }

409: PetscErrorCode PEPMonitorLGAll(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
410: {
411:   PetscDrawLG    lg = (PetscDrawLG)ctx;
412:   PetscInt       i,n = PetscMin(pep->nev,255);
413:   PetscReal      *x,*y;

418:   if (its==1) {
419:     PetscDrawLGReset(lg);
420:     PetscDrawLGSetDimension(lg,n);
421:     PetscDrawLGSetLimits(lg,1,1.0,PetscLog10Real(pep->tol)-2,0.0);
422:   }
423:   PetscMalloc2(n,&x,n,&y);
424:   for (i=0;i<n;i++) {
425:     x[i] = (PetscReal)its;
426:     if (i < nest && errest[i] > 0.0) y[i] = PetscLog10Real(errest[i]);
427:     else y[i] = 0.0;
428:   }
429:   PetscDrawLGAddPoint(lg,x,y);
430:   if (its <= 20 || !(its % 5) || pep->reason) {
431:     PetscDrawLGDraw(lg);
432:     PetscDrawLGSave(lg);
433:   }
434:   PetscFree2(x,y);
435:   return(0);
436: }