Actual source code: pepview.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: The PEP routines related to various viewers
12: */
14: #include <slepc/private/pepimpl.h>
15: #include <slepc/private/bvimpl.h>
16: #include <petscdraw.h>
18: /*@
19: PEPView - Prints the `PEP` data structure.
21: Collective
23: Input Parameters:
24: + pep - the polynomial eigensolver context
25: - viewer - optional visualization context
27: Options Database Key:
28: . -pep_view - calls `PEPView()` at end of `PEPSolve()`
30: Notes:
31: The available visualization contexts include
32: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
33: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
34: first process opens the file; all other processes send their data to the
35: first one to print
37: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
38: to output to a specified file.
40: Level: beginner
42: .seealso: [](ch:pep), `PEPCreate()`, `PEPViewFromOptions()`, `STView()`
43: @*/
44: PetscErrorCode PEPView(PEP pep,PetscViewer viewer)
45: {
46: const char *type=NULL;
47: char str[50];
48: PetscBool isascii,islinear,istrivial;
49: PetscInt i;
50: PetscViewer sviewer;
51: MPI_Comm child;
53: PetscFunctionBegin;
55: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer));
57: PetscCheckSameComm(pep,1,viewer,2);
59: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
60: if (isascii) {
61: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pep,viewer));
62: PetscCall(PetscViewerASCIIPushTab(viewer));
63: PetscTryTypeMethod(pep,view,viewer);
64: PetscCall(PetscViewerASCIIPopTab(viewer));
65: if (pep->problem_type) {
66: switch (pep->problem_type) {
67: case PEP_GENERAL: type = "general polynomial eigenvalue problem"; break;
68: case PEP_HERMITIAN: type = SLEPC_STRING_HERMITIAN " polynomial eigenvalue problem"; break;
69: case PEP_HYPERBOLIC: type = "hyperbolic polynomial eigenvalue problem"; break;
70: case PEP_GYROSCOPIC: type = "gyroscopic polynomial eigenvalue problem"; break;
71: }
72: } else type = "not yet set";
73: PetscCall(PetscViewerASCIIPrintf(viewer," problem type: %s\n",type));
74: PetscCall(PetscViewerASCIIPrintf(viewer," polynomial represented in %s basis\n",PEPBasisTypes[pep->basis]));
75: switch (pep->scale) {
76: case PEP_SCALE_NONE:
77: break;
78: case PEP_SCALE_SCALAR:
79: PetscCall(PetscViewerASCIIPrintf(viewer," parameter scaling enabled, with scaling factor=%g\n",(double)pep->sfactor));
80: break;
81: case PEP_SCALE_DIAGONAL:
82: PetscCall(PetscViewerASCIIPrintf(viewer," diagonal balancing enabled, with its=%" PetscInt_FMT " and lambda=%g\n",pep->sits,(double)pep->slambda));
83: break;
84: case PEP_SCALE_BOTH:
85: PetscCall(PetscViewerASCIIPrintf(viewer," parameter scaling & diagonal balancing enabled, with scaling factor=%g, its=%" PetscInt_FMT " and lambda=%g\n",(double)pep->sfactor,pep->sits,(double)pep->slambda));
86: break;
87: }
88: PetscCall(PetscViewerASCIIPrintf(viewer," selected portion of the spectrum: "));
89: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
90: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),pep->target,PETSC_FALSE));
91: if (!pep->which) PetscCall(PetscViewerASCIIPrintf(viewer,"not yet set\n"));
92: else switch (pep->which) {
93: case PEP_WHICH_USER:
94: PetscCall(PetscViewerASCIIPrintf(viewer,"user defined\n"));
95: break;
96: case PEP_TARGET_MAGNITUDE:
97: PetscCall(PetscViewerASCIIPrintf(viewer,"closest to target: %s (in magnitude)\n",str));
98: break;
99: case PEP_TARGET_REAL:
100: PetscCall(PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the real axis)\n",str));
101: break;
102: case PEP_TARGET_IMAGINARY:
103: PetscCall(PetscViewerASCIIPrintf(viewer,"closest to target: %s (along the imaginary axis)\n",str));
104: break;
105: case PEP_LARGEST_MAGNITUDE:
106: PetscCall(PetscViewerASCIIPrintf(viewer,"largest eigenvalues in magnitude\n"));
107: break;
108: case PEP_SMALLEST_MAGNITUDE:
109: PetscCall(PetscViewerASCIIPrintf(viewer,"smallest eigenvalues in magnitude\n"));
110: break;
111: case PEP_LARGEST_REAL:
112: PetscCall(PetscViewerASCIIPrintf(viewer,"largest real parts\n"));
113: break;
114: case PEP_SMALLEST_REAL:
115: PetscCall(PetscViewerASCIIPrintf(viewer,"smallest real parts\n"));
116: break;
117: case PEP_LARGEST_IMAGINARY:
118: PetscCall(PetscViewerASCIIPrintf(viewer,"largest imaginary parts\n"));
119: break;
120: case PEP_SMALLEST_IMAGINARY:
121: PetscCall(PetscViewerASCIIPrintf(viewer,"smallest imaginary parts\n"));
122: break;
123: case PEP_ALL:
124: if (pep->inta || pep->intb) PetscCall(PetscViewerASCIIPrintf(viewer,"all eigenvalues in interval [%g,%g]\n",(double)pep->inta,(double)pep->intb));
125: else PetscCall(PetscViewerASCIIPrintf(viewer,"all eigenvalues in the region\n"));
126: break;
127: }
128: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
129: PetscCall(PetscViewerASCIIPrintf(viewer," number of eigenvalues (nev): %" PetscInt_FMT "\n",pep->nev));
130: PetscCall(PetscViewerASCIIPrintf(viewer," number of column vectors (ncv): %" PetscInt_FMT "\n",pep->ncv));
131: PetscCall(PetscViewerASCIIPrintf(viewer," maximum dimension of projected problem (mpd): %" PetscInt_FMT "\n",pep->mpd));
132: PetscCall(PetscViewerASCIIPrintf(viewer," maximum number of iterations: %" PetscInt_FMT "\n",pep->max_it));
133: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance: %g\n",(double)pep->tol));
134: PetscCall(PetscViewerASCIIPrintf(viewer," convergence test: "));
135: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
136: switch (pep->conv) {
137: case PEP_CONV_ABS:
138: PetscCall(PetscViewerASCIIPrintf(viewer,"absolute\n"));break;
139: case PEP_CONV_REL:
140: PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the eigenvalue\n"));break;
141: case PEP_CONV_NORM:
142: PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the matrix norms\n"));
143: if (pep->nrma) {
144: PetscCall(PetscViewerASCIIPrintf(viewer," computed matrix norms: %g",(double)pep->nrma[0]));
145: for (i=1;i<pep->nmat;i++) PetscCall(PetscViewerASCIIPrintf(viewer,", %g",(double)pep->nrma[i]));
146: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
147: }
148: break;
149: case PEP_CONV_USER:
150: PetscCall(PetscViewerASCIIPrintf(viewer,"user-defined\n"));break;
151: }
152: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
153: PetscCall(PetscViewerASCIIPrintf(viewer," extraction type: %s\n",PEPExtractTypes[pep->extract]));
154: if (pep->refine) {
155: PetscCall(PetscViewerASCIIPrintf(viewer," iterative refinement: %s, with %s scheme\n",PEPRefineTypes[pep->refine],PEPRefineSchemes[pep->scheme]));
156: PetscCall(PetscViewerASCIIPrintf(viewer," refinement stopping criterion: tol=%g, its=%" PetscInt_FMT "\n",(double)pep->rtol,pep->rits));
157: if (pep->npart>1) PetscCall(PetscViewerASCIIPrintf(viewer," splitting communicator in %" PetscInt_FMT " partitions for refinement\n",pep->npart));
158: }
159: if (pep->nini) PetscCall(PetscViewerASCIIPrintf(viewer," dimension of user-provided initial space: %" PetscInt_FMT "\n",PetscAbs(pep->nini)));
160: } else PetscTryTypeMethod(pep,view,viewer);
161: PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
162: if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
163: PetscCall(BVView(pep->V,viewer));
164: if (!pep->rg) PetscCall(PEPGetRG(pep,&pep->rg));
165: PetscCall(RGIsTrivial(pep->rg,&istrivial));
166: if (!istrivial) PetscCall(RGView(pep->rg,viewer));
167: PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear));
168: if (!islinear) {
169: if (!pep->ds) PetscCall(PEPGetDS(pep,&pep->ds));
170: PetscCall(DSView(pep->ds,viewer));
171: }
172: PetscCall(PetscViewerPopFormat(viewer));
173: if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
174: PetscCall(STView(pep->st,viewer));
175: if (pep->refine!=PEP_REFINE_NONE) {
176: PetscCall(PetscViewerASCIIPushTab(viewer));
177: if (pep->npart>1) {
178: if (pep->refinesubc->color==0) {
179: PetscCall(PetscSubcommGetChild(pep->refinesubc,&child));
180: PetscCall(PetscViewerASCIIGetStdout(child,&sviewer));
181: PetscCall(KSPView(pep->refineksp,sviewer));
182: }
183: } else PetscCall(KSPView(pep->refineksp,viewer));
184: PetscCall(PetscViewerASCIIPopTab(viewer));
185: }
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*@
190: PEPViewFromOptions - View (print) a `PEP` object based on values in the options database.
192: Collective
194: Input Parameters:
195: + pep - the polynomial eigensolver context
196: . obj - optional object that provides the options prefix used to query the options database
197: - name - command line option
199: Level: intermediate
201: .seealso: [](ch:pep), `PEPView()`, `PEPCreate()`, `PetscObjectViewFromOptions()`
202: @*/
203: PetscErrorCode PEPViewFromOptions(PEP pep,PetscObject obj,const char name[])
204: {
205: PetscFunctionBegin;
207: PetscCall(PetscObjectViewFromOptions((PetscObject)pep,obj,name));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@
212: PEPConvergedReasonView - Displays the reason a `PEP` solve converged or diverged.
214: Collective
216: Input Parameters:
217: + pep - the polynomial eigensolver context
218: - viewer - the viewer to display the reason
220: Options Database Key:
221: . -pep_converged_reason - print reason for convergence/divergence, and number of iterations
223: Notes:
224: Use `PEPConvergedReasonViewFromOptions()` to display the reason based on values
225: in the options database.
227: To change the format of the output call `PetscViewerPushFormat()` before this
228: call. Use `PETSC_VIEWER_DEFAULT` for the default, or `PETSC_VIEWER_FAILED` to only
229: display a reason if it fails. The latter can be set in the command line with
230: `-pep_converged_reason ::failed`.
232: Level: intermediate
234: .seealso: [](ch:pep), `PEPSetConvergenceTest()`, `PEPSetTolerances()`, `PEPGetIterationNumber()`, `PEPConvergedReasonViewFromOptions()`
235: @*/
236: PetscErrorCode PEPConvergedReasonView(PEP pep,PetscViewer viewer)
237: {
238: PetscBool isAscii;
239: PetscViewerFormat format;
241: PetscFunctionBegin;
242: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pep));
243: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
244: if (isAscii) {
245: PetscCall(PetscViewerGetFormat(viewer,&format));
246: PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)pep)->tablevel));
247: if (pep->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve converged (%" PetscInt_FMT " eigenpair%s) due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",pep->nconv,(pep->nconv>1)?"s":"",PEPConvergedReasons[pep->reason],pep->its));
248: else if (pep->reason <= 0) PetscCall(PetscViewerASCIIPrintf(viewer,"%s Polynomial eigensolve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)pep)->prefix?((PetscObject)pep)->prefix:"",PEPConvergedReasons[pep->reason],pep->its));
249: PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)pep)->tablevel));
250: }
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: PEPConvergedReasonViewFromOptions - Processes command line options to determine if/how
256: the `PEP` converged reason is to be viewed.
258: Collective
260: Input Parameter:
261: . pep - the polynomial eigensolver context
263: Level: intermediate
265: .seealso: [](ch:pep), `PEPConvergedReasonView()`
266: @*/
267: PetscErrorCode PEPConvergedReasonViewFromOptions(PEP pep)
268: {
269: PetscViewer viewer;
270: PetscBool flg;
271: static PetscBool incall = PETSC_FALSE;
272: PetscViewerFormat format;
274: PetscFunctionBegin;
275: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
276: incall = PETSC_TRUE;
277: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_converged_reason",&viewer,&format,&flg));
278: if (flg) {
279: PetscCall(PetscViewerPushFormat(viewer,format));
280: PetscCall(PEPConvergedReasonView(pep,viewer));
281: PetscCall(PetscViewerPopFormat(viewer));
282: PetscCall(PetscViewerDestroy(&viewer));
283: }
284: incall = PETSC_FALSE;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static PetscErrorCode PEPErrorView_ASCII(PEP pep,PEPErrorType etype,PetscViewer viewer)
289: {
290: PetscReal error;
291: PetscInt i,j,k,nvals;
293: PetscFunctionBegin;
294: nvals = (pep->which==PEP_ALL)? pep->nconv: pep->nev;
295: if (pep->which!=PEP_ALL && pep->nconv<pep->nev) {
296: PetscCall(PetscViewerASCIIPrintf(viewer," Problem: less than %" PetscInt_FMT " eigenvalues converged\n\n",pep->nev));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
299: if (pep->which==PEP_ALL && !nvals) {
300: PetscCall(PetscViewerASCIIPrintf(viewer," No eigenvalues have been found\n\n"));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
303: for (i=0;i<nvals;i++) {
304: PetscCall(PEPComputeError(pep,i,etype,&error));
305: if (error>=5.0*pep->tol) {
306: PetscCall(PetscViewerASCIIPrintf(viewer," Problem: some of the first %" PetscInt_FMT " relative errors are higher than the tolerance\n\n",nvals));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
309: }
310: if (pep->which==PEP_ALL) PetscCall(PetscViewerASCIIPrintf(viewer," Found %" PetscInt_FMT " eigenvalues, all of them computed up to the required tolerance:",nvals));
311: else PetscCall(PetscViewerASCIIPrintf(viewer," All requested eigenvalues computed up to the required tolerance:"));
312: for (i=0;i<=(nvals-1)/8;i++) {
313: PetscCall(PetscViewerASCIIPrintf(viewer,"\n "));
314: for (j=0;j<PetscMin(8,nvals-8*i);j++) {
315: k = pep->perm[8*i+j];
316: PetscCall(SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]));
317: if (8*i+j+1<nvals) PetscCall(PetscViewerASCIIPrintf(viewer,", "));
318: }
319: }
320: PetscCall(PetscViewerASCIIPrintf(viewer,"\n\n"));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode PEPErrorView_DETAIL(PEP pep,PEPErrorType etype,PetscViewer viewer)
325: {
326: PetscReal error,re,im;
327: PetscScalar kr,ki;
328: PetscInt i;
329: char ex[30],sep[]=" ---------------------- --------------------\n";
331: PetscFunctionBegin;
332: if (!pep->nconv) PetscFunctionReturn(PETSC_SUCCESS);
333: switch (etype) {
334: case PEP_ERROR_ABSOLUTE:
335: PetscCall(PetscSNPrintf(ex,sizeof(ex)," ||P(k)x||"));
336: break;
337: case PEP_ERROR_RELATIVE:
338: PetscCall(PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||"));
339: break;
340: case PEP_ERROR_BACKWARD:
341: PetscCall(PetscSNPrintf(ex,sizeof(ex)," eta(x,k)"));
342: break;
343: }
344: PetscCall(PetscViewerASCIIPrintf(viewer,"%s k %s\n%s",sep,ex,sep));
345: for (i=0;i<pep->nconv;i++) {
346: PetscCall(PEPGetEigenpair(pep,i,&kr,&ki,NULL,NULL));
347: PetscCall(PEPComputeError(pep,i,etype,&error));
348: #if defined(PETSC_USE_COMPLEX)
349: re = PetscRealPart(kr);
350: im = PetscImaginaryPart(kr);
351: #else
352: re = kr;
353: im = ki;
354: #endif
355: if (im!=0.0) PetscCall(PetscViewerASCIIPrintf(viewer," % 9f%+9fi %12g\n",(double)re,(double)im,(double)error));
356: else PetscCall(PetscViewerASCIIPrintf(viewer," % 12f %12g\n",(double)re,(double)error));
357: }
358: PetscCall(PetscViewerASCIIPrintf(viewer,"%s",sep));
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode PEPErrorView_MATLAB(PEP pep,PEPErrorType etype,PetscViewer viewer)
363: {
364: PetscReal error;
365: PetscInt i;
366: const char *name;
368: PetscFunctionBegin;
369: PetscCall(PetscObjectGetName((PetscObject)pep,&name));
370: PetscCall(PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name));
371: for (i=0;i<pep->nconv;i++) {
372: PetscCall(PEPComputeError(pep,i,etype,&error));
373: PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error));
374: }
375: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@
380: PEPErrorView - Displays the errors associated with the computed solution
381: (as well as the eigenvalues).
383: Collective
385: Input Parameters:
386: + pep - the polynomial eigensolver context
387: . etype - error type
388: - viewer - optional visualization context
390: Options Database Keys:
391: + -pep_error_absolute - print absolute errors of each eigenpair
392: . -pep_error_relative - print relative errors of each eigenpair
393: - -pep_error_backward - print backward errors of each eigenpair
395: Notes:
396: By default, this function checks the error of all eigenpairs and prints
397: the eigenvalues if all of them are below the requested tolerance.
398: If the viewer has format `PETSC_VIEWER_ASCII_INFO_DETAIL` then a table with
399: eigenvalues and corresponding errors is printed.
401: All the command-line options listed above admit an optional argument
402: specifying the viewer type and options. For instance, use
403: `-pep_error_relative :myerr.m:ascii_matlab` to save the errors in a file
404: that can be executed in Matlab.
406: Level: intermediate
408: .seealso: [](ch:pep), `PEPSolve()`, `PEPValuesView()`, `PEPVectorsView()`
409: @*/
410: PetscErrorCode PEPErrorView(PEP pep,PEPErrorType etype,PetscViewer viewer)
411: {
412: PetscBool isascii;
413: PetscViewerFormat format;
415: PetscFunctionBegin;
417: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer));
419: PetscCheckSameComm(pep,1,viewer,3);
420: PEPCheckSolved(pep,1);
421: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
422: if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
424: PetscCall(PetscViewerGetFormat(viewer,&format));
425: switch (format) {
426: case PETSC_VIEWER_DEFAULT:
427: case PETSC_VIEWER_ASCII_INFO:
428: PetscCall(PEPErrorView_ASCII(pep,etype,viewer));
429: break;
430: case PETSC_VIEWER_ASCII_INFO_DETAIL:
431: PetscCall(PEPErrorView_DETAIL(pep,etype,viewer));
432: break;
433: case PETSC_VIEWER_ASCII_MATLAB:
434: PetscCall(PEPErrorView_MATLAB(pep,etype,viewer));
435: break;
436: default:
437: PetscCall(PetscInfo(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
438: }
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: PEPErrorViewFromOptions - Processes command line options to determine if/how
444: the errors of the computed solution are to be viewed.
446: Collective
448: Input Parameter:
449: . pep - the polynomial eigensolver context
451: Level: developer
453: .seealso: [](ch:pep), `PEPErrorView()`
454: @*/
455: PetscErrorCode PEPErrorViewFromOptions(PEP pep)
456: {
457: PetscViewer viewer;
458: PetscBool flg;
459: static PetscBool incall = PETSC_FALSE;
460: PetscViewerFormat format;
462: PetscFunctionBegin;
463: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
464: incall = PETSC_TRUE;
465: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_absolute",&viewer,&format,&flg));
466: if (flg) {
467: PetscCall(PetscViewerPushFormat(viewer,format));
468: PetscCall(PEPErrorView(pep,PEP_ERROR_ABSOLUTE,viewer));
469: PetscCall(PetscViewerPopFormat(viewer));
470: PetscCall(PetscViewerDestroy(&viewer));
471: }
472: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_relative",&viewer,&format,&flg));
473: if (flg) {
474: PetscCall(PetscViewerPushFormat(viewer,format));
475: PetscCall(PEPErrorView(pep,PEP_ERROR_RELATIVE,viewer));
476: PetscCall(PetscViewerPopFormat(viewer));
477: PetscCall(PetscViewerDestroy(&viewer));
478: }
479: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_error_backward",&viewer,&format,&flg));
480: if (flg) {
481: PetscCall(PetscViewerPushFormat(viewer,format));
482: PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,viewer));
483: PetscCall(PetscViewerPopFormat(viewer));
484: PetscCall(PetscViewerDestroy(&viewer));
485: }
486: incall = PETSC_FALSE;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode PEPValuesView_DRAW(PEP pep,PetscViewer viewer)
491: {
492: PetscDraw draw;
493: PetscDrawSP drawsp;
494: PetscReal re,im;
495: PetscInt i,k;
497: PetscFunctionBegin;
498: if (!pep->nconv) PetscFunctionReturn(PETSC_SUCCESS);
499: PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
500: PetscCall(PetscDrawSetTitle(draw,"Computed Eigenvalues"));
501: PetscCall(PetscDrawSPCreate(draw,1,&drawsp));
502: for (i=0;i<pep->nconv;i++) {
503: k = pep->perm[i];
504: #if defined(PETSC_USE_COMPLEX)
505: re = PetscRealPart(pep->eigr[k]);
506: im = PetscImaginaryPart(pep->eigr[k]);
507: #else
508: re = pep->eigr[k];
509: im = pep->eigi[k];
510: #endif
511: PetscCall(PetscDrawSPAddPoint(drawsp,&re,&im));
512: }
513: PetscCall(PetscDrawSPDraw(drawsp,PETSC_TRUE));
514: PetscCall(PetscDrawSPSave(drawsp));
515: PetscCall(PetscDrawSPDestroy(&drawsp));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode PEPValuesView_BINARY(PEP pep,PetscViewer viewer)
520: {
521: #if defined(PETSC_HAVE_COMPLEX)
522: PetscInt i,k;
523: PetscComplex *ev;
524: #endif
526: PetscFunctionBegin;
527: #if defined(PETSC_HAVE_COMPLEX)
528: PetscCall(PetscMalloc1(pep->nconv,&ev));
529: for (i=0;i<pep->nconv;i++) {
530: k = pep->perm[i];
531: #if defined(PETSC_USE_COMPLEX)
532: ev[i] = pep->eigr[k];
533: #else
534: ev[i] = PetscCMPLX(pep->eigr[k],pep->eigi[k]);
535: #endif
536: }
537: PetscCall(PetscViewerBinaryWrite(viewer,ev,pep->nconv,PETSC_COMPLEX));
538: PetscCall(PetscFree(ev));
539: #endif
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: #if defined(PETSC_HAVE_HDF5)
544: static PetscErrorCode PEPValuesView_HDF5(PEP pep,PetscViewer viewer)
545: {
546: PetscInt i,k,n,N;
547: PetscMPIInt rank;
548: Vec v;
549: char vname[30];
550: const char *ename;
552: PetscFunctionBegin;
553: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rank));
554: N = pep->nconv;
555: n = rank? 0: N;
556: /* create a vector containing the eigenvalues */
557: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),n,N,&v));
558: PetscCall(PetscObjectGetName((PetscObject)pep,&ename));
559: PetscCall(PetscSNPrintf(vname,sizeof(vname),"eigr_%s",ename));
560: PetscCall(PetscObjectSetName((PetscObject)v,vname));
561: if (!rank) {
562: for (i=0;i<pep->nconv;i++) {
563: k = pep->perm[i];
564: PetscCall(VecSetValue(v,i,pep->eigr[k],INSERT_VALUES));
565: }
566: }
567: PetscCall(VecAssemblyBegin(v));
568: PetscCall(VecAssemblyEnd(v));
569: PetscCall(VecView(v,viewer));
570: #if !defined(PETSC_USE_COMPLEX)
571: /* in real scalars write the imaginary part as a separate vector */
572: PetscCall(PetscSNPrintf(vname,sizeof(vname),"eigi_%s",ename));
573: PetscCall(PetscObjectSetName((PetscObject)v,vname));
574: if (!rank) {
575: for (i=0;i<pep->nconv;i++) {
576: k = pep->perm[i];
577: PetscCall(VecSetValue(v,i,pep->eigi[k],INSERT_VALUES));
578: }
579: }
580: PetscCall(VecAssemblyBegin(v));
581: PetscCall(VecAssemblyEnd(v));
582: PetscCall(VecView(v,viewer));
583: #endif
584: PetscCall(VecDestroy(&v));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
587: #endif
589: static PetscErrorCode PEPValuesView_ASCII(PEP pep,PetscViewer viewer)
590: {
591: PetscInt i,k;
593: PetscFunctionBegin;
594: PetscCall(PetscViewerASCIIPrintf(viewer,"Eigenvalues = \n"));
595: for (i=0;i<pep->nconv;i++) {
596: k = pep->perm[i];
597: PetscCall(PetscViewerASCIIPrintf(viewer," "));
598: PetscCall(SlepcPrintEigenvalueASCII(viewer,pep->eigr[k],pep->eigi[k]));
599: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
600: }
601: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode PEPValuesView_MATLAB(PEP pep,PetscViewer viewer)
606: {
607: PetscInt i,k;
608: PetscReal re,im;
609: const char *name;
611: PetscFunctionBegin;
612: PetscCall(PetscObjectGetName((PetscObject)pep,&name));
613: PetscCall(PetscViewerASCIIPrintf(viewer,"Lambda_%s = [\n",name));
614: for (i=0;i<pep->nconv;i++) {
615: k = pep->perm[i];
616: #if defined(PETSC_USE_COMPLEX)
617: re = PetscRealPart(pep->eigr[k]);
618: im = PetscImaginaryPart(pep->eigr[k]);
619: #else
620: re = pep->eigr[k];
621: im = pep->eigi[k];
622: #endif
623: if (im!=0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei\n",(double)re,(double)im));
624: else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)re));
625: }
626: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: /*@
631: PEPValuesView - Displays the computed eigenvalues in a viewer.
633: Collective
635: Input Parameters:
636: + pep - the polynomial eigensolver context
637: - viewer - the viewer
639: Options Database Key:
640: . -pep_view_values - print computed eigenvalues
642: Note:
643: The command-line option listed above admits an optional argument
644: specifying the viewer type and options. For instance, use
645: `-pep_view_values :evals.m:ascii_matlab` to save the values in a file
646: that can be executed in Matlab.
648: Level: intermediate
650: .seealso: [](ch:pep), `PEPSolve()`, `PEPVectorsView()`, `PEPErrorView()`
651: @*/
652: PetscErrorCode PEPValuesView(PEP pep,PetscViewer viewer)
653: {
654: PetscBool isascii,isdraw,isbinary;
655: PetscViewerFormat format;
656: #if defined(PETSC_HAVE_HDF5)
657: PetscBool ishdf5;
658: #endif
660: PetscFunctionBegin;
662: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer));
664: PetscCheckSameComm(pep,1,viewer,2);
665: PEPCheckSolved(pep,1);
666: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
667: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary));
668: #if defined(PETSC_HAVE_HDF5)
669: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5));
670: #endif
671: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
672: if (isdraw) PetscCall(PEPValuesView_DRAW(pep,viewer));
673: else if (isbinary) PetscCall(PEPValuesView_BINARY(pep,viewer));
674: #if defined(PETSC_HAVE_HDF5)
675: else if (ishdf5) PetscCall(PEPValuesView_HDF5(pep,viewer));
676: #endif
677: else if (isascii) {
678: PetscCall(PetscViewerGetFormat(viewer,&format));
679: switch (format) {
680: case PETSC_VIEWER_DEFAULT:
681: case PETSC_VIEWER_ASCII_INFO:
682: case PETSC_VIEWER_ASCII_INFO_DETAIL:
683: PetscCall(PEPValuesView_ASCII(pep,viewer));
684: break;
685: case PETSC_VIEWER_ASCII_MATLAB:
686: PetscCall(PEPValuesView_MATLAB(pep,viewer));
687: break;
688: default:
689: PetscCall(PetscInfo(pep,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
690: }
691: }
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: PEPValuesViewFromOptions - Processes command line options to determine if/how
697: the computed eigenvalues are to be viewed.
699: Collective
701: Input Parameter:
702: . pep - the polynomial eigensolver context
704: Level: developer
706: .seealso: [](ch:pep), `PEPValuesView()`
707: @*/
708: PetscErrorCode PEPValuesViewFromOptions(PEP pep)
709: {
710: PetscViewer viewer;
711: PetscBool flg;
712: static PetscBool incall = PETSC_FALSE;
713: PetscViewerFormat format;
715: PetscFunctionBegin;
716: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
717: incall = PETSC_TRUE;
718: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_values",&viewer,&format,&flg));
719: if (flg) {
720: PetscCall(PetscViewerPushFormat(viewer,format));
721: PetscCall(PEPValuesView(pep,viewer));
722: PetscCall(PetscViewerPopFormat(viewer));
723: PetscCall(PetscViewerDestroy(&viewer));
724: }
725: incall = PETSC_FALSE;
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*@
730: PEPVectorsView - Outputs computed eigenvectors to a viewer.
732: Collective
734: Input Parameters:
735: + pep - the polynomial eigensolver context
736: - viewer - the viewer
738: Options Database Key:
739: . -pep_view_vectors - output eigenvectors
741: Notes:
742: If PETSc was configured with real scalars, complex conjugate eigenvectors
743: will be viewed as two separate real vectors, one containing the real part
744: and another one containing the imaginary part.
746: The command-line option listed above admits an optional argument
747: specifying the viewer type and options. For instance, use
748: `-pep_view_vectors binary:evecs.bin` to save the vectors in a binary file.
750: Level: intermediate
752: .seealso: [](ch:pep), `PEPSolve()`, `PEPValuesView()`, `PEPErrorView()`
753: @*/
754: PetscErrorCode PEPVectorsView(PEP pep,PetscViewer viewer)
755: {
756: PetscInt i,k;
757: Vec xr,xi=NULL;
759: PetscFunctionBegin;
761: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pep),&viewer));
763: PetscCheckSameComm(pep,1,viewer,2);
764: PEPCheckSolved(pep,1);
765: if (pep->nconv) {
766: PetscCall(PEPComputeVectors(pep));
767: PetscCall(BVCreateVec(pep->V,&xr));
768: #if !defined(PETSC_USE_COMPLEX)
769: PetscCall(BVCreateVec(pep->V,&xi));
770: #endif
771: for (i=0;i<pep->nconv;i++) {
772: k = pep->perm[i];
773: PetscCall(BV_GetEigenvector(pep->V,k,pep->eigi[k],xr,xi));
774: PetscCall(SlepcViewEigenvector(viewer,xr,xi,"X",i,(PetscObject)pep));
775: }
776: PetscCall(VecDestroy(&xr));
777: #if !defined(PETSC_USE_COMPLEX)
778: PetscCall(VecDestroy(&xi));
779: #endif
780: }
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: /*@
785: PEPVectorsViewFromOptions - Processes command line options to determine if/how
786: the computed eigenvectors are to be viewed.
788: Collective
790: Input Parameter:
791: . pep - the polynomial eigensolver context
793: Level: developer
795: .seealso: [](ch:pep), `PEPVectorsView()`
796: @*/
797: PetscErrorCode PEPVectorsViewFromOptions(PEP pep)
798: {
799: PetscViewer viewer;
800: PetscBool flg = PETSC_FALSE;
801: static PetscBool incall = PETSC_FALSE;
802: PetscViewerFormat format;
804: PetscFunctionBegin;
805: if (incall) PetscFunctionReturn(PETSC_SUCCESS);
806: incall = PETSC_TRUE;
807: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)pep),((PetscObject)pep)->options,((PetscObject)pep)->prefix,"-pep_view_vectors",&viewer,&format,&flg));
808: if (flg) {
809: PetscCall(PetscViewerPushFormat(viewer,format));
810: PetscCall(PEPVectorsView(pep,viewer));
811: PetscCall(PetscViewerPopFormat(viewer));
812: PetscCall(PetscViewerDestroy(&viewer));
813: }
814: incall = PETSC_FALSE;
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }