Actual source code: svdview.c

slepc-main 2024-12-17
Report Typos and Errors
  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 SVD routines related to various viewers
 12: */

 14: #include <slepc/private/svdimpl.h>
 15: #include <petscdraw.h>

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

 20:    Collective

 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: EPSView()
 43: @*/
 44: PetscErrorCode SVDView(SVD svd,PetscViewer viewer)
 45: {
 46:   const char     *type=NULL;
 47:   PetscBool      isascii,isshell,isexternal;

 49:   PetscFunctionBegin;
 51:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer));
 53:   PetscCheckSameComm(svd,1,viewer,2);

 55:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
 56:   if (isascii) {
 57:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)svd,viewer));
 58:     PetscCall(PetscViewerASCIIPushTab(viewer));
 59:     PetscTryTypeMethod(svd,view,viewer);
 60:     PetscCall(PetscViewerASCIIPopTab(viewer));
 61:     if (svd->problem_type) {
 62:       switch (svd->problem_type) {
 63:         case SVD_STANDARD:    type = "(standard) singular value problem"; break;
 64:         case SVD_GENERALIZED: type = "generalized singular value problem"; break;
 65:         case SVD_HYPERBOLIC:  type = "hyperbolic singular value problem"; break;
 66:       }
 67:     } else type = "not yet set";
 68:     PetscCall(PetscViewerASCIIPrintf(viewer,"  problem type: %s\n",type));
 69:     PetscCall(PetscViewerASCIIPrintf(viewer,"  transpose mode: %s\n",svd->impltrans?"implicit":"explicit"));
 70:     if (svd->which == SVD_LARGEST) PetscCall(PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: largest\n"));
 71:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: smallest\n"));
 72:     if (svd->stop==SVD_STOP_THRESHOLD) PetscCall(PetscViewerASCIIPrintf(viewer,"  computing singular values %s the threshold: %g%s\n",svd->which==SVD_LARGEST?"above":"below",(double)svd->thres,svd->threlative?" (relative)":""));
 73:     if (svd->nsv) PetscCall(PetscViewerASCIIPrintf(viewer,"  number of singular values (nsv): %" PetscInt_FMT "\n",svd->nsv));
 74:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %" PetscInt_FMT "\n",svd->ncv));
 75:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum dimension of projected problem (mpd): %" PetscInt_FMT "\n",svd->mpd));
 76:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %" PetscInt_FMT "\n",svd->max_it));
 77:     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",(double)svd->tol));
 78:     PetscCall(PetscViewerASCIIPrintf(viewer,"  convergence test: "));
 79:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
 80:     switch (svd->conv) {
 81:     case SVD_CONV_ABS:
 82:       PetscCall(PetscViewerASCIIPrintf(viewer,"absolute\n"));break;
 83:     case SVD_CONV_REL:
 84:       PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the singular value\n"));break;
 85:     case SVD_CONV_NORM:
 86:       PetscCall(PetscViewerASCIIPrintf(viewer,"relative to the matrix norms\n"));
 87:       PetscCall(PetscViewerASCIIPrintf(viewer,"  computed matrix norms: norm(A)=%g",(double)svd->nrma));
 88:       if (svd->isgeneralized) PetscCall(PetscViewerASCIIPrintf(viewer,", norm(B)=%g",(double)svd->nrmb));
 89:       PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
 90:       break;
 91:     case SVD_CONV_MAXIT:
 92:       PetscCall(PetscViewerASCIIPrintf(viewer,"maximum number of iterations\n"));break;
 93:     case SVD_CONV_USER:
 94:       PetscCall(PetscViewerASCIIPrintf(viewer,"user-defined\n"));break;
 95:     }
 96:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
 97:     if (svd->nini) PetscCall(PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial space: %" PetscInt_FMT "\n",PetscAbs(svd->nini)));
 98:     if (svd->ninil) PetscCall(PetscViewerASCIIPrintf(viewer,"  dimension of user-provided initial left space: %" PetscInt_FMT "\n",PetscAbs(svd->ninil)));
 99:   } else PetscTryTypeMethod(svd,view,viewer);
100:   PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&isshell,SVDCROSS,SVDCYCLIC,""));
101:   PetscCall(PetscObjectTypeCompareAny((PetscObject)svd,&isexternal,SVDSCALAPACK,SVDKSVD,SVDELEMENTAL,SVDPRIMME,""));
102:   if (!isshell && !isexternal) {
103:     PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO));
104:     if (!svd->V) PetscCall(SVDGetBV(svd,&svd->V,NULL));
105:     PetscCall(BVView(svd->V,viewer));
106:     if (!svd->ds) PetscCall(SVDGetDS(svd,&svd->ds));
107:     PetscCall(DSView(svd->ds,viewer));
108:     PetscCall(PetscViewerPopFormat(viewer));
109:   }
110:   PetscFunctionReturn(PETSC_SUCCESS);
111: }

