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