Actual source code: pepview.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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
 42: @*/
 43: PetscErrorCode PEPView(PEP pep,PetscViewer viewer)
 44: {
 46:   const char     *type=NULL;
 47:   char           str[50];
 48:   PetscBool      isascii,islinear,istrivial;
 49:   PetscInt       i;
 50:   PetscViewer    sviewer;

 54:   if (!viewer) {
 55:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
 56:   }

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

207: /*@C
208:    PEPViewFromOptions - View from options

210:    Collective on PEP

212:    Input Parameters:
213: +  pep  - the eigensolver context
214: .  obj  - optional object
215: -  name - command line option

217:    Level: intermediate

219: .seealso: PEPView(), PEPCreate()
220: @*/
221: PetscErrorCode PEPViewFromOptions(PEP pep,PetscObject obj,const char name[])
222: {

227:   PetscObjectViewFromOptions((PetscObject)pep,obj,name);
228:   return(0);
229: }

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

234:    Collective on pep

236:    Input Parameters:
237: +  pep - the eigensolver context
238: -  viewer - the viewer to display the reason

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

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

249:    Level: intermediate

251: .seealso: PEPSetConvergenceTest(), PEPSetTolerances(), PEPGetIterationNumber(), PEPConvergedReasonViewFromOptions()
252: @*/
253: PetscErrorCode PEPConvergedReasonView(PEP pep,PetscViewer viewer)
254: {
255:   PetscErrorCode    ierr;
256:   PetscBool         isAscii;
257:   PetscViewerFormat format;

260:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pep));
261:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
262:   if (isAscii) {
263:     PetscViewerGetFormat(viewer,&format);
264:     PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
265:     if (pep->reason > 0 && format != PETSC_VIEWER_FAILED) {
266:       PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve converged (%D eigenpair%s) due to %s; iterations %D\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",pep->nconv,(pep->nconv>1)?"s":"",PEPConvergedReasons[pep->reason],pep->its);
267:     } else if (pep->reason <= 0) {
268:       PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve did not converge due to %s; iterations %D\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",PEPConvergedReasons[pep->reason],pep->its);
269:     }
270:     PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
271:   }
272:   return(0);
273: }

275: /*@
276:    PEPConvergedReasonViewFromOptions - Processes command line options to determine if/how
277:    the PEP converged reason is to be viewed.

279:    Collective on pep

281:    Input Parameter:
282: .  pep - the eigensolver context

284:    Level: developer

286: .seealso: PEPConvergedReasonView()
287: @*/
288: PetscErrorCode PEPConvergedReasonViewFromOptions(PEP pep)
289: {
290:   PetscErrorCode    ierr;
291:   PetscViewer       viewer;
292:   PetscBool         flg;
293:   static PetscBool  incall = PETSC_FALSE;
294:   PetscViewerFormat format;

297:   if (incall) return(0);
298:   incall = PETSC_TRUE;
299:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_converged_reason",&viewer,&format,&flg);
300:   if (flg) {
301:     PetscViewerPushFormat(viewer,format);
302:     PEPConvergedReasonView(pep,viewer);
303:     PetscViewerPopFormat(viewer);
304:     PetscViewerDestroy(&viewer);
305:   }
306:   incall = PETSC_FALSE;
307:   return(0);
308: }

310: static PetscErrorCode PEPErrorView_ASCII(PEP pep,PEPErrorType etype,PetscViewer viewer)
311: {
312:   PetscReal      error;
313:   PetscInt       i,j,k,nvals;

317:   nvals = (pep->which==PEP_ALL)? pep->nconv: pep->nev;
318:   if (pep->which!=PEP_ALL && pep->nconv<pep->nev) {
319:     PetscViewerASCIIPrintf(viewer," Problem: less than %D eigenvalues converged\n\n",pep->nev);
320:     return(0);
321:   }
322:   if (pep->which==PEP_ALL && !nvals) {
323:     PetscViewerASCIIPrintf(viewer," No eigenvalues have been found\n\n");
324:     return(0);
325:   }
326:   for (i=0;i<nvals;i++) {
327:     PEPComputeError(pep,i,etype,&error);
328:     if (error>=5.0*pep->tol) {
329:       PetscViewerASCIIPrintf(viewer," Problem: some of the first %D relative errors are higher than the tolerance\n\n",nvals);
330:       return(0);
331:     }
332:   }
333:   if (pep->which==PEP_ALL) {
334:     PetscViewerASCIIPrintf(viewer," Found %D eigenvalues, all of them computed up to the required tolerance:",nvals);
335:   } else {
336:     PetscViewerASCIIPrintf(viewer," All requested eigenvalues computed up to the required tolerance:");
337:   }
338:   for (i=0;i<=(nvals-1)/8;i++) {
339:     PetscViewerASCIIPrintf(viewer,"\n     ");
340:     for (j=0;j<PetscMin(8,nvals-8*i);j++) {
341:       k = pep->perm[8*i+j];
342:       SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]);
343:       if (8*i+j+1<nvals) { PetscViewerASCIIPrintf(viewer,", "); }
344:     }
345:   }
346:   PetscViewerASCIIPrintf(viewer,"\n\n");
347:   return(0);
348: }