113: /*@
114:    SVDViewFromOptions - View from options

116:    Collective

118:    Input Parameters:
119: +  svd  - the singular value solver context
120: .  obj  - optional object
121: -  name - command line option

123:    Level: intermediate

125: .seealso: SVDView(), SVDCreate()
126: @*/
127: PetscErrorCode SVDViewFromOptions(SVD svd,PetscObject obj,const char name[])
128: {
129:   PetscFunctionBegin;
131:   PetscCall(PetscObjectViewFromOptions((PetscObject)svd,obj,name));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: /*@
136:    SVDConvergedReasonView - Displays the reason an SVD solve converged or diverged.

138:    Collective

140:    Input Parameters:
141: +  svd - the singular value solver context
142: -  viewer - the viewer to display the reason

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

147:    Note:
148:    To change the format of the output call PetscViewerPushFormat(viewer,format) before
149:    this call. Use PETSC_VIEWER_DEFAULT for the default, use PETSC_VIEWER_FAILED to only
150:    display a reason if it fails. The latter can be set in the command line with
151:    -svd_converged_reason ::failed

153:    Level: intermediate

155: .seealso: SVDSetTolerances(), SVDGetIterationNumber(), SVDConvergedReasonViewFromOptions()
156: @*/
157: PetscErrorCode SVDConvergedReasonView(SVD svd,PetscViewer viewer)
158: {
159:   PetscBool         isAscii;
160:   PetscViewerFormat format;

162:   PetscFunctionBegin;
163:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)svd));
164:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii));
165:   if (isAscii) {
166:     PetscCall(PetscViewerGetFormat(viewer,&format));
167:     PetscCall(PetscViewerASCIIAddTab(viewer,((PetscObject)svd)->tablevel));
168:     if (svd->reason > 0 && format != PETSC_VIEWER_FAILED) PetscCall(PetscViewerASCIIPrintf(viewer,"%s SVD solve converged (%" PetscInt_FMT " singular triplet%s) due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)svd)->prefix?((PetscObject)svd)->prefix:"",svd->nconv,(svd->nconv>1)?"s":"",SVDConvergedReasons[svd->reason],svd->its));
169:     else if (svd->reason <= 0) PetscCall(PetscViewerASCIIPrintf(viewer,"%s SVD solve did not converge due to %s; iterations %" PetscInt_FMT "\n",((PetscObject)svd)->prefix?((PetscObject)svd)->prefix:"",SVDConvergedReasons[svd->reason],svd->its));
170:     PetscCall(PetscViewerASCIISubtractTab(viewer,((PetscObject)svd)->tablevel));
171:   }
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: /*@
176:    SVDConvergedReasonViewFromOptions - Processes command line options to determine if/how
177:    the SVD converged reason is to be viewed.

179:    Collective

181:    Input Parameter:
182: .  svd - the singular value solver context

184:    Level: developer

186: .seealso: SVDConvergedReasonView()
187: @*/
188: PetscErrorCode SVDConvergedReasonViewFromOptions(SVD svd)
189: {
190:   PetscViewer       viewer;
191:   PetscBool         flg;
192:   static PetscBool  incall = PETSC_FALSE;
193:   PetscViewerFormat format;

195:   PetscFunctionBegin;
196:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
197:   incall = PETSC_TRUE;
198:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_converged_reason",&viewer,&format,&flg));
199:   if (flg) {
200:     PetscCall(PetscViewerPushFormat(viewer,format));
201:     PetscCall(SVDConvergedReasonView(svd,viewer));
202:     PetscCall(PetscViewerPopFormat(viewer));
203:     PetscCall(PetscViewerDestroy(&viewer));
204:   }
205:   incall = PETSC_FALSE;
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: static PetscErrorCode SVDErrorView_ASCII(SVD svd,SVDErrorType etype,PetscViewer viewer)
210: {
211:   PetscReal      error,sigma;
212:   PetscInt       i,j,numsv;

214:   PetscFunctionBegin;
215:   numsv = svd->stop==SVD_STOP_THRESHOLD? svd->nconv: svd->nsv;
216:   if (svd->nconv<numsv) {
217:     PetscCall(PetscViewerASCIIPrintf(viewer," Problem: less than %" PetscInt_FMT " singular values converged\n\n",numsv));
218:     PetscFunctionReturn(PETSC_SUCCESS);
219:   }
220:   for (i=0;i<numsv;i++) {
221:     PetscCall(SVDComputeError(svd,i,etype,&error));
222:     if (error>=5.0*svd->tol) {
223:       PetscCall(PetscViewerASCIIPrintf(viewer," Problem: some of the first %" PetscInt_FMT " relative errors are higher than the tolerance\n\n",numsv));
224:       PetscFunctionReturn(PETSC_SUCCESS);
225:     }
226:   }
227:   PetscCall(PetscViewerASCIIPrintf(viewer," All %s%ssingular values computed up to the required tolerance:",svd->stop==SVD_STOP_THRESHOLD?"":"requested ",svd->isgeneralized?"generalized ":""));
228:   for (i=0;i<=(numsv-1)/8;i++) {
229:     PetscCall(PetscViewerASCIIPrintf(viewer,"\n     "));
230:     for (j=0;j<PetscMin(8,numsv-8*i);j++) {
231:       PetscCall(SVDGetSingularTriplet(svd,8*i+j,&sigma,NULL,NULL));
232:       PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f",(double)sigma));
233:       if (8*i+j+1<numsv) PetscCall(PetscViewerASCIIPrintf(viewer,", "));
234:     }
235:   }
236:   PetscCall(PetscViewerASCIIPrintf(viewer,"\n\n"));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: static PetscErrorCode SVDErrorView_DETAIL(SVD svd,SVDErrorType etype,PetscViewer viewer)
241: {
242:   PetscReal      error,sigma;
243:   PetscInt       i;
244:   char           ex[30],sep[]=" ---------------------- --------------------\n";

246:   PetscFunctionBegin;
247:   if (!svd->nconv) PetscFunctionReturn(PETSC_SUCCESS);
248:   switch (etype) {
249:     case SVD_ERROR_ABSOLUTE:
250:       PetscCall(PetscSNPrintf(ex,sizeof(ex)," absolute error"));
251:       break;
252:     case SVD_ERROR_RELATIVE:
253:       PetscCall(PetscSNPrintf(ex,sizeof(ex)," relative error"));
254:       break;
255:     case SVD_ERROR_NORM:
256:       if (svd->isgeneralized) PetscCall(PetscSNPrintf(ex,sizeof(ex)," ||r||/||[A;B]||"));
257:       else PetscCall(PetscSNPrintf(ex,sizeof(ex),"  ||r||/||A||"));
258:       break;
259:   }
260:   PetscCall(PetscViewerASCIIPrintf(viewer,"%s          sigma           %s\n%s",sep,ex,sep));
261:   for (i=0;i<svd->nconv;i++) {
262:     PetscCall(SVDGetSingularTriplet(svd,i,&sigma,NULL,NULL));
263:     PetscCall(SVDComputeError(svd,i,etype,&error));
264:     PetscCall(PetscViewerASCIIPrintf(viewer,"       % 6f          %12g\n",(double)sigma,(double)error));
265:   }
266:   PetscCall(PetscViewerASCIIPrintf(viewer,"%s",sep));
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

270: static PetscErrorCode SVDErrorView_MATLAB(SVD svd,SVDErrorType etype,PetscViewer viewer)
271: {
272:   PetscReal      error;
273:   PetscInt       i;
274:   const char     *name;

276:   PetscFunctionBegin;
277:   PetscCall(PetscObjectGetName((PetscObject)svd,&name));
278:   PetscCall(PetscViewerASCIIPrintf(viewer,"Error_%s = [\n",name));
279:   for (i=0;i<svd->nconv;i++) {
280:     PetscCall(SVDComputeError(svd,i,etype,&error));
281:     PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)error));
282:   }
283:   PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: /*@
288:    SVDErrorView - Displays the errors associated with the computed solution
289:    (as well as the singular values).

291:    Collective

293:    Input Parameters:
294: +  svd    - the singular value solver context
295: .  etype  - error type
296: -  viewer - optional visualization context

298:    Options Database Keys:
299: +  -svd_error_absolute - print absolute errors of each singular triplet
300: .  -svd_error_relative - print relative errors of each singular triplet
301: -  -svd_error_norm     - print errors relative to the matrix norms of each singular triplet

303:    Notes:
304:    By default, this function checks the error of all singular triplets and prints
305:    the singular values if all of them are below the requested tolerance.
306:    If the viewer has format=PETSC_VIEWER_ASCII_INFO_DETAIL then a table with
307:    singular values and corresponding errors is printed.

309:    Level: intermediate

311: .seealso: SVDSolve(), SVDValuesView(), SVDVectorsView()
312: @*/
313: PetscErrorCode SVDErrorView(SVD svd,SVDErrorType etype,PetscViewer viewer)
314: {
315:   PetscBool         isascii;
316:   PetscViewerFormat format;

318:   PetscFunctionBegin;
320:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer));
322:   PetscCheckSameComm(svd,1,viewer,3);
323:   SVDCheckSolved(svd,1);
324:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
325:   if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);

327:   PetscCall(PetscViewerGetFormat(viewer,&format));
328:   switch (format) {
329:     case PETSC_VIEWER_DEFAULT:
330:     case PETSC_VIEWER_ASCII_INFO:
331:       PetscCall(SVDErrorView_ASCII(svd,etype,viewer));
332:       break;
333:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
334:       PetscCall(SVDErrorView_DETAIL(svd,etype,viewer));
335:       break;
336:     case PETSC_VIEWER_ASCII_MATLAB:
337:       PetscCall(SVDErrorView_MATLAB(svd,etype,viewer));
338:       break;
339:     default:
340:       PetscCall(PetscInfo(svd,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
341:   }
342:   PetscFunctionReturn(PETSC_SUCCESS);
343: }

345: /*@
346:    SVDErrorViewFromOptions - Processes command line options to determine if/how
347:    the errors of the computed solution are to be viewed.

349:    Collective

351:    Input Parameter:
352: .  svd - the singular value solver context

354:    Level: developer

356: .seealso: SVDErrorView()
357: @*/
358: PetscErrorCode SVDErrorViewFromOptions(SVD svd)
359: {
360:   PetscViewer       viewer;
361:   PetscBool         flg;
362:   static PetscBool  incall = PETSC_FALSE;
363:   PetscViewerFormat format;

365:   PetscFunctionBegin;
366:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
367:   incall = PETSC_TRUE;
368:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_error_absolute",&viewer,&format,&flg));
369:   if (flg) {
370:     PetscCall(PetscViewerPushFormat(viewer,format));
371:     PetscCall(SVDErrorView(svd,SVD_ERROR_ABSOLUTE,viewer));
372:     PetscCall(PetscViewerPopFormat(viewer));
373:     PetscCall(PetscViewerDestroy(&viewer));
374:   }
375:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_error_relative",&viewer,&format,&flg));
376:   if (flg) {
377:     PetscCall(PetscViewerPushFormat(viewer,format));
378:     PetscCall(SVDErrorView(svd,SVD_ERROR_RELATIVE,viewer));
379:     PetscCall(PetscViewerPopFormat(viewer));
380:     PetscCall(PetscViewerDestroy(&viewer));
381:   }
382:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_error_norm",&viewer,&format,&flg));
383:   if (flg) {
384:     PetscCall(PetscViewerPushFormat(viewer,format));
385:     PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,viewer));
386:     PetscCall(PetscViewerPopFormat(viewer));
387:     PetscCall(PetscViewerDestroy(&viewer));
388:   }
389:   incall = PETSC_FALSE;
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }

393: static PetscErrorCode SVDValuesView_DRAW(SVD svd,PetscViewer viewer)
394: {
395:   PetscDraw      draw;
396:   PetscDrawSP    drawsp;
397:   PetscReal      re,im=0.0;
398:   PetscInt       i;

400:   PetscFunctionBegin;
401:   if (!svd->nconv) PetscFunctionReturn(PETSC_SUCCESS);
402:   PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
403:   PetscCall(PetscDrawSetTitle(draw,"Computed singular values"));
404:   PetscCall(PetscDrawSPCreate(draw,1,&drawsp));
405:   for (i=0;i<svd->nconv;i++) {
406:     re = svd->sigma[svd->perm[i]];
407:     PetscCall(PetscDrawSPAddPoint(drawsp,&re,&im));
408:   }
409:   PetscCall(PetscDrawSPDraw(drawsp,PETSC_TRUE));
410:   PetscCall(PetscDrawSPSave(drawsp));
411:   PetscCall(PetscDrawSPDestroy(&drawsp));
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: static PetscErrorCode SVDValuesView_BINARY(SVD svd,PetscViewer viewer)
416: {
417:   PetscInt       i,k;
418:   PetscReal      *sv;

420:   PetscFunctionBegin;
421:   PetscCall(PetscMalloc1(svd->nconv,&sv));
422:   for (i=0;i<svd->nconv;i++) {
423:     k = svd->perm[i];
424:     sv[i] = svd->sigma[k];
425:   }
426:   PetscCall(PetscViewerBinaryWrite(viewer,sv,svd->nconv,PETSC_REAL));
427:   PetscCall(PetscFree(sv));
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: #if defined(PETSC_HAVE_HDF5)
432: static PetscErrorCode SVDValuesView_HDF5(SVD svd,PetscViewer viewer)
433: {
434:   PetscInt       i,k,n,N;
435:   PetscMPIInt    rank;
436:   Vec            v;
437:   char           vname[30];
438:   const char     *ename;

440:   PetscFunctionBegin;
441:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)svd),&rank));
442:   N = svd->nconv;
443:   n = rank? 0: N;
444:   /* create a vector containing the singular values */
445:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)svd),n,N,&v));
446:   PetscCall(PetscObjectGetName((PetscObject)svd,&ename));
447:   PetscCall(PetscSNPrintf(vname,sizeof(vname),"sigma_%s",ename));
448:   PetscCall(PetscObjectSetName((PetscObject)v,vname));
449:   if (!rank) {
450:     for (i=0;i<svd->nconv;i++) {
451:       k = svd->perm[i];
452:       PetscCall(VecSetValue(v,i,svd->sigma[k],INSERT_VALUES));
453:     }
454:   }
455:   PetscCall(VecAssemblyBegin(v));
456:   PetscCall(VecAssemblyEnd(v));
457:   PetscCall(VecView(v,viewer));
458:   PetscCall(VecDestroy(&v));
459:   PetscFunctionReturn(PETSC_SUCCESS);
460: }
461: #endif

