Actual source code: svdview.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 SVD routines related to various viewers
 12: */

 14: #include <slepc/private/svdimpl.h>      /*I "slepcsvd.h" I*/
 15: #include <petscdraw.h>

 17: /*@C
 18:    SVDView - Prints the SVD data structure.

 20:    Collective on SVD

 22:    Input Parameters:
 23: +  svd - the singular value solver context
 24: -  viewer - optional visualization context

 26:    Options Database Key:
 27: .  -svd_view -  Calls SVDView() at end of SVDSolve()

 29:    Note:
 30:    The available visualization contexts include
 31: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 32: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 33:          output where only the first processor opens
 34:          the file.  All other processors send their
 35:          data to the first processor to print.

 37:    The user can open an alternative visualization context with
 38:    PetscViewerASCIIOpen() - output to a specified file.

 40:    Level: beginner

 42: .seealso: STView(), PetscViewerASCIIOpen()
 43: @*/
 44: PetscErrorCode SVDView(SVD svd,PetscViewer viewer)
 45: {
 47:   PetscBool      isascii,isshell;

 51:   if (!viewer) {
 52:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer);
 53:   }

 57:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 58:   if (isascii) {
 59:     PetscObjectPrintClassNamePrefixType((PetscObject)svd,viewer);
 60:     if (svd->ops->view) {
 61:       PetscViewerASCIIPushTab(viewer);
 62:       (*svd->ops->view)(svd,viewer);
 63:       PetscViewerASCIIPopTab(viewer);
 64:     }
 65:     PetscViewerASCIIPrintf(viewer,"  transpose mode: %s\n",svd->impltrans?"implicit":"explicit");
 66:     if (svd->which == SVD_LARGEST) {
 67:       PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: largest\n");
 68:     } else {
 69:       PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: smallest\n");
 70:     }
 71:     PetscViewerASCIIPrintf(viewer,"  number of singular values (nsv): %D\n",svd->nsv);
 72:     PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %D\n",svd->ncv);
 73:     PetscViewerASCIIPrintf(viewer,"  maximum dimension of projected problem (mpd): %D\n",svd->mpd);
 74:     PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %D\n",svd->max_it);
 75:     PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)svd->tol);
 76:     PetscViewerASCIIPrintf(viewer,"  convergence test: ");
 77:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 78:     switch (svd->conv) {
 79:     case SVD_CONV_ABS:
 80:       PetscViewerASCIIPrintf(viewer,"absolute\n");break;
 81:     case SVD_CONV_REL:
 82:       PetscViewerASCIIPrintf(viewer,"relative to the singular value\n");break;
 83:     case SVD_CONV_USER:
 84:       PetscViewerASCIIPrintf(viewer,"user-defined\n");break;
 85:     }
 86:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
 87:     if (svd->nini) {
 88:       PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial space: %D\n",PetscAbs(svd->nini));
 89:     }
 90:     if (svd->ninil) {
 91:       PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial left space: %D\n",PetscAbs(svd->ninil));
 92:     }
 93:   } else {
 94:     if (svd->ops->view) {
 95:       (*svd->ops->view)(svd,viewer);
 96:     }
 97:   }
 98:   PetscObjectTypeCompareAny((PetscObject)svd,&isshell,SVDCROSS,SVDCYCLIC,SVDPRIMME,"");
 99:   if (!isshell) {
100:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
101:     if (!svd->V) { SVDGetBV(svd,&svd->V,NULL); }
102:     BVView(svd->V,viewer);
103:     if (!svd->ds) { SVDGetDS(svd,&svd->ds); }
104:     DSView(svd->ds,viewer);
105:     PetscViewerPopFormat(viewer);
106:   }
107:   return(0);
108: }