350: static PetscErrorCode PEPErrorView_DETAIL(PEP pep,PEPErrorType etype,PetscViewer viewer)
351: {
353:   PetscReal      error,re,im;
354:   PetscScalar    kr,ki;
355:   PetscInt       i;
356:   char           ex[30],sep[]=" ---------------------- --------------------\n";

359:   if (!pep->nconv) return(0);
360:   switch (etype) {
361:     case PEP_ERROR_ABSOLUTE:
362:       PetscSNPrintf(ex,sizeof(ex),"   ||P(k)x||");
363:       break;
364:     case PEP_ERROR_RELATIVE:
365:       PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||");
366:       break;
367:     case PEP_ERROR_BACKWARD:
368:       PetscSNPrintf(ex,sizeof(ex),"    eta(x,k)");
369:       break;
370:   }
371:   PetscViewerASCIIPrintf(viewer,"%s            k             %s\n%s",sep,ex,sep);
372:   for (i=0;i<pep->nconv;i++) {
373:     PEPGetEigenpair(pep,i,&kr,&ki,NULL,NULL);
374:     PEPComputeError(pep,i,etype,&error);
375: #if defined(PETSC_USE_COMPLEX)
376:     re = PetscRealPart(kr);
377:     im = PetscImaginaryPart(kr);
378: #else
379:     re = kr;
380:     im = ki;
381: #endif
382:     if (im!=0.0) {
383:       PetscViewerASCIIPrintf(viewer,"  % 9f%+9fi      %12g\n",(double)re,(double)im,(double)error);
384:     } else {
385:       PetscViewerASCIIPrintf(viewer,"    % 12f           %12g\n",(double)re,(double)error);
386:     }
387:   }
388:   PetscViewerASCIIPrintf(viewer,"%s",sep);
389:   return(0);
390: }

392: static PetscErrorCode PEPErrorView_MATLAB(PEP pep,PEPErrorType etype,PetscViewer viewer)
393: {
395:   PetscReal      error;
396:   PetscInt       i;
397:   const char     *name;

400:   PetscObjectGetName((PetscObject)pep,&name);
401:   PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name);
402:   for (i=0;i<pep->nconv;i++) {
403:     PEPComputeError(pep,i,etype,&error);
404:     PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error);
405:   }
406:   PetscViewerASCIIPrintf(viewer,"];\n");
407:   return(0);
408: }

