Actual source code: pepview.c

slepc-3.17.1 2022-04-11
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:    The PEP routines related to various viewers
 12: */

 14: #include <slepc/private/pepimpl.h>
 15: #include <slepc/private/bvimpl.h>
 16: #include <petscdraw.h>

 18: /*@C
 19:    PEPView - Prints the PEP data structure.

 21:    Collective on pep

 23:    Input Parameters:
 24: +  pep - the polynomial eigenproblem solver context
 25: -  viewer - optional visualization context

 27:    Options Database Key:
 28: .  -pep_view -  Calls PEPView() at end of PEPSolve()

 30:    Note:
 31:    The available visualization contexts include
 32: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 33: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 34:          output where only the first processor opens
 35:          the file.  All other processors send their
 36:          data to the first processor to print.

 38:    The user can open an alternative visualization context with
 39:    PetscViewerASCIIOpen() - output to a specified file.

 41:    Level: beginner

 43: .seealso: STView()
 44: @*/
 45: PetscErrorCode PEPView(PEP pep,PetscViewer viewer)
 46: {
 47:   const char     *type=NULL;
 48:   char           str[50];
 49:   PetscBool      isascii,islinear,istrivial;
 50:   PetscInt       i;
 51:   PetscViewer    sviewer;
 52:   MPI_Comm       child;

 55:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);

 59:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 60:   if (isascii) {
 61:     PetscObjectPrintClassNamePrefixType((PetscObject)pep,viewer);
 62:     if (pep->ops->view) {
 63:       PetscViewerASCIIPushTab(viewer);
 64:       (*pep->ops->view)(pep,viewer);
 65:       PetscViewerASCIIPopTab(viewer);
 66:     }
 67:     if (pep->problem_type) {
 68:       switch (pep->problem_type) {
 69:         case PEP_GENERAL:    type = "general polynomial eigenvalue problem"; break;
 70:         case PEP_HERMITIAN:  type = SLEPC_STRING_HERMITIAN " polynomial eigenvalue problem"; break;
 71:         case PEP_HYPERBOLIC: type = "hyperbolic polynomial eigenvalue problem"; break;
 72:         case PEP_GYROSCOPIC: type = "gyroscopic polynomial eigenvalue problem"; break;
 73:       }
 74:     } else type = "not yet set";
 75:     PetscViewerASCIIPrintf(viewer,"  problem type: %s\n",type);
 76:     PetscViewerASCIIPrintf(viewer,"  polynomial represented in %s basis\n",PEPBasisTypes[pep->basis]);
 77:     switch (pep->scale) {
 78:       case PEP_SCALE_NONE:
 79:         break;
 80:       case PEP_SCALE_SCALAR:
 81:         PetscViewerASCIIPrintf(viewer,"  parameter scaling enabled, with scaling factor=%g\n",(double)pep->sfactor);
 82:         break;
 83:       case PEP_SCALE_DIAGONAL:
 84:         PetscViewerASCIIPrintf(viewer,"  diagonal balancing enabled, with its=%" PetscInt_FMT " and lambda=%g\n",pep->sits,(double)pep->slambda);
 85:         break;
 86:       case PEP_SCALE_BOTH:
 87:         PetscViewerASCIIPrintf(viewer,"  parameter scaling & diagonal balancing enabled, with scaling factor=%g, its=%" PetscInt_FMT " and lambda=%g\n",(double)pep->sfactor,pep->sits,(double)pep->slambda);
 88:         break;
 89:     }
 90:     PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: ");
 91:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 92:     SlepcSNPrintfScalar(str,sizeof(str),pep->target,PETSC_FALSE);
 93:     if (!pep->which) PetscViewerASCIIPrintf(viewer,"not yet set\n");
 94:     else switch (pep->which) {
 95:       case PEP_WHICH_USER:
 96:         PetscViewerASCIIPrintf(viewer,"user defined\n");
 97:         break;
 98:       case PEP_TARGET_MAGNITUDE:
 99:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (in magnitude)\n",str);
100:         break;
101:       case PEP_TARGET_REAL:
102:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the real axis)\n",str);
103:         break;
104:       case PEP_TARGET_IMAGINARY:
105:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the imaginary axis)\n",str);
106:         break;
107:       case PEP_LARGEST_MAGNITUDE:
108:         PetscViewerASCIIPrintf(viewer,"largest eigenvalues in magnitude\n");
109:         break;
110:       case PEP_SMALLEST_MAGNITUDE:
111:         PetscViewerASCIIPrintf(viewer,"smallest eigenvalues in magnitude\n");
112:         break;
113:       case PEP_LARGEST_REAL:
114:         PetscViewerASCIIPrintf(viewer,"largest real parts\n");
115:         break;
116:       case PEP_SMALLEST_REAL:
117:         PetscViewerASCIIPrintf(viewer,"smallest real parts\n");
118:         break;
119:       case PEP_LARGEST_IMAGINARY:
120:         PetscViewerASCIIPrintf(viewer,"largest imaginary parts\n");
121:         break;
122:       case PEP_SMALLEST_IMAGINARY:
123:         PetscViewerASCIIPrintf(viewer,"smallest imaginary parts\n");
124:         break;
125:       case PEP_ALL:
126:         if (pep->inta || pep->intb) PetscViewerASCIIPrintf(viewer,"all eigenvalues in interval [%g,%g]\n",(double)pep->inta,(double)pep->intb);
127:         else PetscViewerASCIIPrintf(viewer,"all eigenvalues in the region\n");
128:         break;
129:     }
130:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
131:     PetscViewerASCIIPrintf(viewer,"  number of eigenvalues (nev): %" PetscInt_FMT "\n",pep->nev);
132:     PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",pep->ncv);
133:     PetscViewerASCIIPrintf(viewer,"  maximum dimension of projected problem (mpd): %" PetscInt_FMT "\n",pep->mpd);
134:     PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",pep->max_it);
135:     PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)pep->tol);
136:     PetscViewerASCIIPrintf(viewer,"  convergence test: ");
137:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
138:     switch (pep->conv) {
139:     case PEP_CONV_ABS:
140:       PetscViewerASCIIPrintf(viewer,"absolute\n");break;
141:     case PEP_CONV_REL:
142:       PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue\n");break;
143:     case PEP_CONV_NORM:
144:       PetscViewerASCIIPrintf(viewer,"relative to the matrix norms\n");
145:       if (pep->nrma) {
146:         PetscViewerASCIIPrintf(viewer,"  computed matrix norms: %g",(double)pep->nrma[0]);
147:         for (i=1;i<pep->nmat;i++) PetscViewerASCIIPrintf(viewer,", %g",(double)pep->nrma[i]);
148:         PetscViewerASCIIPrintf(viewer,"\n");
149:       }
150:       break;
151:     case PEP_CONV_USER:
152:       PetscViewerASCIIPrintf(viewer,"user-defined\n");break;
153:     }
154:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
155:     PetscViewerASCIIPrintf(viewer,"  extraction type: %s\n",PEPExtractTypes[pep->extract]);
156:     if (pep->refine) {
157:       PetscViewerASCIIPrintf(viewer,"  iterative refinement: %s, with %s scheme\n",PEPRefineTypes[pep->refine],PEPRefineSchemes[pep->scheme]);
158:       PetscViewerASCIIPrintf(viewer,"  refinement stopping criterion: tol=%g, its=%" PetscInt_FMT "\n",(double)pep->rtol,pep->rits);
159:       if (pep->npart>1) PetscViewerASCIIPrintf(viewer,"  splitting communicator in %" PetscInt_FMT " partitions for refinement\n",pep->npart);
160:     }
161:     if (pep->nini) PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial space: %" PetscInt_FMT "\n",PetscAbs(pep->nini));
162:   } else {
163:     if (pep->ops->view) (*pep->ops->view)(pep,viewer);
164:   }
165:   PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
166:   if (!pep->V) PEPGetBV(pep,&pep->V);
167:   BVView(pep->V,viewer);
168:   if (!pep->rg) PEPGetRG(pep,&pep->rg);
169:   RGIsTrivial(pep->rg,&istrivial);
170:   if (!istrivial) RGView(pep->rg,viewer);
171:   PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear);
172:   if (!islinear) {
173:     if (!pep->ds) PEPGetDS(pep,&pep->ds);
174:     DSView(pep->ds,viewer);
175:   }
176:   PetscViewerPopFormat(viewer);
177:   if (!pep->st) PEPGetST(pep,&pep->st);
178:   STView(pep->st,viewer);
179:   if (pep->refine!=PEP_REFINE_NONE) {
180:     PetscViewerASCIIPushTab(viewer);
181:     if (pep->npart>1) {
182:       if (pep->refinesubc->color==0) {
183:         PetscSubcommGetChild(pep->refinesubc,&child);
184:         PetscViewerASCIIGetStdout(child,&sviewer);
185:         KSPView(pep->refineksp,sviewer);
186:       }
187:     } else KSPView(pep->refineksp,viewer);
188:     PetscViewerASCIIPopTab(viewer);
189:   }
190:   PetscFunctionReturn(0);
191: }