110: /*@C
111:    SVDReasonView - Displays the reason an SVD solve converged or diverged.

113:    Collective on SVD

115:    Parameter:
116: +  svd - the singular value solver context
117: -  viewer - the viewer to display the reason

119:    Options Database Keys:
120: .  -svd_converged_reason - print reason for convergence, and number of iterations

122:    Level: intermediate

124: .seealso: SVDSetTolerances(), SVDGetIterationNumber()
125: @*/
126: PetscErrorCode SVDReasonView(SVD svd,PetscViewer viewer)
127: {
129:   PetscBool      isAscii;

132:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
133:   if (isAscii) {
134:     PetscViewerASCIIAddTab(viewer,((PetscObject)svd)->tablevel);
135:     if (svd->reason > 0) {
136:       PetscViewerASCIIPrintf(viewer,"%s SVD solve converged (%D singular triplet%s) due to %s; iterations %D\n",((PetscObject)svd)->prefix?((PetscObject)svd)->prefix:"",svd->nconv,(svd->nconv>1)?"s":"",SVDConvergedReasons[svd->reason],svd->its);
137:     } else {
138:       PetscViewerASCIIPrintf(viewer,"%s SVD solve did not converge due to %s; iterations %D\n",((PetscObject)svd)->prefix?((PetscObject)svd)->prefix:"",SVDConvergedReasons[svd->reason],svd->its);
139:     }
140:     PetscViewerASCIISubtractTab(viewer,((PetscObject)svd)->tablevel);
141:   }
142:   return(0);
143: }

145: /*@
146:    SVDReasonViewFromOptions - Processes command line options to determine if/how
147:    the SVD converged reason is to be viewed.

149:    Collective on SVD

151:    Input Parameters:
152: .  svd - the singular value solver context

154:    Level: developer
155: @*/
156: PetscErrorCode SVDReasonViewFromOptions(SVD svd)
157: {
158:   PetscErrorCode    ierr;
159:   PetscViewer       viewer;
160:   PetscBool         flg;
161:   static PetscBool  incall = PETSC_FALSE;
162:   PetscViewerFormat format;

165:   if (incall) return(0);
166:   incall = PETSC_TRUE;
167:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_converged_reason",&viewer,&format,&flg);
168:   if (flg) {
169:     PetscViewerPushFormat(viewer,format);
170:     SVDReasonView(svd,viewer);
171:     PetscViewerPopFormat(viewer);
172:     PetscViewerDestroy(&viewer);
173:   }
174:   incall = PETSC_FALSE;
175:   return(0);
176: }

178: static PetscErrorCode SVDErrorView_ASCII(SVD svd,SVDErrorType etype,PetscViewer viewer)
179: {
180:   PetscReal      error,sigma;
181:   PetscInt       i,j;

185:   if (svd->nconv<svd->nsv) {
186:     PetscViewerASCIIPrintf(viewer," Problem: less than %D singular values converged\n\n",svd->nsv);
187:     return(0);
188:   }
189:   for (i=0;i<svd->nsv;i++) {
190:     SVDComputeError(svd,i,etype,&error);
191:     if (error>=5.0*svd->tol) {
192:       PetscViewerASCIIPrintf(viewer," Problem: some of the first %D relative errors are higher than the tolerance\n\n",svd->nsv);
193:       return(0);
194:     }
195:   }
196:   PetscViewerASCIIPrintf(viewer," All requested singular values computed up to the required tolerance:");
197:   for (i=0;i<=(svd->nsv-1)/8;i++) {
198:     PetscViewerASCIIPrintf(viewer,"\n     ");
199:     for (j=0;j<PetscMin(8,svd->nsv-8*i);j++) {
200:       SVDGetSingularTriplet(svd,8*i+j,&sigma,NULL,NULL);
201:       PetscViewerASCIIPrintf(viewer,"%.5f",(double)sigma);
202:       if (8*i+j+1<svd->nsv) { PetscViewerASCIIPrintf(viewer,", "); }
203:     }
204:   }
205:   PetscViewerASCIIPrintf(viewer,"\n\n");
206:   return(0);
207: }

