Actual source code: pepview.c

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

189: /*@C
190:    PEPViewFromOptions - View from options

192:    Collective on PEP

194:    Input Parameters:
195: +  pep  - the eigensolver context
196: .  obj  - optional object
197: -  name - command line option

199:    Level: intermediate

201: .seealso: PEPView(), PEPCreate()
202: @*/
203: PetscErrorCode PEPViewFromOptions(PEP pep,PetscObject obj,const char name[])
204: {
206:   PetscObjectViewFromOptions((PetscObject)pep,obj,name);
207:   return 0;
208: }

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

213:    Collective on pep

215:    Input Parameters:
216: +  pep - the eigensolver context
217: -  viewer - the viewer to display the reason

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

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

228:    Level: intermediate

230: .seealso: PEPSetConvergenceTest(), PEPSetTolerances(), PEPGetIterationNumber(), PEPConvergedReasonViewFromOptions()
231: @*/
232: PetscErrorCode PEPConvergedReasonView(PEP pep,PetscViewer viewer)
233: {
234:   PetscBool         isAscii;
235:   PetscViewerFormat format;

237:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pep));
238:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
239:   if (isAscii) {
240:     PetscViewerGetFormat(viewer,&format);
241:     PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel);
242:     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);
243:     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);
244:     PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel);
245:   }
246:   return 0;
247: }

249: /*@
250:    PEPConvergedReasonViewFromOptions - Processes command line options to determine if/how
251:    the PEP converged reason is to be viewed.

253:    Collective on pep

255:    Input Parameter:
256: .  pep - the eigensolver context

258:    Level: developer

260: .seealso: PEPConvergedReasonView()
261: @*/
262: PetscErrorCode PEPConvergedReasonViewFromOptions(PEP pep)
263: {
264:   PetscViewer       viewer;
265:   PetscBool         flg;
266:   static PetscBool  incall = PETSC_FALSE;
267:   PetscViewerFormat format;

269:   if (incall) return 0;
270:   incall = PETSC_TRUE;
271:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_converged_reason",&viewer,&format,&flg);
272:   if (flg) {
273:     PetscViewerPushFormat(viewer,format);
274:     PEPConvergedReasonView(pep,viewer);
275:     PetscViewerPopFormat(viewer);
276:     PetscViewerDestroy(&viewer);
277:   }
278:   incall = PETSC_FALSE;
279:   return 0;
280: }

282: static PetscErrorCode PEPErrorView_ASCII(PEP pep,PEPErrorType etype,PetscViewer viewer)
283: {
284:   PetscReal      error;
285:   PetscInt       i,j,k,nvals;

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

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

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

354: static PetscErrorCode PEPErrorView_MATLAB(PEP pep,PEPErrorType etype,PetscViewer viewer)
355: {
356:   PetscReal      error;
357:   PetscInt       i;
358:   const char     *name;

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

370: /*@C
371:    PEPErrorView - Displays the errors associated with the computed solution
372:    (as well as the eigenvalues).

374:    Collective on pep

376:    Input Parameters:
377: +  pep    - the eigensolver context
378: .  etype  - error type
379: -  viewer - optional visualization context

381:    Options Database Keys:
382: +  -pep_error_absolute - print absolute errors of each eigenpair
383: .  -pep_error_relative - print relative errors of each eigenpair
384: -  -pep_error_backward - print backward errors of each eigenpair

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

392:    Level: intermediate

394: .seealso: PEPSolve(), PEPValuesView(), PEPVectorsView()
395: @*/
396: PetscErrorCode PEPErrorView(PEP pep,PEPErrorType etype,PetscViewer viewer)
397: {
398:   PetscBool         isascii;
399:   PetscViewerFormat format;

402:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer);
405:   PEPCheckSolved(pep,1);
406:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
407:   if (!isascii) return 0;

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

427: /*@
428:    PEPErrorViewFromOptions - Processes command line options to determine if/how
429:    the errors of the computed solution are to be viewed.

431:    Collective on pep

433:    Input Parameter:
434: .  pep - the eigensolver context

436:    Level: developer

438: .seealso: PEPErrorView()
439: @*/
440: PetscErrorCode PEPErrorViewFromOptions(PEP pep)
441: {
442:   PetscViewer       viewer;
443:   PetscBool         flg;
444:   static PetscBool  incall = PETSC_FALSE;
445:   PetscViewerFormat format;

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

474: static PetscErrorCode PEPValuesView_DRAW(PEP pep,PetscViewer viewer)
475: {
476:   PetscDraw      draw;
477:   PetscDrawSP    drawsp;
478:   PetscReal      re,im;
479:   PetscInt       i,k;

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

502: static PetscErrorCode PEPValuesView_BINARY(PEP pep,PetscViewer viewer)
503: {
504: #if defined(PETSC_HAVE_COMPLEX)
505:   PetscInt       i,k;
506:   PetscComplex   *ev;
507: #endif

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

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

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

570: static PetscErrorCode PEPValuesView_ASCII(PEP pep,PetscViewer viewer)
571: {
572:   PetscInt       i,k;

574:   PetscViewerASCIIPrintf(viewer,"Eigenvalues = \n");
575:   for (i=0;i<pep->nconv;i++) {
576:     k = pep->perm[i];
577:     PetscViewerASCIIPrintf(viewer,"   ");
578:     SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]);
579:     PetscViewerASCIIPrintf(viewer,"\n");
580:   }
581:   PetscViewerASCIIPrintf(viewer,"\n");
582:   return 0;
583: }

585: static PetscErrorCode PEPValuesView_MATLAB(PEP pep,PetscViewer viewer)
586: {
587:   PetscInt       i,k;
588:   PetscReal      re,im;
589:   const char     *name;

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

609: /*@C
610:    PEPValuesView - Displays the computed eigenvalues in a viewer.

612:    Collective on pep

614:    Input Parameters:
615: +  pep    - the eigensolver context
616: -  viewer - the viewer

618:    Options Database Key:
619: .  -pep_view_values - print computed eigenvalues

621:    Level: intermediate

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

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

667: /*@
668:    PEPValuesViewFromOptions - Processes command line options to determine if/how
669:    the computed eigenvalues are to be viewed.

671:    Collective on pep

673:    Input Parameter:
674: .  pep - the eigensolver context

676:    Level: developer

678: .seealso: PEPValuesView()
679: @*/
680: PetscErrorCode PEPValuesViewFromOptions(PEP pep)
681: {
682:   PetscViewer       viewer;
683:   PetscBool         flg;
684:   static PetscBool  incall = PETSC_FALSE;
685:   PetscViewerFormat format;

687:   if (incall) return 0;
688:   incall = PETSC_TRUE;
689:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_values",&viewer,&format,&flg);
690:   if (flg) {
691:     PetscViewerPushFormat(viewer,format);
692:     PEPValuesView(pep,viewer);
693:     PetscViewerPopFormat(viewer);
694:     PetscViewerDestroy(&viewer);
695:   }
696:   incall = PETSC_FALSE;
697:   return 0;
698: }

700: /*@C
701:    PEPVectorsView - Outputs computed eigenvectors to a viewer.

703:    Collective on pep

705:    Input Parameters:
706: +  pep    - the eigensolver context
707: -  viewer - the viewer

709:    Options Database Key:
710: .  -pep_view_vectors - output eigenvectors.

712:    Note:
713:    If PETSc was configured with real scalars, complex conjugate eigenvectors
714:    will be viewed as two separate real vectors, one containing the real part
715:    and another one containing the imaginary part.

717:    Level: intermediate

719: .seealso: PEPSolve(), PEPValuesView(), PEPErrorView()
720: @*/
721: PetscErrorCode PEPVectorsView(PEP pep,PetscViewer viewer)
722: {
723:   PetscInt       i,k;
724:   Vec            xr,xi=NULL;

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

750: /*@
751:    PEPVectorsViewFromOptions - Processes command line options to determine if/how
752:    the computed eigenvectors are to be viewed.

754:    Collective on pep

756:    Input Parameter:
757: .  pep - the eigensolver context

759:    Level: developer

761: .seealso: PEPVectorsView()
762: @*/
763: PetscErrorCode PEPVectorsViewFromOptions(PEP pep)
764: {
765:   PetscViewer       viewer;
766:   PetscBool         flg = PETSC_FALSE;
767:   static PetscBool  incall = PETSC_FALSE;
768:   PetscViewerFormat format;

770:   if (incall) return 0;
771:   incall = PETSC_TRUE;
772:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_vectors",&viewer,&format,&flg);
773:   if (flg) {
774:     PetscViewerPushFormat(viewer,format);
775:     PEPVectorsView(pep,viewer);
776:     PetscViewerPopFormat(viewer);
777:     PetscViewerDestroy(&viewer);
778:   }
779:   incall = PETSC_FALSE;
780:   return 0;
781: }