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

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

 18: /*@C
 19:    EPSView - Prints the EPS data structure.

 21:    Collective on eps

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

 27:    Options Database Key:
 28: .  -eps_view -  Calls EPSView() at end of EPSSolve()

 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 EPSView(EPS eps,PetscViewer viewer)
 46: {
 47:   const char     *type=NULL,*extr=NULL,*bal=NULL;
 48:   char           str[50];
 49:   PetscBool      isascii,isexternal,istrivial;

 52:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer);

 56:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 57:   if (isascii) {
 58:     PetscObjectPrintClassNamePrefixType((PetscObject)eps,viewer);
 59:     PetscViewerASCIIPushTab(viewer);
 60:     PetscTryTypeMethod(eps,view,viewer);
 61:     PetscViewerASCIIPopTab(viewer);
 62:     if (eps->problem_type) {
 63:       switch (eps->problem_type) {
 64:         case EPS_HEP:    type = SLEPC_STRING_HERMITIAN " eigenvalue problem"; break;
 65:         case EPS_GHEP:   type = "generalized " SLEPC_STRING_HERMITIAN " eigenvalue problem"; break;
 66:         case EPS_NHEP:   type = "non-" SLEPC_STRING_HERMITIAN " eigenvalue problem"; break;
 67:         case EPS_GNHEP:  type = "generalized non-" SLEPC_STRING_HERMITIAN " eigenvalue problem"; break;
 68:         case EPS_PGNHEP: type = "generalized non-" SLEPC_STRING_HERMITIAN " eigenvalue problem with " SLEPC_STRING_HERMITIAN " positive definite B"; break;
 69:         case EPS_GHIEP:  type = "generalized " SLEPC_STRING_HERMITIAN "-indefinite eigenvalue problem"; break;
 70:       }
 71:     } else type = "not yet set";
 72:     PetscViewerASCIIPrintf(viewer,"  problem type: %s\n",type);
 73:     if (eps->extraction) {
 74:       switch (eps->extraction) {
 75:         case EPS_RITZ:              extr = "Rayleigh-Ritz"; break;
 76:         case EPS_HARMONIC:          extr = "harmonic Ritz"; break;
 77:         case EPS_HARMONIC_RELATIVE: extr = "relative harmonic Ritz"; break;
 78:         case EPS_HARMONIC_RIGHT:    extr = "right harmonic Ritz"; break;
 79:         case EPS_HARMONIC_LARGEST:  extr = "largest harmonic Ritz"; break;
 80:         case EPS_REFINED:           extr = "refined Ritz"; break;
 81:         case EPS_REFINED_HARMONIC:  extr = "refined harmonic Ritz"; break;
 82:       }
 83:       PetscViewerASCIIPrintf(viewer,"  extraction type: %s\n",extr);
 84:     }
 85:     if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE) {
 86:       switch (eps->balance) {
 87:         case EPS_BALANCE_NONE:    break;
 88:         case EPS_BALANCE_ONESIDE: bal = "one-sided Krylov"; break;
 89:         case EPS_BALANCE_TWOSIDE: bal = "two-sided Krylov"; break;
 90:         case EPS_BALANCE_USER:    bal = "user-defined matrix"; break;
 91:       }
 92:       PetscViewerASCIIPrintf(viewer,"  balancing enabled: %s",bal);
 93:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 94:       if (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE) PetscViewerASCIIPrintf(viewer,", with its=%" PetscInt_FMT,eps->balance_its);
 95:       if (eps->balance==EPS_BALANCE_TWOSIDE && eps->balance_cutoff!=0.0) PetscViewerASCIIPrintf(viewer," and cutoff=%g",(double)eps->balance_cutoff);
 96:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
 97:       PetscViewerASCIIPrintf(viewer,"\n");
 98:     }
 99:     PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: ");
100:     SlepcSNPrintfScalar(str,sizeof(str),eps->target,PETSC_FALSE);
101:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
102:     if (!eps->which) PetscViewerASCIIPrintf(viewer,"not yet set\n");
103:     else switch (eps->which) {
104:       case EPS_WHICH_USER:
105:         PetscViewerASCIIPrintf(viewer,"user defined\n");
106:         break;
107:       case EPS_TARGET_MAGNITUDE:
108:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (in magnitude)\n",str);
109:         break;
110:       case EPS_TARGET_REAL:
111:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the real axis)\n",str);
112:         break;
113:       case EPS_TARGET_IMAGINARY:
114:         PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the imaginary axis)\n",str);
115:         break;
116:       case EPS_LARGEST_MAGNITUDE:
117:         PetscViewerASCIIPrintf(viewer,"largest eigenvalues in magnitude\n");
118:         break;
119:       case EPS_SMALLEST_MAGNITUDE:
120:         PetscViewerASCIIPrintf(viewer,"smallest eigenvalues in magnitude\n");
121:         break;
122:       case EPS_LARGEST_REAL:
123:         PetscViewerASCIIPrintf(viewer,"largest real parts\n");
124:         break;
125:       case EPS_SMALLEST_REAL:
126:         PetscViewerASCIIPrintf(viewer,"smallest real parts\n");
127:         break;
128:       case EPS_LARGEST_IMAGINARY:
129:         PetscViewerASCIIPrintf(viewer,"largest imaginary parts\n");
130:         break;
131:       case EPS_SMALLEST_IMAGINARY:
132:         PetscViewerASCIIPrintf(viewer,"smallest imaginary parts\n");
133:         break;
134:       case EPS_ALL:
135:         if (eps->inta || eps->intb) PetscViewerASCIIPrintf(viewer,"all eigenvalues in interval [%g,%g]\n",(double)eps->inta,(double)eps->intb);
136:         else PetscViewerASCIIPrintf(viewer,"all eigenvalues in the region\n");
137:         break;
138:     }
139:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
140:     if (eps->twosided && eps->problem_type!=EPS_HEP && eps->problem_type!=EPS_GHEP) PetscViewerASCIIPrintf(viewer,"  using two-sided variant (for left eigenvectors)\n");
141:     if (eps->purify) PetscViewerASCIIPrintf(viewer,"  postprocessing eigenvectors with purification\n");
142:     if (eps->trueres) PetscViewerASCIIPrintf(viewer,"  computing true residuals explicitly\n");
143:     if (eps->trackall) PetscViewerASCIIPrintf(viewer,"  computing all residuals (for tracking convergence)\n");
144:     PetscViewerASCIIPrintf(viewer,"  number of eigenvalues (nev): %" PetscInt_FMT "\n",eps->nev);
145:     PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",eps->ncv);
146:     PetscViewerASCIIPrintf(viewer,"  maximum dimension of projected problem (mpd): %" PetscInt_FMT "\n",eps->mpd);
147:     PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",eps->max_it);
148:     PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)eps->tol);
149:     PetscViewerASCIIPrintf(viewer,"  convergence test: ");
150:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
151:     switch (eps->conv) {
152:     case EPS_CONV_ABS:
153:       PetscViewerASCIIPrintf(viewer,"absolute\n");break;
154:     case EPS_CONV_REL:
155:       PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue\n");break;
156:     case EPS_CONV_NORM:
157:       PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue and matrix norms\n");
158:       PetscViewerASCIIPrintf(viewer,"  computed matrix norms: norm(A)=%g",(double)eps->nrma);
159:       if (eps->isgeneralized) PetscViewerASCIIPrintf(viewer,", norm(B)=%g",(double)eps->nrmb);
160:       PetscViewerASCIIPrintf(viewer,"\n");
161:       break;
162:     case EPS_CONV_USER:
163:       PetscViewerASCIIPrintf(viewer,"user-defined\n");break;
164:     }
165:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
166:     if (eps->nini) PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial space: %" PetscInt_FMT "\n",PetscAbs(eps->nini));
167:     if (eps->ninil) PetscViewerASCIIPrintf(viewer,"  dimension of user-provided left initial space: %" PetscInt_FMT "\n",PetscAbs(eps->ninil));
168:     if (eps->nds) PetscViewerASCIIPrintf(viewer,"  dimension of user-provided deflation space: %" PetscInt_FMT "\n",PetscAbs(eps->nds));
169:   } else PetscTryTypeMethod(eps,view,viewer);
170:   PetscObjectTypeCompareAny((PetscObject)eps,&isexternal,EPSARPACK,EPSBLOPEX,EPSELEMENTAL,EPSFEAST,EPSPRIMME,EPSSCALAPACK,EPSELPA,EPSEVSL,EPSTRLAN,"");
171:   if (!isexternal) {
172:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
173:     if (!eps->V) EPSGetBV(eps,&eps->V);
174:     BVView(eps->V,viewer);
175:     if (eps->rg) {
176:       RGIsTrivial(eps->rg,&istrivial);
177:       if (!istrivial) RGView(eps->rg,viewer);
178:     }
179:     if (eps->useds) {
180:       if (!eps->ds) EPSGetDS(eps,&eps->ds);
181:       DSView(eps->ds,viewer);
182:     }
183:     PetscViewerPopFormat(viewer);
184:   }
185:   if (!eps->st) EPSGetST(eps,&eps->st);
186:   STView(eps->st,viewer);
187:   return 0;
188: }

190: /*@C
191:    EPSViewFromOptions - View from options

193:    Collective on EPS

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

200:    Level: intermediate

202: .seealso: EPSView(), EPSCreate()
203: @*/
204: PetscErrorCode EPSViewFromOptions(EPS eps,PetscObject obj,const char name[])
205: {
207:   PetscObjectViewFromOptions((PetscObject)eps,obj,name);
208:   return 0;
209: }