193: /*@C
194:    PEPViewFromOptions - View from options

196:    Collective on PEP

198:    Input Parameters:
199: +  pep  - the eigensolver context
200: .  obj  - optional object
201: -  name - command line option

203:    Level: intermediate

205: .seealso: PEPView(), PEPCreate()
206: @*/
207: PetscErrorCode PEPViewFromOptions(PEP pep,PetscObject obj,const char name[])
208: {
210:   PetscObjectViewFromOptions((PetscObject)pep,obj,name);
211:   PetscFunctionReturn(0);
212: }

214: /*@C
215:    PEPConvergedReasonView - Displays the reason a PEP solve converged or diverged.

217:    Collective on pep

219:    Input Parameters:
220: +  pep - the eigensolver context
221: -  viewer - the viewer to display the reason

223:    Options Database Keys:
224: .  -pep_converged_reason - print reason for convergence, and number of iterations

226:    Note:
227:    To change the format of the output call PetscViewerPushFormat(viewer,format) before
228:    this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
229:    display a reason if it fails. The latter can be set in the command line with
230:    -pep_converged_reason ::failed

232:    Level: intermediate

234: .seealso: PEPSetConvergenceTest(), PEPSetTolerances(), PEPGetIterationNumber(), PEPConvergedReasonViewFromOptions()
235: @*/
236: PetscErrorCode PEPConvergedReasonView(PEP pep,PetscViewer viewer)
237: {
238:   PetscBool         isAscii;
239:   PetscViewerFormat format;

241:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pep));
242:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
243:   if (isAscii) {
244:     PetscViewerGetFormat(viewer,&format);
245:     PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
246:     if (pep->reason > 0 && format != PETSC_VIEWER_FAILED) PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve converged (%" PetscInt_FMT " eigenpair%s) due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",pep->nconv,(pep->nconv>1)?"s":"",PEPConvergedReasons[pep->reason],pep->its);
247:     else if (pep->reason <= 0) PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",PEPConvergedReasons[pep->reason],pep->its);
248:     PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
249:   }
250:   PetscFunctionReturn(0);
251: }

