Actual source code: epsview.c

slepc-main 2024-11-22
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: /*@
 19:    EPSView - Prints the EPS data structure.

 21:    Collective

 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,isstruct=PETSC_FALSE,flg;
 50:   Mat            A;

 52:   PetscFunctionBegin;
 54:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer));
 56:   PetscCheckSameComm(eps,1,viewer,2);

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

199: /*@
200:    EPSViewFromOptions - View from options

202:    Collective

204:    Input Parameters:
205: +  eps  - the eigensolver context
206: .  obj  - optional object
207: -  name - command line option

209:    Level: intermediate

211: .seealso: EPSView(), EPSCreate()
212: @*/
213: PetscErrorCode EPSViewFromOptions(EPS eps,PetscObject obj,const char name[])
214: {
215:   PetscFunctionBegin;
217:   PetscCall(PetscObjectViewFromOptions((PetscObject)eps,obj,name));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /*@
222:    EPSConvergedReasonView - Displays the reason an EPS solve converged or diverged.

224:    Collective

226:    Input Parameters:
227: +  eps - the eigensolver context
228: -  viewer - the viewer to display the reason

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

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

239:    Level: intermediate

241: .seealso: EPSSetConvergenceTest(), EPSSetTolerances(), EPSGetIterationNumber(), EPSConvergedReasonViewFromOptions()
242: @*/
243: PetscErrorCode EPSConvergedReasonView(EPS eps,PetscViewer viewer)
244: {
245:   PetscBool         isAscii;
246:   PetscViewerFormat format;
247:   PetscInt          nconv;

249:   PetscFunctionBegin;
250:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)eps));
251:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
252:   if (isAscii) {
253:     PetscCall(PetscViewerGetFormat(viewer,&format));
254:     PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)eps)->tablevel));
255:     PetscCall(EPS_GetActualConverged(eps,&nconv));
256:     if (eps->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Linear eigensolve converged (%" PetscInt_FMT " eigenpair%s) due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)eps)->prefix?((PetscObject)eps)->prefix:"",nconv,(nconv>1)?"s":"",EPSConvergedReasons[eps->reason],eps->its));
257:     else if (eps->reason <= 0) PetscCall(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));
258:     PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)eps)->tablevel));
259:   }
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: /*@
264:    EPSConvergedReasonViewFromOptions - Processes command line options to determine if/how
265:    the EPS converged reason is to be viewed.

267:    Collective

269:    Input Parameter:
270: .  eps - the eigensolver context

272:    Level: developer

274: .seealso: EPSConvergedReasonView()
275: @*/
276: PetscErrorCode EPSConvergedReasonViewFromOptions(EPS eps)
277: {
278:   PetscViewer       viewer;
279:   PetscBool         flg;
280:   static PetscBool  incall = PETSC_FALSE;
281:   PetscViewerFormat format;

283:   PetscFunctionBegin;
284:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
285:   incall = PETSC_TRUE;
286:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_converged_reason",&viewer,&format,&flg));
287:   if (flg) {
288:     PetscCall(PetscViewerPushFormat(viewer,format));
289:     PetscCall(EPSConvergedReasonView(eps,viewer));
290:     PetscCall(PetscViewerPopFormat(viewer));
291:     PetscCall(PetscViewerDestroy(&viewer));
292:   }
293:   incall = PETSC_FALSE;
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: static PetscErrorCode EPSErrorView_ASCII(EPS eps,EPSErrorType etype,PetscViewer viewer)
298: {
299:   PetscReal      error;
300:   PetscScalar    kr,ki;
301:   PetscInt       i,j,nvals,nconv;

303:   PetscFunctionBegin;
304:   PetscCall(EPS_GetActualConverged(eps,&nconv));
305:   nvals = (eps->which==EPS_ALL)? eps->nconv: eps->nev;
306:   if (eps->which!=EPS_ALL && nconv<nvals) {
307:     PetscCall(PetscViewerASCIIPrintf(viewer," Problem: less than %" PetscInt_FMT " eigenvalues converged\n\n",eps->nev));
308:     PetscFunctionReturn(PETSC_SUCCESS);
309:   }
310:   if (eps->which==EPS_ALL && !nvals) {
311:     PetscCall(PetscViewerASCIIPrintf(viewer," No eigenvalues have been found\n\n"));
312:     PetscFunctionReturn(PETSC_SUCCESS);
313:   }
314:   for (i=0;i<nvals;i++) {
315:     PetscCall(EPSComputeError(eps,i,etype,&error));
316:     if (error>=5.0*eps->tol) {
317:       PetscCall(PetscViewerASCIIPrintf(viewer," Problem: some of the first %" PetscInt_FMT " relative errors are higher than the tolerance\n\n",nvals));
318:       PetscFunctionReturn(PETSC_SUCCESS);
319:     }
320:   }
321:   if (eps->which==EPS_ALL) PetscCall(PetscViewerASCIIPrintf(viewer," Found %" PetscInt_FMT " eigenvalues, all of them computed up to the required tolerance:",nvals));
322:   else PetscCall(PetscViewerASCIIPrintf(viewer," All requested eigenvalues computed up to the required tolerance:"));
323:   for (i=0;i<=(nvals-1)/8;i++) {
324:     PetscCall(PetscViewerASCIIPrintf(viewer,"\n     "));
325:     for (j=0;j<PetscMin(8,nvals-8*i);j++) {
326:       PetscCall(EPSGetEigenvalue(eps,8*i+j,&kr,&ki));
327:       PetscCall(SlepcPrintEigenvalueASCII(viewer,kr,ki));
328:       if (8*i+j+1<nvals) PetscCall(PetscViewerASCIIPrintf(viewer,", "));
329:     }
330:   }
331:   PetscCall(PetscViewerASCIIPrintf(viewer,"\n\n"));
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: static PetscErrorCode EPSErrorView_DETAIL(EPS eps,EPSErrorType etype,PetscViewer viewer)
336: {
337:   PetscReal      error,re,im;
338:   PetscScalar    kr,ki;
339:   PetscInt       i,nconv;
340:   char           ex[30],sep[]=" ---------------------- --------------------\n";

342:   PetscFunctionBegin;
343:   if (!eps->nconv) PetscFunctionReturn(PETSC_SUCCESS);
344:   switch (etype) {
345:     case EPS_ERROR_ABSOLUTE:
346:       PetscCall(PetscSNPrintf(ex,sizeof(ex),"   ||Ax-k%sx||",eps->isgeneralized?"B":""));
347:       break;
348:     case EPS_ERROR_RELATIVE:
349:       PetscCall(PetscSNPrintf(ex,sizeof(ex),"||Ax-k%sx||/||kx||",eps->isgeneralized?"B":""));
350:       break;
351:     case EPS_ERROR_BACKWARD:
352:       PetscCall(PetscSNPrintf(ex,sizeof(ex),"    eta(x,k)"));
353:       break;
354:   }
355:   PetscCall(PetscViewerASCIIPrintf(viewer,"%s            k             %s\n%s",sep,ex,sep));
356:   PetscCall(EPS_GetActualConverged(eps,&nconv));
357:   for (i=0;i<nconv;i++) {
358:     PetscCall(EPSGetEigenvalue(eps,i,&kr,&ki));
359:     PetscCall(EPSComputeError(eps,i,etype,&error));
360: #if defined(PETSC_USE_COMPLEX)
361:     re = PetscRealPart(kr);
362:     im = PetscImaginaryPart(kr);
363: #else
364:     re = kr;
365:     im = ki;
366: #endif
367:     if (im!=0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"  % 9f%+9fi      %12g\n",(double)re,(double)im,(double)error));
368:     else PetscCall(PetscViewerASCIIPrintf(viewer,"    % 12f           %12g\n",(double)re,(double)error));
369:   }
370:   PetscCall(PetscViewerASCIIPrintf(viewer,"%s",sep));
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: static PetscErrorCode EPSErrorView_MATLAB(EPS eps,EPSErrorType etype,PetscViewer viewer)
375: {
376:   PetscReal      error;
377:   PetscInt       i,nconv;
378:   const char     *name;

380:   PetscFunctionBegin;
381:   PetscCall(PetscObjectGetName((PetscObject)eps,&name));
382:   PetscCall(PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name));
383:   PetscCall(EPS_GetActualConverged(eps,&nconv));
384:   for (i=0;i<nconv;i++) {
385:     PetscCall(EPSComputeError(eps,i,etype,&error));
386:     PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error));
387:   }
388:   PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
389:   PetscFunctionReturn(PETSC_SUCCESS);
390: }