410: /*@C
411:    PEPErrorView - Displays the errors associated with the computed solution
412:    (as well as the eigenvalues).

414:    Collective on pep

416:    Input Parameters:
417: +  pep    - the eigensolver context
418: .  etype  - error type
419: -  viewer - optional visualization context

421:    Options Database Key:
422: +  -pep_error_absolute - print absolute errors of each eigenpair
423: .  -pep_error_relative - print relative errors of each eigenpair
424: -  -pep_error_backward - print backward errors of each eigenpair

426:    Notes:
427:    By default, this function checks the error of all eigenpairs and prints
428:    the eigenvalues if all of them are below the requested tolerance.
429:    If the viewer has format=PETSC_VIEWER_ASCII_INFO_DETAIL then a table with
430:    eigenvalues and corresponding errors is printed.

432:    Level: intermediate

434: .seealso: PEPSolve(), PEPValuesView(), PEPVectorsView()
435: @*/
436: PetscErrorCode PEPErrorView(PEP pep,PEPErrorType etype,PetscViewer viewer)
437: {
438:   PetscBool         isascii;
439:   PetscViewerFormat format;
440:   PetscErrorCode    ierr;

444:   if (!viewer) {
445:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
446:   }
449:   PEPCheckSolved(pep,1);
450:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
451:   if (!isascii) return(0);

453:   PetscViewerGetFormat(viewer,&format);
454:   switch (format) {
455:     case PETSC_VIEWER_DEFAULT:
456:     case PETSC_VIEWER_ASCII_INFO:
457:       PEPErrorView_ASCII(pep,etype,viewer);
458:       break;
459:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
460:       PEPErrorView_DETAIL(pep,etype,viewer);
461:       break;
462:     case PETSC_VIEWER_ASCII_MATLAB:
463:       PEPErrorView_MATLAB(pep,etype,viewer);
464:       break;
465:     default:
466:       PetscInfo1(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
467:   }
468:   return(0);
469: }

471: /*@
472:    PEPErrorViewFromOptions - Processes command line options to determine if/how
473:    the errors of the computed solution are to be viewed.

475:    Collective on pep

477:    Input Parameter:
478: .  pep - the eigensolver context

480:    Level: developer
481: @*/
482: PetscErrorCode PEPErrorViewFromOptions(PEP pep)
483: {
484:   PetscErrorCode    ierr;
485:   PetscViewer       viewer;
486:   PetscBool         flg;
487:   static PetscBool  incall = PETSC_FALSE;
488:   PetscViewerFormat format;

491:   if (incall) return(0);
492:   incall = PETSC_TRUE;
493:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_absolute",&viewer,&format,&flg);
494:   if (flg) {
495:     PetscViewerPushFormat(viewer,format);
496:     PEPErrorView(pep,PEP_ERROR_ABSOLUTE,viewer);
497:     PetscViewerPopFormat(viewer);
498:     PetscViewerDestroy(&viewer);
499:   }
500:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_relative",&viewer,&format,&flg);
501:   if (flg) {
502:     PetscViewerPushFormat(viewer,format);
503:     PEPErrorView(pep,PEP_ERROR_RELATIVE,viewer);
504:     PetscViewerPopFormat(viewer);
505:     PetscViewerDestroy(&viewer);
506:   }
507:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_backward",&viewer,&format,&flg);
508:   if (flg) {
509:     PetscViewerPushFormat(viewer,format);
510:     PEPErrorView(pep,PEP_ERROR_BACKWARD,viewer);
511:     PetscViewerPopFormat(viewer);
512:     PetscViewerDestroy(&viewer);
513:   }
514:   incall = PETSC_FALSE;
515:   return(0);
516: }

518: static PetscErrorCode PEPValuesView_DRAW(PEP pep,PetscViewer viewer)
519: {
521:   PetscDraw      draw;
522:   PetscDrawSP    drawsp;
523:   PetscReal      re,im;
524:   PetscInt       i,k;

527:   if (!pep->nconv) return(0);
528:   PetscViewerDrawGetDraw(viewer,0,&draw);
529:   PetscDrawSetTitle(draw,"Computed Eigenvalues");
530:   PetscDrawSPCreate(draw,1,&drawsp);
531:   for (i=0;i<pep->nconv;i++) {
532:     k = pep->perm[i];
533: #if defined(PETSC_USE_COMPLEX)
534:     re = PetscRealPart(pep->eigr[k]);
535:     im = PetscImaginaryPart(pep->eigr[k]);
536: #else
537:     re = pep->eigr[k];
538:     im = pep->eigi[k];
539: #endif
540:     PetscDrawSPAddPoint(drawsp,&re,&im);
541:   }
542:   PetscDrawSPDraw(drawsp,PETSC_TRUE);
543:   PetscDrawSPSave(drawsp);
544:   PetscDrawSPDestroy(&drawsp);
545:   return(0);
546: }

548: static PetscErrorCode PEPValuesView_BINARY(PEP pep,PetscViewer viewer)
549: {
550: #if defined(PETSC_HAVE_COMPLEX)
551:   PetscInt       i,k;
552:   PetscComplex   *ev;
554: #endif

557: #if defined(PETSC_HAVE_COMPLEX)
558:   PetscMalloc1(pep->nconv,&ev);
559:   for (i=0;i<pep->nconv;i++) {
560:     k = pep->perm[i];
561: #if defined(PETSC_USE_COMPLEX)
562:     ev[i] = pep->eigr[k];
563: #else
564:     ev[i] = PetscCMPLX(pep->eigr[k],pep->eigi[k]);
565: #endif
566:   }
567:   PetscViewerBinaryWrite(viewer,ev,pep->nconv,PETSC_COMPLEX);
568:   PetscFree(ev);
569: #endif
570:   return(0);
571: }

573: #if defined(PETSC_HAVE_HDF5)
574: static PetscErrorCode PEPValuesView_HDF5(PEP pep,PetscViewer viewer)
575: {
577:   PetscInt       i,k,n,N;
578:   PetscMPIInt    rank;
579:   Vec            v;
580:   char           vname[30];
581:   const char     *ename;

584:   MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rank);
585:   N = pep->nconv;
586:   n = rank? 0: N;
587:   /* create a vector containing the eigenvalues */
588:   VecCreateMPI(PetscObjectComm((PetscObject)pep),n,N,&v);
589:   PetscObjectGetName((PetscObject)pep,&ename);
590:   PetscSNPrintf(vname,sizeof(vname),"eigr_%s",ename);
591:   PetscObjectSetName((PetscObject)v,vname);
592:   if (!rank) {
593:     for (i=0;i<pep->nconv;i++) {
594:       k = pep->perm[i];
595:       VecSetValue(v,i,pep->eigr[k],INSERT_VALUES);
596:     }
597:   }
598:   VecAssemblyBegin(v);
599:   VecAssemblyEnd(v);
600:   VecView(v,viewer);
601: #if !defined(PETSC_USE_COMPLEX)
602:   /* in real scalars write the imaginary part as a separate vector */
603:   PetscSNPrintf(vname,sizeof(vname),"eigi_%s",ename);
604:   PetscObjectSetName((PetscObject)v,vname);
605:   if (!rank) {
606:     for (i=0;i<pep->nconv;i++) {
607:       k = pep->perm[i];
608:       VecSetValue(v,i,pep->eigi[k],INSERT_VALUES);
609:     }
610:   }
611:   VecAssemblyBegin(v);
612:   VecAssemblyEnd(v);
613:   VecView(v,viewer);
614: #endif
615:   VecDestroy(&v);
616:   return(0);
617: }
618: #endif