463: static PetscErrorCode SVDValuesView_ASCII(SVD svd,PetscViewer viewer)
464: {
465:   PetscInt       i;

467:   PetscFunctionBegin;
468:   PetscCall(PetscViewerASCIIPrintf(viewer,"Singular values = \n"));
469:   for (i=0;i<svd->nconv;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"   %.5f\n",(double)svd->sigma[svd->perm[i]]));
470:   PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
471:   PetscFunctionReturn(PETSC_SUCCESS);
472: }

474: static PetscErrorCode SVDValuesView_MATLAB(SVD svd,PetscViewer viewer)
475: {
476:   PetscInt       i;
477:   const char     *name;

479:   PetscFunctionBegin;
480:   PetscCall(PetscObjectGetName((PetscObject)svd,&name));
481:   PetscCall(PetscViewerASCIIPrintf(viewer,"Sigma_%s = [\n",name));
482:   for (i=0;i<svd->nconv;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)svd->sigma[svd->perm[i]]));
483:   PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
484:   PetscFunctionReturn(PETSC_SUCCESS);
485: }

487: /*@
488:    SVDValuesView - Displays the computed singular values in a viewer.

490:    Collective

492:    Input Parameters:
493: +  svd    - the singular value solver context
494: -  viewer - the viewer

496:    Options Database Key:
497: .  -svd_view_values - print computed singular values

499:    Level: intermediate

501: .seealso: SVDSolve(), SVDVectorsView(), SVDErrorView()
502: @*/
503: PetscErrorCode SVDValuesView(SVD svd,PetscViewer viewer)
504: {
505:   PetscBool         isascii,isdraw,isbinary;
506:   PetscViewerFormat format;
507: #if defined(PETSC_HAVE_HDF5)
508:   PetscBool         ishdf5;
509: #endif

511:   PetscFunctionBegin;
513:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer));
515:   PetscCheckSameComm(svd,1,viewer,2);
516:   SVDCheckSolved(svd,1);
517:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
518:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary));
519: #if defined(PETSC_HAVE_HDF5)
520:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5));
521: #endif
522:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
523:   if (isdraw) PetscCall(SVDValuesView_DRAW(svd,viewer));
524:   else if (isbinary) PetscCall(SVDValuesView_BINARY(svd,viewer));
525: #if defined(PETSC_HAVE_HDF5)
526:   else if (ishdf5) PetscCall(SVDValuesView_HDF5(svd,viewer));
527: #endif
528:   else if (isascii) {
529:     PetscCall(PetscViewerGetFormat(viewer,&format));
530:     switch (format) {
531:       case PETSC_VIEWER_DEFAULT:
532:       case PETSC_VIEWER_ASCII_INFO:
533:       case PETSC_VIEWER_ASCII_INFO_DETAIL:
534:         PetscCall(SVDValuesView_ASCII(svd,viewer));
535:         break;
536:       case PETSC_VIEWER_ASCII_MATLAB:
537:         PetscCall(SVDValuesView_MATLAB(svd,viewer));
538:         break;
539:       default:
540:         PetscCall(PetscInfo(svd,"Unsupported viewer format %s\n",PetscViewerFormats[format]));
541:     }
542:   }
543:   PetscFunctionReturn(PETSC_SUCCESS);
544: }