209: static PetscErrorCode SVDErrorView_DETAIL(SVD svd,SVDErrorType etype,PetscViewer viewer)
210: {
212:   PetscReal      error,sigma;
213:   PetscInt       i;
214: #define EXLEN 30
215:   char           ex[EXLEN],sep[]=" ---------------------- --------------------\n";

218:   if (!svd->nconv) return(0);
219:   switch (etype) {
220:     case SVD_ERROR_ABSOLUTE:
221:       PetscSNPrintf(ex,EXLEN," absolute error");
222:       break;
223:     case SVD_ERROR_RELATIVE:
224:       PetscSNPrintf(ex,EXLEN," relative error");
225:       break;
226:   }
227:   PetscViewerASCIIPrintf(viewer,"%s          sigma           %s\n%s",sep,ex,sep);
228:   for (i=0;i<svd->nconv;i++) {
229:     SVDGetSingularTriplet(svd,i,&sigma,NULL,NULL);
230:     SVDComputeError(svd,i,etype,&error);
231:     PetscViewerASCIIPrintf(viewer,"       % 6f          %12g\n",(double)sigma,(double)error);
232:   }
233:   PetscViewerASCIIPrintf(viewer,sep);
234:   return(0);
235: }

237: static PetscErrorCode SVDErrorView_MATLAB(SVD svd,SVDErrorType etype,PetscViewer viewer)
238: {
240:   PetscReal      error;
241:   PetscInt       i;
242:   const char     *name;

245:   PetscObjectGetName((PetscObject)svd,&name);
246:   PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name);
247:   for (i=0;i<svd->nconv;i++) {
248:     SVDComputeError(svd,i,etype,&error);
249:     PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error);
250:   }
251:   PetscViewerASCIIPrintf(viewer,"];\n");
252:   return(0);
253: }

255: /*@C
256:    SVDErrorView - Displays the errors associated with the computed solution
257:    (as well as the singular values).

259:    Collective on SVD

261:    Input Parameters:
262: +  svd    - the singular value solver context
263: .  etype  - error type
264: -  viewer - optional visualization context

266:    Options Database Key:
267: +  -svd_error_absolute - print absolute errors of each singular triplet
268: -  -svd_error_relative - print relative errors of each singular triplet

270:    Notes:
271:    By default, this function checks the error of all singular triplets and prints
272:    the singular values if all of them are below the requested tolerance.
273:    If the viewer has format=PETSC_VIEWER_ASCII_INFO_DETAIL then a table with
274:    singular values and corresponding errors is printed.

276:    Level: intermediate

278: .seealso: SVDSolve(), SVDValuesView(), SVDVectorsView()
279: @*/
280: PetscErrorCode SVDErrorView(SVD svd,SVDErrorType etype,PetscViewer viewer)
281: {
282:   PetscBool         isascii;
283:   PetscViewerFormat format;
284:   PetscErrorCode    ierr;

288:   if (!viewer) {
289:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer);
290:   }
293:   SVDCheckSolved(svd,1);
294:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
295:   if (!isascii) return(0);

297:   PetscViewerGetFormat(viewer,&format);
298:   switch (format) {
299:     case PETSC_VIEWER_DEFAULT:
300:     case PETSC_VIEWER_ASCII_INFO:
301:       SVDErrorView_ASCII(svd,etype,viewer);
302:       break;
303:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
304:       SVDErrorView_DETAIL(svd,etype,viewer);
305:       break;
306:     case PETSC_VIEWER_ASCII_MATLAB:
307:       SVDErrorView_MATLAB(svd,etype,viewer);
308:       break;
309:     default:
310:       PetscInfo1(svd,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
311:   }
312:   return(0);
313: }

315: /*@
316:    SVDErrorViewFromOptions - Processes command line options to determine if/how
317:    the errors of the computed solution are to be viewed.

319:    Collective on SVD

321:    Input Parameters:
322: .  svd - the singular value solver context

324:    Level: developer
325: @*/
326: PetscErrorCode SVDErrorViewFromOptions(SVD svd)
327: {
328:   PetscErrorCode    ierr;
329:   PetscViewer       viewer;
330:   PetscBool         flg;
331:   static PetscBool  incall = PETSC_FALSE;
332:   PetscViewerFormat format;

335:   if (incall) return(0);
336:   incall = PETSC_TRUE;
337:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_error_absolute",&viewer,&format,&flg);
338:   if (flg) {
339:     PetscViewerPushFormat(viewer,format);
340:     SVDErrorView(svd,SVD_ERROR_ABSOLUTE,viewer);
341:     PetscViewerPopFormat(viewer);
342:     PetscViewerDestroy(&viewer);
343:   }
344:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_error_relative",&viewer,&format,&flg);
345:   if (flg) {
346:     PetscViewerPushFormat(viewer,format);
347:     SVDErrorView(svd,SVD_ERROR_RELATIVE,viewer);
348:     PetscViewerPopFormat(viewer);
349:     PetscViewerDestroy(&viewer);
350:   }
351:   incall = PETSC_FALSE;
352:   return(0);
353: }

355: static PetscErrorCode SVDValuesView_DRAW(SVD svd,PetscViewer viewer)
356: {
358:   PetscDraw      draw;
359:   PetscDrawSP    drawsp;
360:   PetscReal      re,im=0.0;
361:   PetscInt       i;

364:   if (!svd->nconv) return(0);
365:   PetscViewerDrawGetDraw(viewer,0,&draw);
366:   PetscDrawSetTitle(draw,"Computed singular values");
367:   PetscDrawSPCreate(draw,1,&drawsp);
368:   for (i=0;i<svd->nconv;i++) {
369:     re = svd->sigma[svd->perm[i]];
370:     PetscDrawSPAddPoint(drawsp,&re,&im);
371:   }
372:   PetscDrawSPDraw(drawsp,PETSC_TRUE);
373:   PetscDrawSPSave(drawsp);
374:   PetscDrawSPDestroy(&drawsp);
375:   return(0);
376: }

378: static PetscErrorCode SVDValuesView_ASCII(SVD svd,PetscViewer viewer)
379: {
380:   PetscInt       i;

384:   PetscViewerASCIIPrintf(viewer,"Singular values = \n");
385:   for (i=0;i<svd->nconv;i++) {
386:     PetscViewerASCIIPrintf(viewer,"   %.5f\n",(double)svd->sigma[svd->perm[i]]);
387:   }
388:   PetscViewerASCIIPrintf(viewer,"\n");
389:   return(0);
390: }

392: static PetscErrorCode SVDValuesView_MATLAB(SVD svd,PetscViewer viewer)
393: {
395:   PetscInt       i;
396:   const char     *name;

399:   PetscObjectGetName((PetscObject)svd,&name);
400:   PetscViewerASCIIPrintf(viewer,"Sigma_%s = [\n",name);
401:   for (i=0;i<svd->nconv;i++) {
402:     PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)svd->sigma[svd->perm[i]]);
403:   }
404:   PetscViewerASCIIPrintf(viewer,"];\n");
405:   return(0);
406: }

408: /*@C
409:    SVDValuesView - Displays the computed singular values in a viewer.

411:    Collective on SVD

413:    Input Parameters:
414: +  svd    - the singular value solver context
415: -  viewer - the viewer

417:    Options Database Key:
418: .  -svd_view_values - print computed singular values

420:    Level: intermediate

422: .seealso: SVDSolve(), SVDVectorsView(), SVDErrorView()
423: @*/
424: PetscErrorCode SVDValuesView(SVD svd,PetscViewer viewer)
425: {
426:   PetscBool         isascii,isdraw;
427:   PetscViewerFormat format;
428:   PetscErrorCode    ierr;

432:   if (!viewer) {
433:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer);
434:   }
437:   SVDCheckSolved(svd,1);
438:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
439:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
440:   if (isdraw) {
441:     SVDValuesView_DRAW(svd,viewer);
442:   } else if (isascii) {
443:     PetscViewerGetFormat(viewer,&format);
444:     switch (format) {
445:       case PETSC_VIEWER_DEFAULT:
446:       case PETSC_VIEWER_ASCII_INFO:
447:       case PETSC_VIEWER_ASCII_INFO_DETAIL:
448:         SVDValuesView_ASCII(svd,viewer);
449:         break;
450:       case PETSC_VIEWER_ASCII_MATLAB:
451:         SVDValuesView_MATLAB(svd,viewer);
452:         break;
453:       default:
454:         PetscInfo1(svd,"Unsupported viewer format %s\n",PetscViewerFormats[format]);
455:     }
456:   }
457:   return(0);
458: }