211: /*@C
212:    EPSConvergedReasonView - Displays the reason an EPS solve converged or diverged.

214:    Collective on eps

216:    Input Parameters:
217: +  eps - the eigensolver context
218: -  viewer - the viewer to display the reason

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

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

229:    Level: intermediate

231: .seealso: EPSSetConvergenceTest(), EPSSetTolerances(), EPSGetIterationNumber(), EPSConvergedReasonViewFromOptions()
232: @*/
233: PetscErrorCode EPSConvergedReasonView(EPS eps,PetscViewer viewer)
234: {
235:   PetscBool         isAscii;
236:   PetscViewerFormat format;

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

250: /*@
251:    EPSConvergedReasonViewFromOptions - Processes command line options to determine if/how
252:    the EPS converged reason is to be viewed.

254:    Collective on eps

256:    Input Parameter:
257: .  eps - the eigensolver context

259:    Level: developer

261: .seealso: EPSConvergedReasonView()
262: @*/
263: PetscErrorCode EPSConvergedReasonViewFromOptions(EPS eps)
264: {
265:   PetscViewer       viewer;
266:   PetscBool         flg;
267:   static PetscBool  incall = PETSC_FALSE;
268:   PetscViewerFormat format;

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

283: static PetscErrorCode EPSErrorView_ASCII(EPS eps,EPSErrorType etype,PetscViewer viewer)
284: {
285:   PetscReal      error;
286:   PetscInt       i,j,k,nvals;

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

318: static PetscErrorCode EPSErrorView_DETAIL(EPS eps,EPSErrorType etype,PetscViewer viewer)
319: {
320:   PetscReal      error,re,im;
321:   PetscScalar    kr,ki;
322:   PetscInt       i;
323:   char           ex[30],sep[]=" ---------------------- --------------------\n";

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

355: static PetscErrorCode EPSErrorView_MATLAB(EPS eps,EPSErrorType etype,PetscViewer viewer)
356: {
357:   PetscReal      error;
358:   PetscInt       i;
359:   const char     *name;

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

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

375:    Collective on eps

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

382:    Options Database Keys:
383: +  -eps_error_absolute - print absolute errors of each eigenpair
384: .  -eps_error_relative - print relative errors of each eigenpair
385: -  -eps_error_backward - print backward errors of each eigenpair

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

393:    Level: intermediate

395: .seealso: EPSSolve(), EPSValuesView(), EPSVectorsView()
396: @*/
397: PetscErrorCode EPSErrorView(EPS eps,EPSErrorType etype,PetscViewer viewer)
398: {
399:   PetscBool         isascii;
400:   PetscViewerFormat format;

403:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer);
406:   EPSCheckSolved(eps,1);
407:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
408:   if (!isascii) return 0;

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

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

432:    Collective on eps

434:    Input Parameter:
435: .  eps - the eigensolver context

437:    Level: developer

439: .seealso: EPSErrorView()
440: @*/
441: PetscErrorCode EPSErrorViewFromOptions(EPS eps)
442: {
443:   PetscViewer       viewer;
444:   PetscBool         flg;
445:   static PetscBool  incall = PETSC_FALSE;
446:   PetscViewerFormat format;

448:   if (incall) return 0;
449:   incall = PETSC_TRUE;
450:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_absolute",&viewer,&format,&flg);
451:   if (flg) {
452:     PetscViewerPushFormat(viewer,format);
453:     EPSErrorView(eps,EPS_ERROR_ABSOLUTE,viewer);
454:     PetscViewerPopFormat(viewer);
455:     PetscViewerDestroy(&viewer);
456:   }
457:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_relative",&viewer,&format,&flg);
458:   if (flg) {
459:     PetscViewerPushFormat(viewer,format);
460:     EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer);
461:     PetscViewerPopFormat(viewer);
462:     PetscViewerDestroy(&viewer);
463:   }
464:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_backward",&viewer,&format,&flg);
465:   if (flg) {
466:     PetscViewerPushFormat(viewer,format);
467:     EPSErrorView(eps,EPS_ERROR_BACKWARD,viewer);
468:     PetscViewerPopFormat(viewer);
469:     PetscViewerDestroy(&viewer);
470:   }
471:   incall = PETSC_FALSE;
472:   return 0;
473: }

475: static PetscErrorCode EPSValuesView_DRAW(EPS eps,PetscViewer viewer)
476: {
477:   PetscDraw      draw;
478:   PetscDrawSP    drawsp;
479:   PetscReal      re,im;
480:   PetscInt       i,k;

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

503: static PetscErrorCode EPSValuesView_BINARY(EPS eps,PetscViewer viewer)
504: {
505: #if defined(PETSC_HAVE_COMPLEX)
506:   PetscInt       i,k;
507:   PetscComplex   *ev;
508: #endif

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

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

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

571: static PetscErrorCode EPSValuesView_ASCII(EPS eps,PetscViewer viewer)
572: {
573:   PetscInt       i,k;

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

586: static PetscErrorCode EPSValuesView_MATLAB(EPS eps,PetscViewer viewer)
587: {
588:   PetscInt       i,k;
589:   PetscReal      re,im;
590:   const char     *name;

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

610: /*@C
611:    EPSValuesView - Displays the computed eigenvalues in a viewer.

613:    Collective on eps

615:    Input Parameters:
616: +  eps    - the eigensolver context
617: -  viewer - the viewer

619:    Options Database Key:
620: .  -eps_view_values - print computed eigenvalues

622:    Level: intermediate

624: .seealso: EPSSolve(), EPSVectorsView(), EPSErrorView()
625: @*/
626: PetscErrorCode EPSValuesView(EPS eps,PetscViewer viewer)
627: {
628:   PetscBool         isascii,isdraw,isbinary;
629:   PetscViewerFormat format;
630: #if defined(PETSC_HAVE_HDF5)
631:   PetscBool         ishdf5;
632: #endif

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

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

672:    Collective on eps

674:    Input Parameters:
675: .  eps - the eigensolver context

677:    Level: developer

679: .seealso: EPSValuesView()
680: @*/
681: PetscErrorCode EPSValuesViewFromOptions(EPS eps)
682: {
683:   PetscViewer       viewer;
684:   PetscBool         flg;
685:   static PetscBool  incall = PETSC_FALSE;
686:   PetscViewerFormat format;

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

701: /*@C
702:    EPSVectorsView - Outputs computed eigenvectors to a viewer.

704:    Collective on eps

706:    Input Parameters:
707: +  eps    - the eigensolver context
708: -  viewer - the viewer

710:    Options Database Key:
711: .  -eps_view_vectors - output eigenvectors.

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

718:    If left eigenvectors were computed with a two-sided eigensolver, the right
719:    and left eigenvectors are interleaved, that is, the vectors are output in
720:    the following order X0, Y0, X1, Y1, X2, Y2, ...

722:    Level: intermediate

724: .seealso: EPSSolve(), EPSValuesView(), EPSErrorView()
725: @*/
726: PetscErrorCode EPSVectorsView(EPS eps,PetscViewer viewer)
727: {
728:   PetscInt       i,k;
729:   Vec            xr,xi=NULL;

732:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer);
735:   EPSCheckSolved(eps,1);
736:   if (eps->nconv) {
737:     EPSComputeVectors(eps);
738:     BVCreateVec(eps->V,&xr);
739: #if !defined(PETSC_USE_COMPLEX)
740:     BVCreateVec(eps->V,&xi);
741: #endif
742:     for (i=0;i<eps->nconv;i++) {
743:       k = eps->perm[i];
744:       BV_GetEigenvector(eps->V,k,eps->eigi[k],xr,xi);
745:       SlepcViewEigenvector(viewer,xr,xi,"X",i,(PetscObject)eps);
746:       if (eps->twosided) {
747:         BV_GetEigenvector(eps->W,k,eps->eigi[k],xr,xi);
748:         SlepcViewEigenvector(viewer,xr,xi,"Y",i,(PetscObject)eps);
749:       }
750:     }
751:     VecDestroy(&xr);
752: #if !defined(PETSC_USE_COMPLEX)
753:     VecDestroy(&xi);
754: #endif
755:   }
756:   return 0;
757: }

759: /*@
760:    EPSVectorsViewFromOptions - Processes command line options to determine if/how
761:    the computed eigenvectors are to be viewed.

763:    Collective on eps

765:    Input Parameter:
766: .  eps - the eigensolver context

768:    Level: developer

770: .seealso: EPSVectorsView()
771: @*/
772: PetscErrorCode EPSVectorsViewFromOptions(EPS eps)
773: {
774:   PetscViewer       viewer;
775:   PetscBool         flg = PETSC_FALSE;
776:   static PetscBool  incall = PETSC_FALSE;
777:   PetscViewerFormat format;

779:   if (incall) return 0;
780:   incall = PETSC_TRUE;
781:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_vectors",&viewer,&format,&flg);
782:   if (flg) {
783:     PetscViewerPushFormat(viewer,format);
784:     EPSVectorsView(eps,viewer);
785:     PetscViewerPopFormat(viewer);
786:     PetscViewerDestroy(&viewer);
787:   }
788:   incall = PETSC_FALSE;
789:   return 0;
790: }