546: /*@
547:    SVDValuesViewFromOptions - Processes command line options to determine if/how
548:    the computed singular values are to be viewed.

550:    Collective

552:    Input Parameter:
553: .  svd - the singular value solver context

555:    Level: developer

557: .seealso: SVDValuesView()
558: @*/
559: PetscErrorCode SVDValuesViewFromOptions(SVD svd)
560: {
561:   PetscViewer       viewer;
562:   PetscBool         flg;
563:   static PetscBool  incall = PETSC_FALSE;
564:   PetscViewerFormat format;

566:   PetscFunctionBegin;
567:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
568:   incall = PETSC_TRUE;
569:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_view_values",&viewer,&format,&flg));
570:   if (flg) {
571:     PetscCall(PetscViewerPushFormat(viewer,format));
572:     PetscCall(SVDValuesView(svd,viewer));
573:     PetscCall(PetscViewerPopFormat(viewer));
574:     PetscCall(PetscViewerDestroy(&viewer));
575:   }
576:   incall = PETSC_FALSE;
577:   PetscFunctionReturn(PETSC_SUCCESS);
578: }

580: /*@
581:    SVDVectorsView - Outputs computed singular vectors to a viewer.

583:    Collective

585:    Input Parameters:
586: +  svd    - the singular value solver context
587: -  viewer - the viewer

589:    Options Database Key:
590: .  -svd_view_vectors - output singular vectors

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

596:    Level: intermediate

598: .seealso: SVDSolve(), SVDValuesView(), SVDErrorView()
599: @*/
600: PetscErrorCode SVDVectorsView(SVD svd,PetscViewer viewer)
601: {
602:   PetscInt       i,k;
603:   Vec            x;
604:   char           vname[30];
605:   const char     *ename;

607:   PetscFunctionBegin;
609:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)svd),&viewer));
611:   PetscCheckSameComm(svd,1,viewer,2);
612:   SVDCheckSolved(svd,1);
613:   if (svd->nconv) {
614:     PetscCall(PetscObjectGetName((PetscObject)svd,&ename));
615:     PetscCall(SVDComputeVectors(svd));
616:     for (i=0;i<svd->nconv;i++) {
617:       k = svd->perm[i];
618:       PetscCall(PetscSNPrintf(vname,sizeof(vname),"V%" PetscInt_FMT "_%s",i,ename));
619:       PetscCall(BVGetColumn(svd->V,k,&x));
620:       PetscCall(PetscObjectSetName((PetscObject)x,vname));
621:       PetscCall(VecView(x,viewer));
622:       PetscCall(BVRestoreColumn(svd->V,k,&x));
623:       PetscCall(PetscSNPrintf(vname,sizeof(vname),"U%" PetscInt_FMT "_%s",i,ename));
624:       PetscCall(BVGetColumn(svd->U,k,&x));
625:       PetscCall(PetscObjectSetName((PetscObject)x,vname));
626:       PetscCall(VecView(x,viewer));
627:       PetscCall(BVRestoreColumn(svd->U,k,&x));
628:     }
629:   }
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: /*@
634:    SVDVectorsViewFromOptions - Processes command line options to determine if/how
635:    the computed singular vectors are to be viewed.

637:    Collective

639:    Input Parameter:
640: .  svd - the singular value solver context

642:    Level: developer

644: .seealso: SVDVectorsView()
645: @*/
646: PetscErrorCode SVDVectorsViewFromOptions(SVD svd)
647: {
648:   PetscViewer       viewer;
649:   PetscBool         flg = PETSC_FALSE;
650:   static PetscBool  incall = PETSC_FALSE;
651:   PetscViewerFormat format;

653:   PetscFunctionBegin;
654:   if (incall) PetscFunctionReturn(PETSC_SUCCESS);
655:   incall = PETSC_TRUE;
656:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)svd),((PetscObject)svd)->options,((PetscObject)svd)->prefix,"-svd_view_vectors",&viewer,&format,&flg));
657:   if (flg) {
658:     PetscCall(PetscViewerPushFormat(viewer,format));
659:     PetscCall(SVDVectorsView(svd,viewer));
660:     PetscCall(PetscViewerPopFormat(viewer));
661:     PetscCall(PetscViewerDestroy(&viewer));
662:   }
663:   incall = PETSC_FALSE;
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }