Actual source code: epsview.c

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

203: /*@
204:    EPSViewFromOptions - View from options

206:    Collective

208:    Input Parameters:
209: +  eps  - the eigensolver context
210: .  obj  - optional object
211: -  name - command line option

213:    Level: intermediate

215: .seealso: `EPSView()`, `EPSCreate()`
216: @*/
217: PetscErrorCode EPSViewFromOptions(EPS eps,PetscObject obj,const char name[])
218: {
219:   PetscFunctionBegin;
221:   PetscCall(PetscObjectViewFromOptions((PetscObject)eps,obj,name));
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

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

228:    Collective

230:    Input Parameters:
231: +  eps - the eigensolver context
232: -  viewer - the viewer to display the reason

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

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

243:    Level: intermediate

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

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

267: /*@
268:    EPSConvergedReasonViewFromOptions - Processes command line options to determine if/how
269:    the EPS converged reason is to be viewed.

271:    Collective

273:    Input Parameter:
274: .  eps - the eigensolver context

276:    Level: developer

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

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

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

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

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

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

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

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

396: /*@
397:    EPSErrorView - Displays the errors associated with the computed solution
398:    (as well as the eigenvalues).

400:    Collective

402:    Input Parameters:
403: +  eps    - the eigensolver context
404: .  etype  - error type
405: -  viewer - optional visualization context

407:    Options Database Keys:
408: +  -eps_error_absolute - print absolute errors of each eigenpair
409: .  -eps_error_relative - print relative errors of each eigenpair
410: -  -eps_error_backward - print backward errors of each eigenpair

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

418:    Level: intermediate

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

427:   PetscFunctionBegin;
429:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)eps),&viewer));
431:   PetscCheckSameComm(eps,1,viewer,3);
432:   EPSCheckSolved(eps,1);
433:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
434:   if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);

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

454: /*@
455:    EPSErrorViewFromOptions - Processes command line options to determine if/how
456:    the errors of the computed solution are to be viewed.

458:    Collective

460:    Input Parameter:
461: .  eps - the eigensolver context

463:    Level: developer

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

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

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

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

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

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

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

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

607: static PetscErrorCode EPSValuesView_ASCII(EPS eps,PetscViewer viewer)
608: {
609:   PetscInt       i,nconv;
610:   PetscScalar    kr,ki;

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

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

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

652: /*@
653:    EPSValuesView - Displays the computed eigenvalues in a viewer.

655:    Collective

657:    Input Parameters:
658: +  eps    - the eigensolver context
659: -  viewer - the viewer

661:    Options Database Key:
662: .  -eps_view_values - print computed eigenvalues

664:    Level: intermediate

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

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

711: /*@
712:    EPSValuesViewFromOptions - Processes command line options to determine if/how
713:    the computed eigenvalues are to be viewed.

715:    Collective

717:    Input Parameters:
718: .  eps - the eigensolver context

720:    Level: developer

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

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

745: /*@
746:    EPSVectorsView - Outputs computed eigenvectors to a viewer.

748:    Collective

750:    Input Parameters:
751: +  eps    - the eigensolver context
752: -  viewer - the viewer

754:    Options Database Key:
755: .  -eps_view_vectors - output eigenvectors.

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

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

766:    Level: intermediate

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

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

803: /*@
804:    EPSVectorsViewFromOptions - Processes command line options to determine if/how
805:    the computed eigenvectors are to be viewed.

807:    Collective

809:    Input Parameter:
810: .  eps - the eigensolver context

812:    Level: developer

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

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