392: /*@
393:    EPSErrorView - Displays the errors associated with the computed solution
394:    (as well as the eigenvalues).

396:    Collective

398:    Input Parameters:
399: +  eps    - the eigensolver context
400: .  etype  - error type
401: -  viewer - optional visualization context

403:    Options Database Keys:
404: +  -eps_error_absolute - print absolute errors of each eigenpair
405: .  -eps_error_relative - print relative errors of each eigenpair
406: -  -eps_error_backward - print backward errors of each eigenpair

408:    Notes:
409:    By default, this function checks the error of all eigenpairs and prints
410:    the eigenvalues if all of them are below the requested tolerance.
411:    If the viewer has format=PETSC_VIEWER_ASCII_INFO_DETAIL then a table with
412:    eigenvalues and corresponding errors is printed.

414:    Level: intermediate

416: .seealso: EPSSolve(), EPSValuesView(), EPSVectorsView()
417: @*/
418: PetscErrorCode EPSErrorView(EPS eps,EPSErrorType etype,PetscViewer viewer)
419: {
420:   PetscBool         isascii;
421:   PetscViewerFormat format;

423:   PetscFunctionBegin;
425:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer));
427:   PetscCheckSameComm(eps,1,viewer,3);
428:   EPSCheckSolved(eps,1);
429:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
430:   if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);

432:   PetscCall(PetscViewerGetFormat(viewer,&format));
433:   switch (format) {
434:     case PETSC_VIEWER_DEFAULT:
435:     case PETSC_VIEWER_ASCII_INFO:
436:       PetscCall(EPSErrorView_ASCII(eps,etype,viewer));
437:       break;
438:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
439:       PetscCall(EPSErrorView_DETAIL(eps,etype,viewer));
440:       break;
441:     case PETSC_VIEWER_ASCII_MATLAB:
442:       PetscCall(EPSErrorView_MATLAB(eps,etype,viewer));
443:       break;
444:     default:
445:       PetscCall(PetscInfo(eps,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
446:   }
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: /*@
451:    EPSErrorViewFromOptions - Processes command line options to determine if/how
452:    the errors of the computed solution are to be viewed.

454:    Collective

456:    Input Parameter:
457: .  eps - the eigensolver context

459:    Level: developer

461: .seealso: EPSErrorView()
462: @*/
463: PetscErrorCode EPSErrorViewFromOptions(EPS eps)
464: {
465:   PetscViewer       viewer;
466:   PetscBool         flg;
467:   static PetscBool  incall = PETSC_FALSE;
468:   PetscViewerFormat format;

470:   PetscFunctionBegin;
471:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
472:   incall = PETSC_TRUE;
473:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_absolute",&viewer,&format,&flg));
474:   if (flg) {
475:     PetscCall(PetscViewerPushFormat(viewer,format));
476:     PetscCall(EPSErrorView(eps,EPS_ERROR_ABSOLUTE,viewer));
477:     PetscCall(PetscViewerPopFormat(viewer));
478:     PetscCall(PetscViewerDestroy(&viewer));
479:   }
480:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_relative",&viewer,&format,&flg));
481:   if (flg) {
482:     PetscCall(PetscViewerPushFormat(viewer,format));
483:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
484:     PetscCall(PetscViewerPopFormat(viewer));
485:     PetscCall(PetscViewerDestroy(&viewer));
486:   }
487:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_error_backward",&viewer,&format,&flg));
488:   if (flg) {
489:     PetscCall(PetscViewerPushFormat(viewer,format));
490:     PetscCall(EPSErrorView(eps,EPS_ERROR_BACKWARD,viewer));
491:     PetscCall(PetscViewerPopFormat(viewer));
492:     PetscCall(PetscViewerDestroy(&viewer));
493:   }
494:   incall = PETSC_FALSE;
495:   PetscFunctionReturn(PETSC_SUCCESS);
496: }

498: static PetscErrorCode EPSValuesView_DRAW(EPS eps,PetscViewer viewer)
499: {
500:   PetscDraw      draw;
501:   PetscDrawSP    drawsp;
502:   PetscReal      re,im;
503:   PetscScalar    kr,ki;
504:   PetscInt       i,nconv;

506:   PetscFunctionBegin;
507:   if (!eps->nconv) PetscFunctionReturn(PETSC_SUCCESS);
508:   PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
509:   PetscCall(PetscDrawSetTitle(draw,"Computed Eigenvalues"));
510:   PetscCall(PetscDrawSPCreate(draw,1,&drawsp));
511:   PetscCall(EPS_GetActualConverged(eps,&nconv));
512:   for (i=0;i<nconv;i++) {
513:     PetscCall(EPSGetEigenvalue(eps,i,&kr,&ki));
514: #if defined(PETSC_USE_COMPLEX)
515:     re = PetscRealPart(kr);
516:     im = PetscImaginaryPart(kr);
517: #else
518:     re = kr;
519:     im = ki;
520: #endif
521:     PetscCall(PetscDrawSPAddPoint(drawsp,&re,&im));
522:   }
523:   PetscCall(PetscDrawSPDraw(drawsp,PETSC_TRUE));
524:   PetscCall(PetscDrawSPSave(drawsp));
525:   PetscCall(PetscDrawSPDestroy(&drawsp));
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }

529: static PetscErrorCode EPSValuesView_BINARY(EPS eps,PetscViewer viewer)
530: {
531:   PetscScalar    kr,ki;
532: #if defined(PETSC_HAVE_COMPLEX)
533:   PetscInt       i,nconv;
534:   PetscComplex   *ev;
535: #endif

537:   PetscFunctionBegin;
538: #if defined(PETSC_HAVE_COMPLEX)
539:   PetscCall(EPS_GetActualConverged(eps,&nconv));
540:   PetscCall(PetscMalloc1(nconv,&ev));
541:   for (i=0;i<nconv;i++) {
542:     PetscCall(EPSGetEigenvalue(eps,i,&kr,&ki));
543: #if defined(PETSC_USE_COMPLEX)
544:     ev[i] = kr;
545: #else
546:     ev[i] = PetscCMPLX(kr,ki);
547: #endif
548:   }
549:   PetscCall(PetscViewerBinaryWrite(viewer,ev,nconv,PETSC_COMPLEX));
550:   PetscCall(PetscFree(ev));
551: #endif
552:   PetscFunctionReturn(PETSC_SUCCESS);
553: }

555: #if defined(PETSC_HAVE_HDF5)
556: static PetscErrorCode EPSValuesView_HDF5(EPS eps,PetscViewer viewer)
557: {
558:   PetscInt       i,n,N,nconv;
559:   PetscScalar    eig;
560:   PetscMPIInt    rank;
561:   Vec            v;
562:   char           vname[30];
563:   const char     *ename;

565:   PetscFunctionBegin;
566:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
567:   PetscCall(EPS_GetActualConverged(eps,&nconv));
568:   N = nconv;
569:   n = rank? 0: N;
570:   /* create a vector containing the eigenvalues */
571:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)eps),n,N,&v));
572:   PetscCall(PetscObjectGetName((PetscObject)eps,&ename));
573:   PetscCall(PetscSNPrintf(vname,sizeof(vname),"eigr_%s",ename));
574:   PetscCall(PetscObjectSetName((PetscObject)v,vname));
575:   if (!rank) {
576:     for (i=0;i<nconv;i++) {
577:       PetscCall(EPSGetEigenvalue(eps,i,&eig,NULL));
578:       PetscCall(VecSetValue(v,i,eig,INSERT_VALUES));
579:     }
580:   }
581:   PetscCall(VecAssemblyBegin(v));
582:   PetscCall(VecAssemblyEnd(v));
583:   PetscCall(VecView(v,viewer));
584: #if !defined(PETSC_USE_COMPLEX)
585:   /* in real scalars write the imaginary part as a separate vector */
586:   PetscCall(PetscSNPrintf(vname,sizeof(vname),"eigi_%s",ename));
587:   PetscCall(PetscObjectSetName((PetscObject)v,vname));
588:   if (!rank) {
589:     for (i=0;i<nconv;i++) {
590:       PetscCall(EPSGetEigenvalue(eps,i,NULL,&eig));
591:       PetscCall(VecSetValue(v,i,eig,INSERT_VALUES));
592:     }
593:   }
594:   PetscCall(VecAssemblyBegin(v));
595:   PetscCall(VecAssemblyEnd(v));
596:   PetscCall(VecView(v,viewer));
597: #endif
598:   PetscCall(VecDestroy(&v));
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }
601: #endif