253: /*@
254:    PEPConvergedReasonViewFromOptions - Processes command line options to determine if/how
255:    the PEP converged reason is to be viewed.

257:    Collective on pep

259:    Input Parameter:
260: .  pep - the eigensolver context

262:    Level: developer

264: .seealso: PEPConvergedReasonView()
265: @*/
266: PetscErrorCode PEPConvergedReasonViewFromOptions(PEP pep)
267: {
268:   PetscViewer       viewer;
269:   PetscBool         flg;
270:   static PetscBool  incall = PETSC_FALSE;
271:   PetscViewerFormat format;

273:   if (incall) PetscFunctionReturn(0);
274:   incall = PETSC_TRUE;
275:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_converged_reason",&viewer,&format,&flg);
276:   if (flg) {
277:     PetscViewerPushFormat(viewer,format);
278:     PEPConvergedReasonView(pep,viewer);
279:     PetscViewerPopFormat(viewer);
280:     PetscViewerDestroy(&viewer);
281:   }
282:   incall = PETSC_FALSE;
283:   PetscFunctionReturn(0);
284: }

286: static PetscErrorCode PEPErrorView_ASCII(PEP pep,PEPErrorType etype,PetscViewer viewer)
287: {
288:   PetscReal      error;
289:   PetscInt       i,j,k,nvals;

291:   nvals = (pep->which==PEP_ALL)? pep->nconv: pep->nev;
292:   if (pep->which!=PEP_ALL && pep->nconv<pep->nev) {
293:     PetscViewerASCIIPrintf(viewer," Problem: less than %" PetscInt_FMT " eigenvalues converged\n\n",pep->nev);
294:     PetscFunctionReturn(0);
295:   }
296:   if (pep->which==PEP_ALL && !nvals) {
297:     PetscViewerASCIIPrintf(viewer," No eigenvalues have been found\n\n");
298:     PetscFunctionReturn(0);
299:   }
300:   for (i=0;i<nvals;i++) {
301:     PEPComputeError(pep,i,etype,&error);
302:     if (error>=5.0*pep->tol) {
303:       PetscViewerASCIIPrintf(viewer," Problem: some of the first %" PetscInt_FMT " relative errors are higher than the tolerance\n\n",nvals);
304:       PetscFunctionReturn(0);
305:     }
306:   }
307:   if (pep->which==PEP_ALL) PetscViewerASCIIPrintf(viewer," Found %" PetscInt_FMT " eigenvalues, all of them computed up to the required tolerance:",nvals);
308:   else PetscViewerASCIIPrintf(viewer," All requested eigenvalues computed up to the required tolerance:");
309:   for (i=0;i<=(nvals-1)/8;i++) {
310:     PetscViewerASCIIPrintf(viewer,"\n     ");
311:     for (j=0;j<PetscMin(8,nvals-8*i);j++) {
312:       k = pep->perm[8*i+j];
313:       SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]);
314:       if (8*i+j+1<nvals) PetscViewerASCIIPrintf(viewer,", ");
315:     }
316:   }
317:   PetscViewerASCIIPrintf(viewer,"\n\n");
318:   PetscFunctionReturn(0);
319: }