620: static PetscErrorCode PEPValuesView_ASCII(PEP pep,PetscViewer viewer)
621: {
622:   PetscInt       i,k;

626:   PetscViewerASCIIPrintf(viewer,"Eigenvalues = \n");
627:   for (i=0;i<pep->nconv;i++) {
628:     k = pep->perm[i];
629:     PetscViewerASCIIPrintf(viewer,"   ");
630:     SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]);
631:     PetscViewerASCIIPrintf(viewer,"\n");
632:   }
633:   PetscViewerASCIIPrintf(viewer,"\n");
634:   return(0);
635: }

637: static PetscErrorCode PEPValuesView_MATLAB(PEP pep,PetscViewer viewer)
638: {
640:   PetscInt       i,k;
641:   PetscReal      re,im;
642:   const char     *name;

645:   PetscObjectGetName((PetscObject)pep,&name);
646:   PetscViewerASCIIPrintf(viewer,"Lambda_%s = [\n",name);
647:   for (i=0;i<pep->nconv;i++) {
648:     k = pep->perm[i];
649: #if defined(PETSC_USE_COMPLEX)
650:     re = PetscRealPart(pep->eigr[k]);
651:     im = PetscImaginaryPart(pep->eigr[k]);
652: #else
653:     re = pep->eigr[k];
654:     im = pep->eigi[k];
655: #endif
656:     if (im!=0.0) {
657:       PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei\n",(double)re,(double)im);
658:     } else {
659:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)re);
660:     }
661:   }
662:   PetscViewerASCIIPrintf(viewer,"];\n");
663:   return(0);
664: }

666: /*@C
667:    PEPValuesView - Displays the computed eigenvalues in a viewer.

669:    Collective on pep

671:    Input Parameters:
672: +  pep    - the eigensolver context
673: -  viewer - the viewer

675:    Options Database Key:
676: .  -pep_view_values - print computed eigenvalues

678:    Level: intermediate

680: .seealso: PEPSolve(), PEPVectorsView(), PEPErrorView()
681: @*/
682: PetscErrorCode PEPValuesView(PEP pep,PetscViewer viewer)
683: {
684:   PetscBool         isascii,isdraw,isbinary;
685:   PetscViewerFormat format;
686:   PetscErrorCode    ierr;
687: #if defined(PETSC_HAVE_HDF5)
688:   PetscBool         ishdf5;
689: #endif

693:   if (!viewer) {
694:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
695:   }
698:   PEPCheckSolved(pep,1);
699:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
700:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
701: #if defined(PETSC_HAVE_HDF5)
702:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
703: #endif
704:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
705:   if (isdraw) {
706:     PEPValuesView_DRAW(pep,viewer);
707:   } else if (isbinary) {
708:     PEPValuesView_BINARY(pep,viewer);
709: #if defined(PETSC_HAVE_HDF5)
710:   } else if (ishdf5) {
711:     PEPValuesView_HDF5(pep,viewer);
712: #endif
713:   } else if (isascii) {
714:     PetscViewerGetFormat(viewer,&format);
715:     switch (format) {
716:       case PETSC_VIEWER_DEFAULT:
717:       case PETSC_VIEWER_ASCII_INFO:
718:       case PETSC_VIEWER_ASCII_INFO_DETAIL:
719:         PEPValuesView_ASCII(pep,viewer);
720:         break;
721:       case PETSC_VIEWER_ASCII_MATLAB:
722:         PEPValuesView_MATLAB(pep,viewer);
723:         break;
724:       default:
725:         PetscInfo1(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
726:     }
727:   }
728:   return(0);
729: }

731: /*@
732:    PEPValuesViewFromOptions - Processes command line options to determine if/how
733:    the computed eigenvalues are to be viewed.

735:    Collective on pep

737:    Input Parameter:
738: .  pep - the eigensolver context

740:    Level: developer
741: @*/
742: PetscErrorCode PEPValuesViewFromOptions(PEP pep)
743: {
744:   PetscErrorCode    ierr;
745:   PetscViewer       viewer;
746:   PetscBool         flg;
747:   static PetscBool  incall = PETSC_FALSE;
748:   PetscViewerFormat format;

751:   if (incall) return(0);
752:   incall = PETSC_TRUE;
753:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_values",&viewer,&format,&flg);
754:   if (flg) {
755:     PetscViewerPushFormat(viewer,format);
756:     PEPValuesView(pep,viewer);
757:     PetscViewerPopFormat(viewer);
758:     PetscViewerDestroy(&viewer);
759:   }
760:   incall = PETSC_FALSE;
761:   return(0);
762: }

764: /*@C
765:    PEPVectorsView - Outputs computed eigenvectors to a viewer.

767:    Collective on pep

769:    Input Parameters:
770: +  pep    - the eigensolver context
771: -  viewer - the viewer

773:    Options Database Keys:
774: .  -pep_view_vectors - output eigenvectors.

776:    Note:
777:    If PETSc was configured with real scalars, complex conjugate eigenvectors
778:    will be viewed as two separate real vectors, one containing the real part
779:    and another one containing the imaginary part.

781:    Level: intermediate

783: .seealso: PEPSolve(), PEPValuesView(), PEPErrorView()
784: @*/
785: PetscErrorCode PEPVectorsView(PEP pep,PetscViewer viewer)
786: {
788:   PetscInt       i,k;
789:   Vec            xr,xi=NULL;

793:   if (!viewer) {
794:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
795:   }
798:   PEPCheckSolved(pep,1);
799:   if (pep->nconv) {
800:     PEPComputeVectors(pep);
801:     BVCreateVec(pep->V,&xr);
802: #if !defined(PETSC_USE_COMPLEX)
803:     BVCreateVec(pep->V,&xi);
804: #endif
805:     for (i=0;i<pep->nconv;i++) {
806:       k = pep->perm[i];
807:       BV_GetEigenvector(pep->V,k,pep->eigi[k],xr,xi);
808:       SlepcViewEigenvector(viewer,xr,xi,"X",i,(PetscObject)pep);
809:     }
810:     VecDestroy(&xr);
811: #if !defined(PETSC_USE_COMPLEX)
812:     VecDestroy(&xi);
813: #endif
814:   }
815:   return(0);
816: }

818: /*@
819:    PEPVectorsViewFromOptions - Processes command line options to determine if/how
820:    the computed eigenvectors are to be viewed.

822:    Collective on pep

824:    Input Parameter:
825: .  pep - the eigensolver context

827:    Level: developer
828: @*/
829: PetscErrorCode PEPVectorsViewFromOptions(PEP pep)
830: {
831:   PetscErrorCode    ierr;
832:   PetscViewer       viewer;
833:   PetscBool         flg = PETSC_FALSE;
834:   static PetscBool  incall = PETSC_FALSE;
835:   PetscViewerFormat format;

838:   if (incall) return(0);
839:   incall = PETSC_TRUE;
840:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_vectors",&viewer,&format,&flg);
841:   if (flg) {
842:     PetscViewerPushFormat(viewer,format);
843:     PEPVectorsView(pep,viewer);
844:     PetscViewerPopFormat(viewer);
845:     PetscViewerDestroy(&viewer);
846:   }
847:   incall = PETSC_FALSE;
848:   return(0);
849: }