603: static PetscErrorCode EPSValuesView_ASCII(EPS eps,PetscViewer viewer)
604: {
605:   PetscInt       i,nconv;
606:   PetscScalar    kr,ki;

608:   PetscFunctionBegin;
609:   PetscCall(PetscViewerASCIIPrintf(viewer,"Eigenvalues = \n"));
610:   PetscCall(EPS_GetActualConverged(eps,&nconv));
611:   for (i=0;i<nconv;i++) {
612:     PetscCall(EPSGetEigenvalue(eps,i,&kr,&ki));
613:     PetscCall(PetscViewerASCIIPrintf(viewer,"   "));
614:     PetscCall(SlepcPrintEigenvalueASCII(viewer,kr,ki));
615:     PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
616:   }
617:   PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
618:   PetscFunctionReturn(PETSC_SUCCESS);
619: }

621: static PetscErrorCode EPSValuesView_MATLAB(EPS eps,PetscViewer viewer)
622: {
623:   PetscInt       i,nconv;
624:   PetscReal      re,im;
625:   PetscScalar    kr,ki;
626:   const char     *name;

628:   PetscFunctionBegin;
629:   PetscCall(PetscObjectGetName((PetscObject)eps,&name));
630:   PetscCall(PetscViewerASCIIPrintf(viewer,"Lambda_%s = [\n",name));
631:   PetscCall(EPS_GetActualConverged(eps,&nconv));
632:   for (i=0;i<nconv;i++) {
633:     PetscCall(EPSGetEigenvalue(eps,i,&kr,&ki));
634: #if defined(PETSC_USE_COMPLEX)
635:     re = PetscRealPart(kr);
636:     im = PetscImaginaryPart(kr);
637: #else
638:     re = kr;
639:     im = ki;
640: #endif
641:     if (im!=0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei\n",(double)re,(double)im));
642:     else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)re));
643:   }
644:   PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
645:   PetscFunctionReturn(PETSC_SUCCESS);
646: }