321: static PetscErrorCode PEPErrorView_DETAIL(PEP pep,PEPErrorType etype,PetscViewer viewer)
322: {
323:   PetscReal      error,re,im;
324:   PetscScalar    kr,ki;
325:   PetscInt       i;
326:   char           ex[30],sep[]=" ---------------------- --------------------\n";

328:   if (!pep->nconv) PetscFunctionReturn(0);
329:   switch (etype) {
330:     case PEP_ERROR_ABSOLUTE:
331:       PetscSNPrintf(ex,sizeof(ex),"   ||P(k)x||");
332:       break;
333:     case PEP_ERROR_RELATIVE:
334:       PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||");
335:       break;
336:     case PEP_ERROR_BACKWARD:
337:       PetscSNPrintf(ex,sizeof(ex),"    eta(x,k)");
338:       break;
339:   }
340:   PetscViewerASCIIPrintf(viewer,"%s            k             %s\n%s",sep,ex,sep);
341:   for (i=0;i<pep->nconv;i++) {
342:     PEPGetEigenpair(pep,i,&kr,&ki,NULL,NULL);
343:     PEPComputeError(pep,i,etype,&error);
344: #if defined(PETSC_USE_COMPLEX)
345:     re = PetscRealPart(kr);
346:     im = PetscImaginaryPart(kr);
347: #else
348:     re = kr;
349:     im = ki;
350: #endif
351:     if (im!=0.0) PetscViewerASCIIPrintf(viewer,"  % 9f%+9fi      %12g\n",(double)re,(double)im,(double)error);
352:     else PetscViewerASCIIPrintf(viewer,"    % 12f           %12g\n",(double)re,(double)error);
353:   }
354:   PetscViewerASCIIPrintf(viewer,"%s",sep);
355:   PetscFunctionReturn(0);
356: }

