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