Actual source code: epsview.c

slepc-main 2025-01-19
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:     if (eps->stop==EPS_STOP_THRESHOLD) PetscCall(PetscViewerASCIIPrintf(viewer,"  computing eigenvalues %s the threshold: %g%s\n",(eps->which==EPS_SMALLEST_MAGNITUDE||eps->which==EPS_SMALLEST_REAL)?"below":"above",(double)eps->thres,eps->threlative?" (relative)":""));
154:     if (eps->nev) PetscCall(PetscViewerASCIIPrintf(viewer,"  number of eigenvalues (nev): %" PetscInt_FMT "\n",eps->nev));
155:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",eps->ncv));
156:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum dimension of projected problem (mpd): %" PetscInt_FMT "\n",eps->mpd));
157:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",eps->max_it));
158:     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)eps->tol));
159:     PetscCall(PetscViewerASCIIPrintf(viewer,"  convergence test: "));
160:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
161:     switch (eps->conv) {
162:     case EPS_CONV_ABS:
163:       PetscCall(PetscViewerASCIIPrintf(viewer,"absolute\n"));break;
164:     case EPS_CONV_REL:
165:       PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue\n"));break;
166:     case EPS_CONV_NORM:
167:       PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue and matrix norms\n"));
168:       PetscCall(PetscViewerASCIIPrintf(viewer,"  computed matrix norms: norm(A)=%g",(double)eps->nrma));
169:       if (eps->isgeneralized) PetscCall(PetscViewerASCIIPrintf(viewer,", norm(B)=%g",(double)eps->nrmb));
170:       PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
171:       break;
172:     case EPS_CONV_USER:
173:       PetscCall(PetscViewerASCIIPrintf(viewer,"user-defined\n"));break;
174:     }
175:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
176:     if (eps->nini) PetscCall(PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial space: %" PetscInt_FMT "\n",PetscAbs(eps->nini)));
177:     if (eps->ninil) PetscCall(PetscViewerASCIIPrintf(viewer,"  dimension of user-provided left initial space: %" PetscInt_FMT "\n",PetscAbs(eps->ninil)));
178:     if (eps->nds) PetscCall(PetscViewerASCIIPrintf(viewer,"  dimension of user-provided deflation space: %" PetscInt_FMT "\n",PetscAbs(eps->nds)));
179:   } else PetscTryTypeMethod(eps,view,viewer);
180:   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&isexternal,EPSARPACK,EPSBLOPEX,EPSELEMENTAL,EPSFEAST,EPSPRIMME,EPSSCALAPACK,EPSELPA,EPSEVSL,EPSTRLAN,""));
181:   if (!isexternal) {
182:     PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
183:     if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
184:     PetscCall(BVView(eps->V,viewer));
185:     if (eps->rg) {
186:       PetscCall(RGIsTrivial(eps->rg,&istrivial));
187:       if (!istrivial) PetscCall(RGView(eps->rg,viewer));
188:     }
189:     if (eps->useds) {
190:       if (!eps->ds) PetscCall(EPSGetDS(eps,&eps->ds));
191:       PetscCall(DSView(eps->ds,viewer));
192:     }
193:     PetscCall(PetscViewerPopFormat(viewer));
194:   }
195:   if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
196:   PetscCall(STView(eps->st,viewer));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

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

203:    Collective

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

210:    Level: intermediate

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

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

225:    Collective

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

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

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

240:    Level: intermediate

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

250:   PetscFunctionBegin;
251:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)eps));
252:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
253:   if (isAscii) {
254:     PetscCall(PetscViewerGetFormat(viewer,&format));
255:     PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)eps)->tablevel));
256:     PetscCall(EPS_GetActualConverged(eps,&nconv));
257:     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));
258:     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));
259:     PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)eps)->tablevel));
260:   }
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

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

268:    Collective

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

273:    Level: developer

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

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

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

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

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

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

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

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

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

397:    Collective

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

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

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

415:    Level: intermediate

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

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

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

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

455:    Collective

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

460:    Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

652:    Collective

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

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

661:    Level: intermediate

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

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

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

712:    Collective

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

717:    Level: developer

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

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

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

745:    Collective

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

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

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

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

763:    Level: intermediate

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

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

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

804:    Collective

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

809:    Level: developer

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

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