358: static PetscErrorCode PEPErrorView_MATLAB(PEP pep,PEPErrorType etype,PetscViewer viewer)
359: {
360:   PetscReal      error;
361:   PetscInt       i;
362:   const char     *name;

364:   PetscObjectGetName((PetscObject)pep,&name);
365:   PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name);
366:   for (i=0;i<pep->nconv;i++) {
367:     PEPComputeError(pep,i,etype,&error);
368:     PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error);
369:   }
370:   PetscViewerASCIIPrintf(viewer,"];\n");
371:   PetscFunctionReturn(0);
372: }

374: /*@C
375:    PEPErrorView - Displays the errors associated with the computed solution
376:    (as well as the eigenvalues).

378:    Collective on pep

380:    Input Parameters:
381: +  pep    - the eigensolver context
382: .  etype  - error type
383: -  viewer - optional visualization context

385:    Options Database Keys:
386: +  -pep_error_absolute - print absolute errors of each eigenpair
387: .  -pep_error_relative - print relative errors of each eigenpair
388: -  -pep_error_backward - print backward errors of each eigenpair

390:    Notes:
391:    By default, this function checks the error of all eigenpairs and prints
392:    the eigenvalues if all of them are below the requested tolerance.
393:    If the viewer has format=PETSC_VIEWER_ASCII_INFO_DETAIL then a table with
394:    eigenvalues and corresponding errors is printed.

396:    Level: intermediate

398: .seealso: PEPSolve(), PEPValuesView(), PEPVectorsView()
399: @*/
400: PetscErrorCode PEPErrorView(PEP pep,PEPErrorType etype,PetscViewer viewer)
401: {
402:   PetscBool         isascii;
403:   PetscViewerFormat format;

406:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
409:   PEPCheckSolved(pep,1);
410:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
411:   if (!isascii) PetscFunctionReturn(0);

413:   PetscViewerGetFormat(viewer,&format);
414:   switch (format) {
415:     case PETSC_VIEWER_DEFAULT:
416:     case PETSC_VIEWER_ASCII_INFO:
417:       PEPErrorView_ASCII(pep,etype,viewer);
418:       break;
419:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
420:       PEPErrorView_DETAIL(pep,etype,viewer);
421:       break;
422:     case PETSC_VIEWER_ASCII_MATLAB:
423:       PEPErrorView_MATLAB(pep,etype,viewer);
424:       break;
425:     default:
426:       PetscInfo(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
427:   }
428:   PetscFunctionReturn(0);
429: }

431: /*@
432:    PEPErrorViewFromOptions - Processes command line options to determine if/how
433:    the errors of the computed solution are to be viewed.

435:    Collective on pep

437:    Input Parameter:
438: .  pep - the eigensolver context

440:    Level: developer

442: .seealso: PEPErrorView()
443: @*/
444: PetscErrorCode PEPErrorViewFromOptions(PEP pep)
445: {
446:   PetscViewer       viewer;
447:   PetscBool         flg;
448:   static PetscBool  incall = PETSC_FALSE;
449:   PetscViewerFormat format;

451:   if (incall) PetscFunctionReturn(0);
452:   incall = PETSC_TRUE;
453:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_absolute",&viewer,&format,&flg);
454:   if (flg) {
455:     PetscViewerPushFormat(viewer,format);
456:     PEPErrorView(pep,PEP_ERROR_ABSOLUTE,viewer);
457:     PetscViewerPopFormat(viewer);
458:     PetscViewerDestroy(&viewer);
459:   }
460:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_relative",&viewer,&format,&flg);
461:   if (flg) {
462:     PetscViewerPushFormat(viewer,format);
463:     PEPErrorView(pep,PEP_ERROR_RELATIVE,viewer);
464:     PetscViewerPopFormat(viewer);
465:     PetscViewerDestroy(&viewer);
466:   }
467:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_backward",&viewer,&format,&flg);
468:   if (flg) {
469:     PetscViewerPushFormat(viewer,format);
470:     PEPErrorView(pep,PEP_ERROR_BACKWARD,viewer);
471:     PetscViewerPopFormat(viewer);
472:     PetscViewerDestroy(&viewer);
473:   }
474:   incall = PETSC_FALSE;
475:   PetscFunctionReturn(0);
476: }

478: static PetscErrorCode PEPValuesView_DRAW(PEP pep,PetscViewer viewer)
479: {
480:   PetscDraw      draw;
481:   PetscDrawSP    drawsp;
482:   PetscReal      re,im;
483:   PetscInt       i,k;

485:   if (!pep->nconv) PetscFunctionReturn(0);
486:   PetscViewerDrawGetDraw(viewer,0,&draw);
487:   PetscDrawSetTitle(draw,"Computed Eigenvalues");
488:   PetscDrawSPCreate(draw,1,&drawsp);
489:   for (i=0;i<pep->nconv;i++) {
490:     k = pep->perm[i];
491: #if defined(PETSC_USE_COMPLEX)
492:     re = PetscRealPart(pep->eigr[k]);
493:     im = PetscImaginaryPart(pep->eigr[k]);
494: #else
495:     re = pep->eigr[k];
496:     im = pep->eigi[k];
497: #endif
498:     PetscDrawSPAddPoint(drawsp,&re,&im);
499:   }
500:   PetscDrawSPDraw(drawsp,PETSC_TRUE);
501:   PetscDrawSPSave(drawsp);
502:   PetscDrawSPDestroy(&drawsp);
503:   PetscFunctionReturn(0);
504: }

506: static PetscErrorCode PEPValuesView_BINARY(PEP pep,PetscViewer viewer)
507: {
508: #if defined(PETSC_HAVE_COMPLEX)
509:   PetscInt       i,k;
510:   PetscComplex   *ev;
511: #endif

513: #if defined(PETSC_HAVE_COMPLEX)
514:   PetscMalloc1(pep->nconv,&ev);
515:   for (i=0;i<pep->nconv;i++) {
516:     k = pep->perm[i];
517: #if defined(PETSC_USE_COMPLEX)
518:     ev[i] = pep->eigr[k];
519: #else
520:     ev[i] = PetscCMPLX(pep->eigr[k],pep->eigi[k]);
521: #endif
522:   }
523:   PetscViewerBinaryWrite(viewer,ev,pep->nconv,PETSC_COMPLEX);
524:   PetscFree(ev);
525: #endif
526:   PetscFunctionReturn(0);
527: }

529: #if defined(PETSC_HAVE_HDF5)
530: static PetscErrorCode PEPValuesView_HDF5(PEP pep,PetscViewer viewer)
531: {
532:   PetscInt       i,k,n,N;
533:   PetscMPIInt    rank;
534:   Vec            v;
535:   char           vname[30];
536:   const char     *ename;

538:   MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rank);
539:   N = pep->nconv;
540:   n = rank? 0: N;
541:   /* create a vector containing the eigenvalues */
542:   VecCreateMPI(PetscObjectComm((PetscObject)pep),n,N,&v);
543:   PetscObjectGetName((PetscObject)pep,&ename);
544:   PetscSNPrintf(vname,sizeof(vname),"eigr_%s",ename);
545:   PetscObjectSetName((PetscObject)v,vname);
546:   if (!rank) {
547:     for (i=0;i<pep->nconv;i++) {
548:       k = pep->perm[i];
549:       VecSetValue(v,i,pep->eigr[k],INSERT_VALUES);
550:     }
551:   }
552:   VecAssemblyBegin(v);
553:   VecAssemblyEnd(v);
554:   VecView(v,viewer);
555: #if !defined(PETSC_USE_COMPLEX)
556:   /* in real scalars write the imaginary part as a separate vector */
557:   PetscSNPrintf(vname,sizeof(vname),"eigi_%s",ename);
558:   PetscObjectSetName((PetscObject)v,vname);
559:   if (!rank) {
560:     for (i=0;i<pep->nconv;i++) {
561:       k = pep->perm[i];
562:       VecSetValue(v,i,pep->eigi[k],INSERT_VALUES);
563:     }
564:   }
565:   VecAssemblyBegin(v);
566:   VecAssemblyEnd(v);
567:   VecView(v,viewer);
568: #endif
569:   VecDestroy(&v);
570:   PetscFunctionReturn(0);
571: }
572: #endif