648: /*@
649:    EPSValuesView - Displays the computed eigenvalues in a viewer.

651:    Collective

653:    Input Parameters:
654: +  eps    - the eigensolver context
655: -  viewer - the viewer

657:    Options Database Key:
658: .  -eps_view_values - print computed eigenvalues

660:    Level: intermediate

662: .seealso: EPSSolve(), EPSVectorsView(), EPSErrorView()
663: @*/
664: PetscErrorCode EPSValuesView(EPS eps,PetscViewer viewer)
665: {
666:   PetscBool         isascii,isdraw,isbinary;
667:   PetscViewerFormat format;
668: #if defined(PETSC_HAVE_HDF5)
669:   PetscBool         ishdf5;
670: #endif

672:   PetscFunctionBegin;
674:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer));
676:   PetscCheckSameComm(eps,1,viewer,2);
677:   EPSCheckSolved(eps,1);
678:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
679:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary));
680: #if defined(PETSC_HAVE_HDF5)
681:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5));
682: #endif
683:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
684:   if (isdraw) PetscCall(EPSValuesView_DRAW(eps,viewer));
685:   else if (isbinary) PetscCall(EPSValuesView_BINARY(eps,viewer));
686: #if defined(PETSC_HAVE_HDF5)
687:   else if (ishdf5) PetscCall(EPSValuesView_HDF5(eps,viewer));
688: #endif
689:   else if (isascii) {
690:     PetscCall(PetscViewerGetFormat(viewer,&format));
691:     switch (format) {
692:       case PETSC_VIEWER_DEFAULT:
693:       case PETSC_VIEWER_ASCII_INFO:
694:       case PETSC_VIEWER_ASCII_INFO_DETAIL:
695:         PetscCall(EPSValuesView_ASCII(eps,viewer));
696:         break;
697:       case PETSC_VIEWER_ASCII_MATLAB:
698:         PetscCall(EPSValuesView_MATLAB(eps,viewer));
699:         break;
700:       default:
701:         PetscCall(PetscInfo(eps,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
702:     }
703:   }
704:   PetscFunctionReturn(PETSC_SUCCESS);
705: }

707: /*@
708:    EPSValuesViewFromOptions - Processes command line options to determine if/how
709:    the computed eigenvalues are to be viewed.

711:    Collective

713:    Input Parameters:
714: .  eps - the eigensolver context

716:    Level: developer

718: .seealso: EPSValuesView()
719: @*/
720: PetscErrorCode EPSValuesViewFromOptions(EPS eps)
721: {
722:   PetscViewer       viewer;
723:   PetscBool         flg;
724:   static PetscBool  incall = PETSC_FALSE;
725:   PetscViewerFormat format;

727:   PetscFunctionBegin;
728:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
729:   incall = PETSC_TRUE;
730:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_values",&viewer,&format,&flg));
731:   if (flg) {
732:     PetscCall(PetscViewerPushFormat(viewer,format));
733:     PetscCall(EPSValuesView(eps,viewer));
734:     PetscCall(PetscViewerPopFormat(viewer));
735:     PetscCall(PetscViewerDestroy(&viewer));
736:   }
737:   incall = PETSC_FALSE;
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: /*@
742:    EPSVectorsView - Outputs computed eigenvectors to a viewer.

744:    Collective

746:    Input Parameters:
747: +  eps    - the eigensolver context
748: -  viewer - the viewer

750:    Options Database Key:
751: .  -eps_view_vectors - output eigenvectors.

753:    Notes:
754:    If PETSc was configured with real scalars, complex conjugate eigenvectors
755:    will be viewed as two separate real vectors, one containing the real part
756:    and another one containing the imaginary part.

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

762:    Level: intermediate

764: .seealso: EPSSolve(), EPSValuesView(), EPSErrorView()
765: @*/
766: PetscErrorCode EPSVectorsView(EPS eps,PetscViewer viewer)
767: {
768:   PetscInt       i,nconv;
769:   Vec            xr,xi=NULL;

771:   PetscFunctionBegin;
773:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer));
775:   PetscCheckSameComm(eps,1,viewer,2);
776:   EPSCheckSolved(eps,1);
777:   PetscCall(EPS_GetActualConverged(eps,&nconv));
778:   if (nconv) {
779:     PetscCall(BVCreateVec(eps->V,&xr));
780: #if !defined(PETSC_USE_COMPLEX)
781:     PetscCall(BVCreateVec(eps->V,&xi));
782: #endif
783:     for (i=0;i<nconv;i++) {
784:       PetscCall(EPSGetEigenvector(eps,i,xr,xi));
785:       PetscCall(SlepcViewEigenvector(viewer,xr,xi,"X",i,(PetscObject)eps));
786:       if (eps->twosided || eps->problem_type==EPS_BSE) {
787:         PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
788:         PetscCall(SlepcViewEigenvector(viewer,xr,xi,"Y",i,(PetscObject)eps));
789:       }
790:     }
791:     PetscCall(VecDestroy(&xr));
792: #if !defined(PETSC_USE_COMPLEX)
793:     PetscCall(VecDestroy(&xi));
794: #endif
795:   }
796:   PetscFunctionReturn(PETSC_SUCCESS);
797: }

799: /*@
800:    EPSVectorsViewFromOptions - Processes command line options to determine if/how
801:    the computed eigenvectors are to be viewed.

803:    Collective

805:    Input Parameter:
806: .  eps - the eigensolver context

808:    Level: developer

810: .seealso: EPSVectorsView()
811: @*/
812: PetscErrorCode EPSVectorsViewFromOptions(EPS eps)
813: {
814:   PetscViewer       viewer;
815:   PetscBool         flg = PETSC_FALSE;
816:   static PetscBool  incall = PETSC_FALSE;
817:   PetscViewerFormat format;

819:   PetscFunctionBegin;
820:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
821:   incall = PETSC_TRUE;
822:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)eps),((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_vectors",&viewer,&format,&flg));
823:   if (flg) {
824:     PetscCall(PetscViewerPushFormat(viewer,format));
825:     PetscCall(EPSVectorsView(eps,viewer));
826:     PetscCall(PetscViewerPopFormat(viewer));
827:     PetscCall(PetscViewerDestroy(&viewer));
828:   }
829:   incall = PETSC_FALSE;
830:   PetscFunctionReturn(PETSC_SUCCESS);
831: }