460: /*@
461:    SVDValuesViewFromOptions - Processes command line options to determine if/how
462:    the computed singular values are to be viewed.

464:    Collective on SVD

466:    Input Parameters:
467: .  svd - the singular value solver context

469:    Level: developer
470: @*/
471: PetscErrorCode SVDValuesViewFromOptions(SVD svd)
472: {
473:   PetscErrorCode    ierr;
474:   PetscViewer       viewer;
475:   PetscBool         flg;
476:   static PetscBool  incall = PETSC_FALSE;
477:   PetscViewerFormat format;

480:   if (incall) return(0);
481:   incall = PETSC_TRUE;
482:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_view_values",&viewer,&format,&flg);
483:   if (flg) {
484:     PetscViewerPushFormat(viewer,format);
485:     SVDValuesView(svd,viewer);
486:     PetscViewerPopFormat(viewer);
487:     PetscViewerDestroy(&viewer);
488:   }
489:   incall = PETSC_FALSE;
490:   return(0);
491: }

493: /*@C
494:    SVDVectorsView - Outputs computed singular vectors to a viewer.

496:    Collective on SVD

498:    Parameter:
499: +  svd    - the singular value solver context
500: -  viewer - the viewer

502:    Options Database Keys:
503: .  -svd_view_vectors - output singular vectors

505:    Note:
506:    Right and left singular vectors are interleaved, that is, the vectors are
507:    output in the following order: V0, U0, V1, U1, V2, U2, ...

509:    Level: intermediate

511: .seealso: SVDSolve(), SVDValuesView(), SVDErrorView()
512: @*/
513: PetscErrorCode SVDVectorsView(SVD svd,PetscViewer viewer)
514: {
516:   PetscInt       i,k;
517:   Vec            x;
518: #define NMLEN 30
519:   char           vname[NMLEN];
520:   const char     *ename;

524:   if (!viewer) {
525:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer);
526:   }
529:   SVDCheckSolved(svd,1);
530:   if (svd->nconv) {
531:     PetscObjectGetName((PetscObject)svd,&ename);
532:     SVDComputeVectors(svd);
533:     for (i=0;i<svd->nconv;i++) {
534:       k = svd->perm[i];
535:       PetscSNPrintf(vname,NMLEN,"V%d_%s",(int)i,ename);
536:       BVGetColumn(svd->V,k,&x);
537:       PetscObjectSetName((PetscObject)x,vname);
538:       VecView(x,viewer);
539:       BVRestoreColumn(svd->V,k,&x);
540:       PetscSNPrintf(vname,NMLEN,"U%d_%s",(int)i,ename);
541:       BVGetColumn(svd->U,k,&x);
542:       PetscObjectSetName((PetscObject)x,vname);
543:       VecView(x,viewer);
544:       BVRestoreColumn(svd->U,k,&x);
545:     }
546:   }
547:   return(0);
548: }

550: /*@
551:    SVDVectorsViewFromOptions - Processes command line options to determine if/how
552:    the computed singular vectors are to be viewed.

554:    Collective on SVD

556:    Input Parameters:
557: .  svd - the singular value solver context

559:    Level: developer
560: @*/
561: PetscErrorCode SVDVectorsViewFromOptions(SVD svd)
562: {
563:   PetscErrorCode    ierr;
564:   PetscViewer       viewer;
565:   PetscBool         flg = PETSC_FALSE;
566:   static PetscBool  incall = PETSC_FALSE;
567:   PetscViewerFormat format;

570:   if (incall) return(0);
571:   incall = PETSC_TRUE;
572:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_view_vectors",&viewer,&format,&flg);
573:   if (flg) {
574:     PetscViewerPushFormat(viewer,format);
575:     SVDVectorsView(svd,viewer);
576:     PetscViewerPopFormat(viewer);
577:     PetscViewerDestroy(&viewer);
578:   }
579:   incall = PETSC_FALSE;
580:   return(0);
581: }