574: static PetscErrorCode PEPValuesView_ASCII(PEP pep,PetscViewer viewer)
575: {
576:   PetscInt       i,k;

578:   PetscViewerASCIIPrintf(viewer,"Eigenvalues = \n");
579:   for (i=0;i<pep->nconv;i++) {
580:     k = pep->perm[i];
581:     PetscViewerASCIIPrintf(viewer,"   ");
582:     SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]);
583:     PetscViewerASCIIPrintf(viewer,"\n");
584:   }
585:   PetscViewerASCIIPrintf(viewer,"\n");
586:   PetscFunctionReturn(0);
587: }

589: static PetscErrorCode PEPValuesView_MATLAB(PEP pep,PetscViewer viewer)
590: {
591:   PetscInt       i,k;
592:   PetscReal      re,im;
593:   const char     *name;

595:   PetscObjectGetName((PetscObject)pep,&name);
596:   PetscViewerASCIIPrintf(viewer,"Lambda_%s = [\n",name);
597:   for (i=0;i<pep->nconv;i++) {
598:     k = pep->perm[i];
599: #if defined(PETSC_USE_COMPLEX)
600:     re = PetscRealPart(pep->eigr[k]);
601:     im = PetscImaginaryPart(pep->eigr[k]);
602: #else
603:     re = pep->eigr[k];
604:     im = pep->eigi[k];
605: #endif
606:     if (im!=0.0) PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei\n",(double)re,(double)im);
607:     else PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)re);
608:   }
609:   PetscViewerASCIIPrintf(viewer,"];\n");
610:   PetscFunctionReturn(0);
611: }

613: /*@C
614:    PEPValuesView - Displays the computed eigenvalues in a viewer.

616:    Collective on pep

618:    Input Parameters:
619: +  pep    - the eigensolver context
620: -  viewer - the viewer

622:    Options Database Key:
623: .  -pep_view_values - print computed eigenvalues

625:    Level: intermediate

627: .seealso: PEPSolve(), PEPVectorsView(), PEPErrorView()
628: @*/
629: PetscErrorCode PEPValuesView(PEP pep,PetscViewer viewer)
630: {
631:   PetscBool         isascii,isdraw,isbinary;
632:   PetscViewerFormat format;
633: #if defined(PETSC_HAVE_HDF5)
634:   PetscBool         ishdf5;
635: #endif

638:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
641:   PEPCheckSolved(pep,1);
642:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
643:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
644: #if defined(PETSC_HAVE_HDF5)
645:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
646: #endif
647:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
648:   if (isdraw) PEPValuesView_DRAW(pep,viewer);
649:   else if (isbinary) PEPValuesView_BINARY(pep,viewer);
650: #if defined(PETSC_HAVE_HDF5)
651:   else if (ishdf5) PEPValuesView_HDF5(pep,viewer);
652: #endif
653:   else if (isascii) {
654:     PetscViewerGetFormat(viewer,&format);
655:     switch (format) {
656:       case PETSC_VIEWER_DEFAULT:
657:       case PETSC_VIEWER_ASCII_INFO:
658:       case PETSC_VIEWER_ASCII_INFO_DETAIL:
659:         PEPValuesView_ASCII(pep,viewer);
660:         break;
661:       case PETSC_VIEWER_ASCII_MATLAB:
662:         PEPValuesView_MATLAB(pep,viewer);
663:         break;
664:       default:
665:         PetscInfo(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
666:     }
667:   }
668:   PetscFunctionReturn(0);
669: }

671: /*@
672:    PEPValuesViewFromOptions - Processes command line options to determine if/how
673:    the computed eigenvalues are to be viewed.

675:    Collective on pep

677:    Input Parameter:
678: .  pep - the eigensolver context

680:    Level: developer

682: .seealso: PEPValuesView()
683: @*/
684: PetscErrorCode PEPValuesViewFromOptions(PEP pep)
685: {
686:   PetscViewer       viewer;
687:   PetscBool         flg;
688:   static PetscBool  incall = PETSC_FALSE;
689:   PetscViewerFormat format;

691:   if (incall) PetscFunctionReturn(0);
692:   incall = PETSC_TRUE;
693:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_values",&viewer,&format,&flg);
694:   if (flg) {
695:     PetscViewerPushFormat(viewer,format);
696:     PEPValuesView(pep,viewer);
697:     PetscViewerPopFormat(viewer);
698:     PetscViewerDestroy(&viewer);
699:   }
700:   incall = PETSC_FALSE;
701:   PetscFunctionReturn(0);
702: }

704: /*@C
705:    PEPVectorsView - Outputs computed eigenvectors to a viewer.

707:    Collective on pep

709:    Input Parameters:
710: +  pep    - the eigensolver context
711: -  viewer - the viewer

713:    Options Database Key:
714: .  -pep_view_vectors - output eigenvectors.

716:    Note:
717:    If PETSc was configured with real scalars, complex conjugate eigenvectors
718:    will be viewed as two separate real vectors, one containing the real part
719:    and another one containing the imaginary part.

721:    Level: intermediate

723: .seealso: PEPSolve(), PEPValuesView(), PEPErrorView()
724: @*/
725: PetscErrorCode PEPVectorsView(PEP pep,PetscViewer viewer)
726: {
727:   PetscInt       i,k;
728:   Vec            xr,xi=NULL;

731:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
734:   PEPCheckSolved(pep,1);
735:   if (pep->nconv) {
736:     PEPComputeVectors(pep);
737:     BVCreateVec(pep->V,&xr);
738: #if !defined(PETSC_USE_COMPLEX)
739:     BVCreateVec(pep->V,&xi);
740: #endif
741:     for (i=0;i<pep->nconv;i++) {
742:       k = pep->perm[i];
743:       BV_GetEigenvector(pep->V,k,pep->eigi[k],xr,xi);
744:       SlepcViewEigenvector(viewer,xr,xi,"X",i,(PetscObject)pep);
745:     }
746:     VecDestroy(&xr);
747: #if !defined(PETSC_USE_COMPLEX)
748:     VecDestroy(&xi);
749: #endif
750:   }
751:   PetscFunctionReturn(0);
752: }

754: /*@
755:    PEPVectorsViewFromOptions - Processes command line options to determine if/how
756:    the computed eigenvectors are to be viewed.

758:    Collective on pep

760:    Input Parameter:
761: .  pep - the eigensolver context

763:    Level: developer

765: .seealso: PEPVectorsView()
766: @*/
767: PetscErrorCode PEPVectorsViewFromOptions(PEP pep)
768: {
769:   PetscViewer       viewer;
770:   PetscBool         flg = PETSC_FALSE;
771:   static PetscBool  incall = PETSC_FALSE;
772:   PetscViewerFormat format;

774:   if (incall) PetscFunctionReturn(0);
775:   incall = PETSC_TRUE;
776:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_vectors",&viewer,&format,&flg);
777:   if (flg) {
778:     PetscViewerPushFormat(viewer,format);
779:     PEPVectorsView(pep,viewer);
780:     PetscViewerPopFormat(viewer);
781:     PetscViewerDestroy(&viewer);
782:   }
783:   incall = PETSC_FALSE;
784:   PetscFunctionReturn(0);
785: }