Actual source code: hdf5v.c

petsc-3.11.1 2019-04-17
Report Typos and Errors
  1: #include <petsc/private/viewerimpl.h>    /*I   "petscsys.h"   I*/
  2: #include <petscviewerhdf5.h>    /*I   "petscviewerhdf5.h"   I*/
  3: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE < 10800)
  4: #error "PETSc needs HDF5 version >= 1.8.0"
  5: #endif

  7: static PetscErrorCode PetscViewerHDF5Traverse_Internal(PetscViewer, const char[], PetscBool, PetscBool*, H5O_type_t*);
  8: static PetscErrorCode PetscViewerHDF5HasAttribute_Internal(PetscViewer, const char[], const char[], PetscBool*);

 10: typedef struct GroupList {
 11:   const char       *name;
 12:   struct GroupList *next;
 13: } GroupList;

 15: typedef struct {
 16:   char          *filename;
 17:   PetscFileMode btype;
 18:   hid_t         file_id;
 19:   PetscInt      timestep;
 20:   GroupList     *groups;
 21:   PetscBool     basedimension2;  /* save vectors and DMDA vectors with a dimension of at least 2 even if the bs/dof is 1 */
 22:   PetscBool     spoutput;  /* write data in single precision even if PETSc is compiled with double precision PetscReal */
 23:   char          *mataij_iname;
 24:   char          *mataij_jname;
 25:   char          *mataij_aname;
 26:   char          *mataij_cname;
 27:   PetscBool     mataij_names_set;
 28: } PetscViewer_HDF5;

 30: struct _n_HDF5ReadCtx {
 31:   hid_t file, group, dataset, dataspace, plist;
 32:   PetscInt timestep;
 33:   PetscBool complexVal, dim2, horizontal;
 34: };
 35: typedef struct _n_HDF5ReadCtx* HDF5ReadCtx;

 37: static PetscErrorCode PetscViewerHDF5GetAbsolutePath_Internal(PetscViewer viewer, const char objname[], char **fullpath)
 38: {
 39:   const char *group;
 40:   char buf[PETSC_MAX_PATH_LEN]="";

 44:   PetscViewerHDF5GetGroup(viewer, &group);
 45:   PetscStrcat(buf, group);
 46:   PetscStrcat(buf, "/");
 47:   PetscStrcat(buf, objname);
 48:   PetscStrallocpy(buf, fullpath);
 49:   return(0);
 50: }

 52: static PetscErrorCode PetscViewerHDF5CheckNamedObject_Internal(PetscViewer viewer, PetscObject obj)
 53: {
 54:   PetscBool has;
 55:   const char *group;

 59:   if (!obj->name) SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Object must be named");
 60:   PetscViewerHDF5HasObject(viewer, obj, &has);
 61:   if (!has) {
 62:     PetscViewerHDF5GetGroup(viewer, &group);
 63:     SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) %s not stored in group %s", obj->name, group);
 64:   }
 65:   return(0);
 66: }

 68: static PetscErrorCode PetscViewerSetFromOptions_HDF5(PetscOptionItems *PetscOptionsObject,PetscViewer v)
 69: {
 70:   PetscErrorCode   ierr;
 71:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*)v->data;

 74:   PetscOptionsHead(PetscOptionsObject,"HDF5 PetscViewer Options");
 75:   PetscOptionsBool("-viewer_hdf5_base_dimension2","1d Vectors get 2 dimensions in HDF5","PetscViewerHDF5SetBaseDimension2",hdf5->basedimension2,&hdf5->basedimension2,NULL);
 76:   PetscOptionsBool("-viewer_hdf5_sp_output","Force data to be written in single precision","PetscViewerHDF5SetSPOutput",hdf5->spoutput,&hdf5->spoutput,NULL);
 77:   PetscOptionsTail();
 78:   return(0);
 79: }

 81: static PetscErrorCode PetscViewerFileClose_HDF5(PetscViewer viewer)
 82: {
 83:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*)viewer->data;
 84:   PetscErrorCode   ierr;

 87:   PetscFree(hdf5->filename);
 88:   if (hdf5->file_id) PetscStackCallHDF5(H5Fclose,(hdf5->file_id));
 89:   return(0);
 90: }

 92: static PetscErrorCode PetscViewerDestroy_HDF5(PetscViewer viewer)
 93: {
 94:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
 95:   PetscErrorCode   ierr;

 98:   PetscViewerFileClose_HDF5(viewer);
 99:   while (hdf5->groups) {
100:     GroupList *tmp = hdf5->groups->next;

102:     PetscFree(hdf5->groups->name);
103:     PetscFree(hdf5->groups);
104:     hdf5->groups = tmp;
105:   }
106:   PetscFree(hdf5->mataij_iname);
107:   PetscFree(hdf5->mataij_jname);
108:   PetscFree(hdf5->mataij_aname);
109:   PetscFree(hdf5->mataij_cname);
110:   PetscFree(hdf5);
111:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);
112:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);
113:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);
114:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerHDF5SetBaseDimension2_C",NULL);
115:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerHDF5SetSPOutput_C",NULL);
116:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerHDF5SetAIJNames_C",NULL);
117:   PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerHDF5GetAIJNames_C",NULL);
118:   return(0);
119: }

121: static PetscErrorCode  PetscViewerFileSetMode_HDF5(PetscViewer viewer, PetscFileMode type)
122: {
123:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

126:   hdf5->btype = type;
127:   return(0);
128: }

130: static PetscErrorCode  PetscViewerFileGetMode_HDF5(PetscViewer viewer, PetscFileMode *type)
131: {
132:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

135:   *type = hdf5->btype;
136:   return(0);
137: }

139: static PetscErrorCode  PetscViewerHDF5SetBaseDimension2_HDF5(PetscViewer viewer, PetscBool flg)
140: {
141:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

144:   hdf5->basedimension2 = flg;
145:   return(0);
146: }

148: /*@
149:      PetscViewerHDF5SetBaseDimension2 - Vectors of 1 dimension (i.e. bs/dof is 1) will be saved in the HDF5 file with a
150:        dimension of 2.

152:     Logically Collective on PetscViewer

154:   Input Parameters:
155: +  viewer - the PetscViewer; if it is not hdf5 then this command is ignored
156: -  flg - if PETSC_TRUE the vector will always have at least a dimension of 2 even if that first dimension is of size 1

158:   Options Database:
159: .  -viewer_hdf5_base_dimension2 - turns on (true) or off (false) using a dimension of 2 in the HDF5 file even if the bs/dof of the vector is 1


162:   Notes:
163:     Setting this option allegedly makes code that reads the HDF5 in easier since they do not have a "special case" of a bs/dof
164:          of one when the dimension is lower. Others think the option is crazy.

166:   Level: intermediate

168: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()

170: @*/
171: PetscErrorCode PetscViewerHDF5SetBaseDimension2(PetscViewer viewer,PetscBool flg)
172: {

177:   PetscTryMethod(viewer,"PetscViewerHDF5SetBaseDimension2_C",(PetscViewer,PetscBool),(viewer,flg));
178:   return(0);
179: }

181: /*@
182:      PetscViewerHDF5GetBaseDimension2 - Vectors of 1 dimension (i.e. bs/dof is 1) will be saved in the HDF5 file with a
183:        dimension of 2.

185:     Logically Collective on PetscViewer

187:   Input Parameter:
188: .  viewer - the PetscViewer, must be of type HDF5

190:   Output Parameter:
191: .  flg - if PETSC_TRUE the vector will always have at least a dimension of 2 even if that first dimension is of size 1

193:   Notes:
194:     Setting this option allegedly makes code that reads the HDF5 in easier since they do not have a "special case" of a bs/dof
195:          of one when the dimension is lower. Others think the option is crazy.

197:   Level: intermediate

199: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()

201: @*/
202: PetscErrorCode PetscViewerHDF5GetBaseDimension2(PetscViewer viewer,PetscBool *flg)
203: {
204:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

208:   *flg = hdf5->basedimension2;
209:   return(0);
210: }

212: static PetscErrorCode  PetscViewerHDF5SetSPOutput_HDF5(PetscViewer viewer, PetscBool flg)
213: {
214:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

217:   hdf5->spoutput = flg;
218:   return(0);
219: }

221: /*@
222:      PetscViewerHDF5SetSPOutput - Data is written to disk in single precision even if PETSc is
223:        compiled with double precision PetscReal.

225:     Logically Collective on PetscViewer

227:   Input Parameters:
228: +  viewer - the PetscViewer; if it is not hdf5 then this command is ignored
229: -  flg - if PETSC_TRUE the data will be written to disk with single precision

231:   Options Database:
232: .  -viewer_hdf5_sp_output - turns on (true) or off (false) output in single precision


235:   Notes:
236:     Setting this option does not make any difference if PETSc is compiled with single precision
237:          in the first place. It does not affect reading datasets (HDF5 handle this internally).

239:   Level: intermediate

241: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen(),
242:           PetscReal

244: @*/
245: PetscErrorCode PetscViewerHDF5SetSPOutput(PetscViewer viewer,PetscBool flg)
246: {

251:   PetscTryMethod(viewer,"PetscViewerHDF5SetSPOutput_C",(PetscViewer,PetscBool),(viewer,flg));
252:   return(0);
253: }

255: /*@
256:      PetscViewerHDF5GetSPOutput - Data is written to disk in single precision even if PETSc is
257:        compiled with double precision PetscReal.

259:     Logically Collective on PetscViewer

261:   Input Parameter:
262: .  viewer - the PetscViewer, must be of type HDF5

264:   Output Parameter:
265: .  flg - if PETSC_TRUE the data will be written to disk with single precision

267:   Notes:
268:     Setting this option does not make any difference if PETSc is compiled with single precision
269:          in the first place. It does not affect reading datasets (HDF5 handle this internally).

271:   Level: intermediate

273: .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen(),
274:           PetscReal

276: @*/
277: PetscErrorCode PetscViewerHDF5GetSPOutput(PetscViewer viewer,PetscBool *flg)
278: {
279:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

283:   *flg = hdf5->spoutput;
284:   return(0);
285: }

287: static PetscErrorCode  PetscViewerFileSetName_HDF5(PetscViewer viewer, const char name[])
288: {
289:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
290: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
291:   MPI_Info          info = MPI_INFO_NULL;
292: #endif
293:   hid_t             plist_id;
294:   PetscErrorCode    ierr;

297:   if (hdf5->file_id) PetscStackCallHDF5(H5Fclose,(hdf5->file_id));
298:   if (hdf5->filename) {PetscFree(hdf5->filename);}
299:   PetscStrallocpy(name, &hdf5->filename);
300:   /* Set up file access property list with parallel I/O access */
301:   PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_FILE_ACCESS));
302: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
303:   PetscStackCallHDF5(H5Pset_fapl_mpio,(plist_id, PetscObjectComm((PetscObject)viewer), info));
304: #endif
305:   /* Create or open the file collectively */
306:   switch (hdf5->btype) {
307:   case FILE_MODE_READ:
308:     PetscStackCallHDF5Return(hdf5->file_id,H5Fopen,(name, H5F_ACC_RDONLY, plist_id));
309:     break;
310:   case FILE_MODE_APPEND:
311:     PetscStackCallHDF5Return(hdf5->file_id,H5Fopen,(name, H5F_ACC_RDWR, plist_id));
312:     break;
313:   case FILE_MODE_WRITE:
314:     PetscStackCallHDF5Return(hdf5->file_id,H5Fcreate,(name, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id));
315:     break;
316:   default:
317:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
318:   }
319:   if (hdf5->file_id < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB, "H5Fcreate failed for %s", name);
320:   PetscStackCallHDF5(H5Pclose,(plist_id));
321:   return(0);
322: }

324: static PetscErrorCode PetscViewerFileGetName_HDF5(PetscViewer viewer,const char **name)
325: {
326:   PetscViewer_HDF5 *vhdf5 = (PetscViewer_HDF5*)viewer->data;

329:   *name = vhdf5->filename;
330:   return(0);
331: }

333: static PetscErrorCode  PetscViewerHDF5SetAIJNames_HDF5(PetscViewer viewer, const char iname[], const char jname[], const char aname[], const char cname[])
334: {
335:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

339:   PetscFree(hdf5->mataij_iname);
340:   PetscFree(hdf5->mataij_jname);
341:   PetscFree(hdf5->mataij_aname);
342:   PetscFree(hdf5->mataij_cname);
343:   PetscStrallocpy(iname,&hdf5->mataij_iname);
344:   PetscStrallocpy(jname,&hdf5->mataij_jname);
345:   PetscStrallocpy(aname,&hdf5->mataij_aname);
346:   PetscStrallocpy(cname,&hdf5->mataij_cname);
347:   hdf5->mataij_names_set = PETSC_TRUE;
348:   return(0);
349: }

351: /*@C
352:   PetscViewerHDF5SetAIJNames - Set the names of the datasets representing the three AIJ (CRS) arrays and the name of the attribute storing the number of columns within the HDF5 file.

354:   Collective on PetscViewer

356:   Input Parameters:
357: +  viewer - the PetscViewer; either ASCII or binary
358: .  iname - name of dataset i representing row pointers; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix
359: .  jname - name of dataset j representing column indices
360: .  aname - name of dataset a representing matrix values
361: -  cname - name of attribute stoting column count

363:   Level: advanced

365:   Notes:
366:   Current defaults are (iname, jname, aname, cname) = ("jc", "ir", "data", "MATLAB_sparse") so that MAT files can be readily loaded.

368: .seealso: MatLoad(), PetscViewerCreate(), PetscViewerSetType(), PETSCVIEWERHDF5, PetscViewerHDF5GetAIJNames()
369: @*/
370: /* TODO Once corresponding MatView is implemented, add this:
371:   Current defaults are (iname, jname, aname, cname) = ("i", "j", "a", "ncols").
372:   For PetscViewerFormat PETSC_VIEWER_HDF5_MAT they are ("jc", "ir", "data", "MATLAB_sparse") so that MAT files can be loaded.
373: */
374: PetscErrorCode  PetscViewerHDF5SetAIJNames(PetscViewer viewer, const char iname[], const char jname[], const char aname[], const char cname[])
375: {

384:   PetscTryMethod(viewer,"PetscViewerHDF5SetAIJNames_C",(PetscViewer,const char[],const char[],const char[],const char[]),(viewer,iname,jname,aname,cname));
385:   return(0);
386: }

388: static PetscErrorCode  PetscViewerHDF5GetAIJNames_HDF5(PetscViewer viewer, const char *iname[], const char *jname[], const char *aname[], const char *cname[])
389: {
390:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

393:   *iname = hdf5->mataij_iname;
394:   *jname = hdf5->mataij_jname;
395:   *aname = hdf5->mataij_aname;
396:   *cname = hdf5->mataij_cname;
397:   return(0);
398: }

400: /*@C
401:   PetscViewerHDF5GetAIJNames - Get the names of the datasets representing the three AIJ (CRS) arrays and the name of the attribute storing the number of columns within the HDF5 file.

403:   Collective on PetscViewer

405:   Input Parameters:
406: .  viewer - the PetscViewer; either ASCII or binary

408:   Output Parameters:
409: +  iname - name of dataset i representing row pointers; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix
410: .  jname - name of dataset j representing column indices
411: .  aname - name of dataset a representing matrix values
412: -  cname - name of attribute stoting column count

414:   Level: advanced

416:   Notes:
417:   Current defaults are (iname, jname, aname, cname) = ("jc", "ir", "data", "MATLAB_sparse") so that MAT files can be readily loaded.

419: .seealso: MatLoad(), PetscViewerCreate(), PetscViewerSetType(), PETSCVIEWERHDF5, PetscViewerHDF5SetAIJNames()
420: @*/
421: /* TODO Once corresponding MatView is implemented, add this:
422:   Current defaults are (iname, jname, aname, cname) = ("i", "j", "a", "ncols").
423:   For PetscViewerFormat PETSC_VIEWER_HDF5_MAT they are ("jc", "ir", "data", "MATLAB_sparse") so that MAT files can be loaded.
424: */
425: PetscErrorCode  PetscViewerHDF5GetAIJNames(PetscViewer viewer, const char *iname[], const char *jname[], const char *aname[], const char *cname[])
426: {

435:   PetscUseMethod(viewer,"PetscViewerHDF5GetAIJNames_C",(PetscViewer,const char*[],const char*[],const char*[],const char*[]),(viewer,iname,jname,aname,cname));
436:   return(0);
437: }

439: static PetscErrorCode PetscViewerSetUp_HDF5(PetscViewer viewer)
440: {
441:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
442:   PetscErrorCode   ierr;

445:   if (!hdf5->mataij_names_set) {
446: /* TODO Once corresponding MatView is implemented, uncomment */
447: #if 0
448:     if (viewer->format == PETSC_VIEWER_HDF5_MAT) {
449: #endif
450:       PetscViewerHDF5SetAIJNames_HDF5(viewer,"jc","ir","data","MATLAB_sparse");
451: #if 0
452:     } else {
453:       PetscViewerHDF5SetAIJNames_HDF5(viewer,"i","j","a","ncols");
454:     }
455: #endif
456:   }
457:   return(0);
458: }

460: /*MC
461:    PETSCVIEWERHDF5 - A viewer that writes to an HDF5 file


464: .seealso:  PetscViewerHDF5Open(), PetscViewerStringSPrintf(), PetscViewerSocketOpen(), PetscViewerDrawOpen(), PETSCVIEWERSOCKET,
465:            PetscViewerCreate(), PetscViewerASCIIOpen(), PetscViewerBinaryOpen(), PETSCVIEWERBINARY, PETSCVIEWERDRAW, PETSCVIEWERSTRING,
466:            PetscViewerMatlabOpen(), VecView(), DMView(), PetscViewerMatlabPutArray(), PETSCVIEWERASCII, PETSCVIEWERMATLAB,
467:            PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType()

469:   Level: beginner
470: M*/

472: PETSC_EXTERN PetscErrorCode PetscViewerCreate_HDF5(PetscViewer v)
473: {
474:   PetscViewer_HDF5 *hdf5;
475:   PetscErrorCode   ierr;

478:   PetscNewLog(v,&hdf5);

480:   v->data                = (void*) hdf5;
481:   v->ops->destroy        = PetscViewerDestroy_HDF5;
482:   v->ops->setfromoptions = PetscViewerSetFromOptions_HDF5;
483:   v->ops->setup          = PetscViewerSetUp_HDF5;
484:   v->ops->flush          = 0;
485:   hdf5->btype            = (PetscFileMode) -1;
486:   hdf5->filename         = 0;
487:   hdf5->timestep         = -1;
488:   hdf5->groups           = NULL;

490:   hdf5->mataij_iname     = NULL;
491:   hdf5->mataij_jname     = NULL;
492:   hdf5->mataij_aname     = NULL;
493:   hdf5->mataij_cname     = NULL;
494:   hdf5->mataij_names_set = PETSC_FALSE;

496:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_HDF5);
497:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_HDF5);
498:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_HDF5);
499:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_HDF5);
500:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerHDF5SetBaseDimension2_C",PetscViewerHDF5SetBaseDimension2_HDF5);
501:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerHDF5SetSPOutput_C",PetscViewerHDF5SetSPOutput_HDF5);
502:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerHDF5SetAIJNames_C",PetscViewerHDF5SetAIJNames_HDF5);
503:   PetscObjectComposeFunction((PetscObject)v,"PetscViewerHDF5GetAIJNames_C",PetscViewerHDF5GetAIJNames_HDF5);
504:   return(0);
505: }

507: /*@C
508:    PetscViewerHDF5Open - Opens a file for HDF5 input/output.

510:    Collective on MPI_Comm

512:    Input Parameters:
513: +  comm - MPI communicator
514: .  name - name of file
515: -  type - type of file
516: $    FILE_MODE_WRITE - create new file for binary output
517: $    FILE_MODE_READ - open existing file for binary input
518: $    FILE_MODE_APPEND - open existing file for binary output

520:    Output Parameter:
521: .  hdf5v - PetscViewer for HDF5 input/output to use with the specified file

523:   Options Database:
524: .  -viewer_hdf5_base_dimension2 - turns on (true) or off (false) using a dimension of 2 in the HDF5 file even if the bs/dof of the vector is 1
525: .  -viewer_hdf5_sp_output - forces (if true) the viewer to write data in single precision independent on the precision of PetscReal

527:    Level: beginner

529:    Note:
530:    This PetscViewer should be destroyed with PetscViewerDestroy().

532:    Concepts: HDF5 files
533:    Concepts: PetscViewerHDF5^creating

535: .seealso: PetscViewerASCIIOpen(), PetscViewerPushFormat(), PetscViewerDestroy(), PetscViewerHDF5SetBaseDimension2(),
536:           PetscViewerHDF5SetSPOutput(), PetscViewerHDF5GetBaseDimension2(), VecView(), MatView(), VecLoad(),
537:           MatLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName()
538: @*/
539: PetscErrorCode  PetscViewerHDF5Open(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *hdf5v)
540: {

544:   PetscViewerCreate(comm, hdf5v);
545:   PetscViewerSetType(*hdf5v, PETSCVIEWERHDF5);
546:   PetscViewerFileSetMode(*hdf5v, type);
547:   PetscViewerFileSetName(*hdf5v, name);
548:   return(0);
549: }

551: /*@C
552:   PetscViewerHDF5GetFileId - Retrieve the file id, this file ID then can be used in direct HDF5 calls

554:   Not collective

556:   Input Parameter:
557: . viewer - the PetscViewer

559:   Output Parameter:
560: . file_id - The file id

562:   Level: intermediate

564: .seealso: PetscViewerHDF5Open()
565: @*/
566: PetscErrorCode  PetscViewerHDF5GetFileId(PetscViewer viewer, hid_t *file_id)
567: {
568:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

572:   if (file_id) *file_id = hdf5->file_id;
573:   return(0);
574: }

576: /*@C
577:   PetscViewerHDF5PushGroup - Set the current HDF5 group for output

579:   Not collective

581:   Input Parameters:
582: + viewer - the PetscViewer
583: - name - The group name

585:   Level: intermediate

587: .seealso: PetscViewerHDF5Open(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup(),PetscViewerHDF5OpenGroup()
588: @*/
589: PetscErrorCode  PetscViewerHDF5PushGroup(PetscViewer viewer, const char *name)
590: {
591:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
592:   GroupList        *groupNode;
593:   PetscErrorCode   ierr;

598:   PetscNew(&groupNode);
599:   PetscStrallocpy(name, (char**) &groupNode->name);

601:   groupNode->next = hdf5->groups;
602:   hdf5->groups    = groupNode;
603:   return(0);
604: }

606: /*@
607:   PetscViewerHDF5PopGroup - Return the current HDF5 group for output to the previous value

609:   Not collective

611:   Input Parameter:
612: . viewer - the PetscViewer

614:   Level: intermediate

616: .seealso: PetscViewerHDF5Open(),PetscViewerHDF5PushGroup(),PetscViewerHDF5GetGroup(),PetscViewerHDF5OpenGroup()
617: @*/
618: PetscErrorCode  PetscViewerHDF5PopGroup(PetscViewer viewer)
619: {
620:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
621:   GroupList        *groupNode;
622:   PetscErrorCode   ierr;

626:   if (!hdf5->groups) SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONGSTATE, "HDF5 group stack is empty, cannot pop");
627:   groupNode    = hdf5->groups;
628:   hdf5->groups = hdf5->groups->next;
629:   PetscFree(groupNode->name);
630:   PetscFree(groupNode);
631:   return(0);
632: }

634: /*@C
635:   PetscViewerHDF5GetGroup - Get the current HDF5 group name (full path), set with PetscViewerHDF5PushGroup()/PetscViewerHDF5PopGroup().
636:   If none has been assigned, returns NULL.

638:   Not collective

640:   Input Parameter:
641: . viewer - the PetscViewer

643:   Output Parameter:
644: . name - The group name

646:   Level: intermediate

648: .seealso: PetscViewerHDF5Open(),PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5OpenGroup()
649: @*/
650: PetscErrorCode  PetscViewerHDF5GetGroup(PetscViewer viewer, const char **name)
651: {
652:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *) viewer->data;

657:   if (hdf5->groups) *name = hdf5->groups->name;
658:   else *name = NULL;
659:   return(0);
660: }

662: /*@
663:   PetscViewerHDF5OpenGroup - Open the HDF5 group with the name (full path) returned by PetscViewerHDF5GetGroup(),
664:   and return this group's ID and file ID.
665:   If PetscViewerHDF5GetGroup() yields NULL, then group ID is file ID.

667:   Not collective

669:   Input Parameter:
670: . viewer - the PetscViewer

672:   Output Parameter:
673: + fileId - The HDF5 file ID
674: - groupId - The HDF5 group ID

676:   Notes:
677:   If the viewer is writable, the group is created if it doesn't exist yet.

679:   Level: intermediate

681: .seealso: PetscViewerHDF5Open(),PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
682: @*/
683: PetscErrorCode PetscViewerHDF5OpenGroup(PetscViewer viewer, hid_t *fileId, hid_t *groupId)
684: {
685:   hid_t          file_id;
686:   H5O_type_t     type;
687:   const char     *groupName = NULL;
688:   PetscBool      create;

692:   PetscViewerWritable(viewer, &create);
693:   PetscViewerHDF5GetFileId(viewer, &file_id);
694:   PetscViewerHDF5GetGroup(viewer, &groupName);
695:   PetscViewerHDF5Traverse_Internal(viewer, groupName, create, NULL, &type);
696:   if (type != H5O_TYPE_GROUP) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Path %s resolves to something which is not a group", groupName);
697:   PetscStackCallHDF5Return(*groupId,H5Gopen2,(file_id, groupName ? groupName : "/", H5P_DEFAULT));
698:   *fileId  = file_id;
699:   return(0);
700: }

702: /*@
703:   PetscViewerHDF5IncrementTimestep - Increments the current timestep for the HDF5 output. Fields are stacked in time.

705:   Not collective

707:   Input Parameter:
708: . viewer - the PetscViewer

710:   Level: intermediate

712: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5SetTimestep(), PetscViewerHDF5GetTimestep()
713: @*/
714: PetscErrorCode PetscViewerHDF5IncrementTimestep(PetscViewer viewer)
715: {
716:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

720:   ++hdf5->timestep;
721:   return(0);
722: }

724: /*@
725:   PetscViewerHDF5SetTimestep - Set the current timestep for the HDF5 output. Fields are stacked in time. A timestep
726:   of -1 disables blocking with timesteps.

728:   Not collective

730:   Input Parameters:
731: + viewer - the PetscViewer
732: - timestep - The timestep number

734:   Level: intermediate

736: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5IncrementTimestep(), PetscViewerHDF5GetTimestep()
737: @*/
738: PetscErrorCode  PetscViewerHDF5SetTimestep(PetscViewer viewer, PetscInt timestep)
739: {
740:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

744:   hdf5->timestep = timestep;
745:   return(0);
746: }

748: /*@
749:   PetscViewerHDF5GetTimestep - Get the current timestep for the HDF5 output. Fields are stacked in time.

751:   Not collective

753:   Input Parameter:
754: . viewer - the PetscViewer

756:   Output Parameter:
757: . timestep - The timestep number

759:   Level: intermediate

761: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5IncrementTimestep(), PetscViewerHDF5SetTimestep()
762: @*/
763: PetscErrorCode  PetscViewerHDF5GetTimestep(PetscViewer viewer, PetscInt *timestep)
764: {
765:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;

770:   *timestep = hdf5->timestep;
771:   return(0);
772: }

774: /*@C
775:   PetscDataTypeToHDF5DataType - Converts the PETSc name of a datatype to its HDF5 name.

777:   Not collective

779:   Input Parameter:
780: . ptype - the PETSc datatype name (for example PETSC_DOUBLE)

782:   Output Parameter:
783: . mtype - the MPI datatype (for example MPI_DOUBLE, ...)

785:   Level: advanced

787: .seealso: PetscDataType, PetscHDF5DataTypeToPetscDataType()
788: @*/
789: PetscErrorCode PetscDataTypeToHDF5DataType(PetscDataType ptype, hid_t *htype)
790: {
792:   if (ptype == PETSC_INT)
793: #if defined(PETSC_USE_64BIT_INDICES)
794:                                        *htype = H5T_NATIVE_LLONG;
795: #else
796:                                        *htype = H5T_NATIVE_INT;
797: #endif
798:   else if (ptype == PETSC_DOUBLE)      *htype = H5T_NATIVE_DOUBLE;
799:   else if (ptype == PETSC_LONG)        *htype = H5T_NATIVE_LONG;
800:   else if (ptype == PETSC_SHORT)       *htype = H5T_NATIVE_SHORT;
801:   else if (ptype == PETSC_ENUM)        *htype = H5T_NATIVE_DOUBLE;
802:   else if (ptype == PETSC_BOOL)        *htype = H5T_NATIVE_INT;
803:   else if (ptype == PETSC_FLOAT)       *htype = H5T_NATIVE_FLOAT;
804:   else if (ptype == PETSC_CHAR)        *htype = H5T_NATIVE_CHAR;
805:   else if (ptype == PETSC_BIT_LOGICAL) *htype = H5T_NATIVE_UCHAR;
806:   else if (ptype == PETSC_STRING)      *htype = H5Tcopy(H5T_C_S1);
807:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported PETSc datatype");
808:   return(0);
809: }

811: /*@C
812:   PetscHDF5DataTypeToPetscDataType - Finds the PETSc name of a datatype from its HDF5 name

814:   Not collective

816:   Input Parameter:
817: . htype - the HDF5 datatype (for example H5T_NATIVE_DOUBLE, ...)

819:   Output Parameter:
820: . ptype - the PETSc datatype name (for example PETSC_DOUBLE)

822:   Level: advanced

824: .seealso: PetscDataType, PetscHDF5DataTypeToPetscDataType()
825: @*/
826: PetscErrorCode PetscHDF5DataTypeToPetscDataType(hid_t htype, PetscDataType *ptype)
827: {
829: #if defined(PETSC_USE_64BIT_INDICES)
830:   if      (htype == H5T_NATIVE_INT)    *ptype = PETSC_LONG;
831:   else if (htype == H5T_NATIVE_LLONG)  *ptype = PETSC_INT;
832: #else
833:   if      (htype == H5T_NATIVE_INT)    *ptype = PETSC_INT;
834: #endif
835:   else if (htype == H5T_NATIVE_DOUBLE) *ptype = PETSC_DOUBLE;
836:   else if (htype == H5T_NATIVE_LONG)   *ptype = PETSC_LONG;
837:   else if (htype == H5T_NATIVE_SHORT)  *ptype = PETSC_SHORT;
838:   else if (htype == H5T_NATIVE_FLOAT)  *ptype = PETSC_FLOAT;
839:   else if (htype == H5T_NATIVE_CHAR)   *ptype = PETSC_CHAR;
840:   else if (htype == H5T_NATIVE_UCHAR)  *ptype = PETSC_CHAR;
841:   else if (htype == H5T_C_S1)          *ptype = PETSC_STRING;
842:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported HDF5 datatype");
843:   return(0);
844: }

846: /*@C
847:  PetscViewerHDF5WriteAttribute - Write an attribute

849:   Input Parameters:
850: + viewer - The HDF5 viewer
851: . dataset - The parent dataset name, relative to the current group. NULL means a group-wise attribute.
852: . name   - The attribute name
853: . datatype - The attribute type
854: - value    - The attribute value

856:   Level: advanced

858: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5WriteObjectAttribute(), PetscViewerHDF5ReadAttribute(), PetscViewerHDF5HasAttribute(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
859: @*/
860: PetscErrorCode PetscViewerHDF5WriteAttribute(PetscViewer viewer, const char dataset[], const char name[], PetscDataType datatype, const void *value)
861: {
862:   char           *parent;
863:   hid_t          h5, dataspace, obj, attribute, dtype;
864:   PetscBool      has;

872:   PetscViewerHDF5GetAbsolutePath_Internal(viewer, dataset, &parent);
873:   PetscViewerHDF5Traverse_Internal(viewer, parent, PETSC_TRUE, NULL, NULL);
874:   PetscViewerHDF5HasAttribute_Internal(viewer, parent, name, &has);
875:   PetscDataTypeToHDF5DataType(datatype, &dtype);
876:   if (datatype == PETSC_STRING) {
877:     size_t len;
878:     PetscStrlen((const char *) value, &len);
879:     PetscStackCallHDF5(H5Tset_size,(dtype, len+1));
880:   }
881:   PetscViewerHDF5GetFileId(viewer, &h5);
882:   PetscStackCallHDF5Return(dataspace,H5Screate,(H5S_SCALAR));
883:   PetscStackCallHDF5Return(obj,H5Oopen,(h5, parent, H5P_DEFAULT));
884:   if (has) {
885:     PetscStackCallHDF5Return(attribute,H5Aopen_name,(obj, name));
886:   } else {
887:     PetscStackCallHDF5Return(attribute,H5Acreate2,(obj, name, dtype, dataspace, H5P_DEFAULT, H5P_DEFAULT));
888:   }
889:   PetscStackCallHDF5(H5Awrite,(attribute, dtype, value));
890:   if (datatype == PETSC_STRING) PetscStackCallHDF5(H5Tclose,(dtype));
891:   PetscStackCallHDF5(H5Aclose,(attribute));
892:   PetscStackCallHDF5(H5Oclose,(obj));
893:   PetscStackCallHDF5(H5Sclose,(dataspace));
894:   PetscFree(parent);
895:   return(0);
896: }

898: /*@C
899:  PetscViewerHDF5WriteObjectAttribute - Write an attribute to the dataset matching the given PetscObject by name

901:   Input Parameters:
902: + viewer   - The HDF5 viewer
903: . obj      - The object whose name is used to lookup the parent dataset, relative to the current group.
904: . name     - The attribute name
905: . datatype - The attribute type
906: - value    - The attribute value

908:   Notes:
909:   This fails if current_group/object_name doesn't resolve to a dataset (the path doesn't exist or is not a dataset).
910:   You might want to check first if it does using PetscViewerHDF5HasObject().

912:   Level: advanced

914: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5WriteAttribute(), PetscViewerHDF5ReadObjectAttribute(), PetscViewerHDF5HasObjectAttribute(), PetscViewerHDF5HasObject(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
915: @*/
916: PetscErrorCode PetscViewerHDF5WriteObjectAttribute(PetscViewer viewer, PetscObject obj, const char name[], PetscDataType datatype, const void *value)
917: {

925:   PetscViewerHDF5CheckNamedObject_Internal(viewer, obj);
926:   PetscViewerHDF5WriteAttribute(viewer, obj->name, name, datatype, value);
927:   return(0);
928: }

930: /*@C
931:  PetscViewerHDF5ReadAttribute - Read an attribute

933:   Input Parameters:
934: + viewer - The HDF5 viewer
935: . dataset - The parent dataset name, relative to the current group. NULL means a group-wise attribute.
936: . name   - The attribute name
937: - datatype - The attribute type

939:   Output Parameter:
940: . value    - The attribute value

942:   Level: advanced

944: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5ReadObjectAttribute(), PetscViewerHDF5WriteAttribute(), PetscViewerHDF5HasAttribute(), PetscViewerHDF5HasObject(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
945: @*/
946: PetscErrorCode PetscViewerHDF5ReadAttribute(PetscViewer viewer, const char dataset[], const char name[], PetscDataType datatype, void *value)
947: {
948:   char           *parent;
949:   hid_t          h5, obj, attribute, atype, dtype;
950:   PetscBool      has;

958:   PetscViewerHDF5GetAbsolutePath_Internal(viewer, dataset, &parent);
959:   PetscViewerHDF5Traverse_Internal(viewer, parent, PETSC_FALSE, &has, NULL);
960:   if (has) {PetscViewerHDF5HasAttribute_Internal(viewer, parent, name, &has);}
961:   if (!has) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Attribute %s/%s does not exist", parent, name);
962:   PetscDataTypeToHDF5DataType(datatype, &dtype);
963:   PetscViewerHDF5GetFileId(viewer, &h5);
964:   PetscStackCallHDF5Return(obj,H5Oopen,(h5, parent, H5P_DEFAULT));
965:   PetscStackCallHDF5Return(attribute,H5Aopen_name,(obj, name));
966:   if (datatype == PETSC_STRING) {
967:     size_t len;
968:     PetscStackCallHDF5Return(atype,H5Aget_type,(attribute));
969:     PetscStackCall("H5Tget_size",len = H5Tget_size(atype));
970:     PetscStackCallHDF5(H5Tclose,(atype));
971:     PetscMalloc((len+1) * sizeof(char *), &value);
972:   }
973:   PetscStackCallHDF5(H5Aread,(attribute, dtype, value));
974:   PetscStackCallHDF5(H5Aclose,(attribute));
975:   /* H5Oclose can be used to close groups, datasets, or committed datatypes */
976:   PetscStackCallHDF5(H5Oclose,(obj));
977:   PetscFree(parent);
978:   return(0);
979: }

981: /*@C
982:  PetscViewerHDF5ReadObjectAttribute - Read an attribute from the dataset matching the given PetscObject by name

984:   Input Parameters:
985: + viewer   - The HDF5 viewer
986: . obj      - The object whose name is used to lookup the parent dataset, relative to the current group.
987: . name     - The attribute name
988: - datatype - The attribute type

990:   Output Parameter:
991: . value    - The attribute value

993:   Notes:
994:   This fails if current_group/object_name doesn't resolve to a dataset (the path doesn't exist or is not a dataset).
995:   You might want to check first if it does using PetscViewerHDF5HasObject().

997:   Level: advanced

999: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5ReadAttribute() PetscViewerHDF5WriteObjectAttribute(), PetscViewerHDF5HasObjectAttribute(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
1000: @*/
1001: PetscErrorCode PetscViewerHDF5ReadObjectAttribute(PetscViewer viewer, PetscObject obj, const char name[], PetscDataType datatype, void *value)
1002: {

1010:   PetscViewerHDF5CheckNamedObject_Internal(viewer, obj);
1011:   PetscViewerHDF5ReadAttribute(viewer, obj->name, name, datatype, value);
1012:   return(0);
1013: }

1015: PETSC_STATIC_INLINE PetscErrorCode PetscViewerHDF5Traverse_Inner_Internal(hid_t h5, const char name[], PetscBool createGroup, PetscBool *exists_)
1016: {
1017:   htri_t exists;
1018:   hid_t group;

1021:   PetscStackCallHDF5Return(exists,H5Lexists,(h5, name, H5P_DEFAULT));
1022:   if (exists) PetscStackCallHDF5Return(exists,H5Oexists_by_name,(h5, name, H5P_DEFAULT));
1023:   if (!exists && createGroup) {
1024:     PetscStackCallHDF5Return(group,H5Gcreate2,(h5, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
1025:     PetscStackCallHDF5(H5Gclose,(group));
1026:     exists = PETSC_TRUE;
1027:   }
1028:   *exists_ = (PetscBool) exists;
1029:   return(0);
1030: }

1032: static PetscErrorCode PetscViewerHDF5Traverse_Internal(PetscViewer viewer, const char name[], PetscBool createGroup, PetscBool *has, H5O_type_t *otype)
1033: {
1034:   const char     rootGroupName[] = "/";
1035:   hid_t          h5;
1036:   PetscBool      exists=PETSC_FALSE;
1037:   PetscInt       i;
1038:   int            n;
1039:   char           **hierarchy;
1040:   char           buf[PETSC_MAX_PATH_LEN]="";

1046:   else name = rootGroupName;
1047:   if (has) {
1049:     *has = PETSC_FALSE;
1050:   }
1051:   if (otype) {
1053:     *otype = H5O_TYPE_UNKNOWN;
1054:   }
1055:   PetscViewerHDF5GetFileId(viewer, &h5);

1057:   /*
1058:      Unfortunately, H5Oexists_by_name() fails if any object in hierarchy is missing.
1059:      Hence, each of them needs to be tested separately:
1060:      1) whether it's a valid link
1061:      2) whether this link resolves to an object
1062:      See H5Oexists_by_name() documentation.
1063:   */
1064:   PetscStrToArray(name,'/',&n,&hierarchy);
1065:   if (!n) {
1066:     /*  Assume group "/" always exists in accordance with HDF5 >= 1.10.0. See H5Lexists() documentation. */
1067:     if (has)   *has   = PETSC_TRUE;
1068:     if (otype) *otype = H5O_TYPE_GROUP;
1069:     PetscStrToArrayDestroy(n,hierarchy);
1070:     return(0);
1071:   }
1072:   for (i=0; i<n; i++) {
1073:     PetscStrcat(buf,"/");
1074:     PetscStrcat(buf,hierarchy[i]);
1075:     PetscViewerHDF5Traverse_Inner_Internal(h5, buf, createGroup, &exists);
1076:     if (!exists) break;
1077:   }
1078:   PetscStrToArrayDestroy(n,hierarchy);

1080:   /* If the object exists, get its type */
1081:   if (exists && otype) {
1082:     H5O_info_t info;

1084:     /* We could use H5Iget_type() here but that would require opening the object. This way we only need its name. */
1085:     PetscStackCallHDF5(H5Oget_info_by_name,(h5, name, &info, H5P_DEFAULT));
1086:     *otype = info.type;
1087:   }
1088:   if (has) *has = exists;
1089:   return(0);
1090: }

1092: /*@
1093:  PetscViewerHDF5HasGroup - Check whether the current (pushed) group exists in the HDF5 file

1095:   Input Parameters:
1096: . viewer - The HDF5 viewer

1098:   Output Parameter:
1099: . has    - Flag for group existence

1101:   Notes:
1102:   If the path exists but is not a group, this returns PETSC_FALSE as well.

1104:   Level: advanced

1106: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5PushGroup(), PetscViewerHDF5PopGroup(), PetscViewerHDF5OpenGroup()
1107: @*/
1108: PetscErrorCode PetscViewerHDF5HasGroup(PetscViewer viewer, PetscBool *has)
1109: {
1110:   H5O_type_t type;
1111:   const char *name;

1117:   PetscViewerHDF5GetGroup(viewer, &name);
1118:   PetscViewerHDF5Traverse_Internal(viewer, name, PETSC_FALSE, has, &type);
1119:   *has = (type == H5O_TYPE_GROUP) ? PETSC_TRUE : PETSC_FALSE;
1120:   return(0);
1121: }

1123: /*@
1124:  PetscViewerHDF5HasObject - Check whether a dataset with the same name as given object exists in the HDF5 file under current group

1126:   Input Parameters:
1127: + viewer - The HDF5 viewer
1128: - obj    - The named object

1130:   Output Parameter:
1131: . has    - Flag for dataset existence; PETSC_FALSE for unnamed object

1133:   Notes:
1134:   If the path exists but is not a dataset, this returns PETSC_FALSE as well.

1136:   Level: advanced

1138: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5HasAttribute(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
1139: @*/
1140: PetscErrorCode PetscViewerHDF5HasObject(PetscViewer viewer, PetscObject obj, PetscBool *has)
1141: {
1142:   H5O_type_t type;
1143:   char *path;

1150:   *has = PETSC_FALSE;
1151:   if (!obj->name) return(0);
1152:   PetscViewerHDF5GetAbsolutePath_Internal(viewer, obj->name, &path);
1153:   PetscViewerHDF5Traverse_Internal(viewer, path, PETSC_FALSE, has, &type);
1154:   *has = (type == H5O_TYPE_DATASET) ? PETSC_TRUE : PETSC_FALSE;
1155:   PetscFree(path);
1156:   return(0);
1157: }

1159: /*@C
1160:  PetscViewerHDF5HasAttribute - Check whether an attribute exists

1162:   Input Parameters:
1163: + viewer - The HDF5 viewer
1164: . dataset - The parent dataset name, relative to the current group. NULL means a group-wise attribute.
1165: - name   - The attribute name

1167:   Output Parameter:
1168: . has    - Flag for attribute existence

1170:   Level: advanced

1172: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5HasObjectAttribute(), PetscViewerHDF5WriteAttribute(), PetscViewerHDF5ReadAttribute(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
1173: @*/
1174: PetscErrorCode PetscViewerHDF5HasAttribute(PetscViewer viewer, const char dataset[], const char name[], PetscBool *has)
1175: {
1176:   char           *parent;

1184:   PetscViewerHDF5GetAbsolutePath_Internal(viewer, dataset, &parent);
1185:   PetscViewerHDF5Traverse_Internal(viewer, parent, PETSC_FALSE, has, NULL);
1186:   if (*has) {PetscViewerHDF5HasAttribute_Internal(viewer, parent, name, has);}
1187:   PetscFree(parent);
1188:   return(0);
1189: }

1191: /*@C
1192:  PetscViewerHDF5HasObjectAttribute - Check whether an attribute is attached to the dataset matching the given PetscObject by name

1194:   Input Parameters:
1195: + viewer - The HDF5 viewer
1196: . obj    - The object whose name is used to lookup the parent dataset, relative to the current group.
1197: - name   - The attribute name

1199:   Output Parameter:
1200: . has    - Flag for attribute existence

1202:   Notes:
1203:   This fails if current_group/object_name doesn't resolve to a dataset (the path doesn't exist or is not a dataset).
1204:   You might want to check first if it does using PetscViewerHDF5HasObject().

1206:   Level: advanced

1208: .seealso: PetscViewerHDF5Open(), PetscViewerHDF5HasAttribute(), PetscViewerHDF5WriteObjectAttribute(), PetscViewerHDF5ReadObjectAttribute(), PetscViewerHDF5HasObject(), PetscViewerHDF5PushGroup(),PetscViewerHDF5PopGroup(),PetscViewerHDF5GetGroup()
1209: @*/
1210: PetscErrorCode PetscViewerHDF5HasObjectAttribute(PetscViewer viewer, PetscObject obj, const char name[], PetscBool *has)
1211: {

1219:   PetscViewerHDF5CheckNamedObject_Internal(viewer, obj);
1220:   PetscViewerHDF5HasAttribute(viewer, obj->name, name, has);
1221:   return(0);
1222: }

1224: static PetscErrorCode PetscViewerHDF5HasAttribute_Internal(PetscViewer viewer, const char parent[], const char name[], PetscBool *has)
1225: {
1226:   hid_t          h5;
1227:   htri_t         hhas;

1231:   PetscViewerHDF5GetFileId(viewer, &h5);
1232:   PetscStackCallHDF5Return(hhas,H5Aexists_by_name,(h5, parent, name, H5P_DEFAULT));
1233:   *has = hhas ? PETSC_TRUE : PETSC_FALSE;
1234:   return(0);
1235: }

1237: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
1238: {
1239:   HDF5ReadCtx    h=NULL;

1243:   PetscNew(&h);
1244:   PetscViewerHDF5OpenGroup(viewer, &h->file, &h->group);
1245:   PetscStackCallHDF5Return(h->dataset,H5Dopen2,(h->group, name, H5P_DEFAULT));
1246:   PetscStackCallHDF5Return(h->dataspace,H5Dget_space,(h->dataset));
1247:   PetscViewerHDF5GetTimestep(viewer, &h->timestep);
1248:   PetscViewerHDF5HasAttribute(viewer,name,"complex",&h->complexVal);
1249:   if (h->complexVal) {PetscViewerHDF5ReadAttribute(viewer,name,"complex",PETSC_BOOL,&h->complexVal);}
1250:   /* MATLAB stores column vectors horizontally */
1251:   PetscViewerHDF5HasAttribute(viewer,name,"MATLAB_class",&h->horizontal);
1252:   /* Create property list for collective dataset read */
1253:   PetscStackCallHDF5Return(h->plist,H5Pcreate,(H5P_DATASET_XFER));
1254: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
1255:   PetscStackCallHDF5(H5Pset_dxpl_mpio,(h->plist, H5FD_MPIO_COLLECTIVE));
1256: #endif
1257:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
1258:   *ctx = h;
1259:   return(0);
1260: }

1262: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
1263: {
1264:   HDF5ReadCtx    h;

1268:   h = *ctx;
1269:   PetscStackCallHDF5(H5Pclose,(h->plist));
1270:   PetscStackCallHDF5(H5Gclose,(h->group));
1271:   PetscStackCallHDF5(H5Sclose,(h->dataspace));
1272:   PetscStackCallHDF5(H5Dclose,(h->dataset));
1273:   PetscFree(*ctx);
1274:   return(0);
1275: }

1277: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout *map_)
1278: {
1279:   int            rdim, dim;
1280:   hsize_t        dims[4];
1281:   PetscInt       bsInd, lenInd, bs, len, N;
1282:   PetscLayout    map;

1286:   if (!(*map_)) {
1287:     PetscLayoutCreate(PetscObjectComm((PetscObject)viewer),map_);
1288:   }
1289:   map = *map_;
1290:   /* calculate expected number of dimensions */
1291:   dim = 0;
1292:   if (ctx->timestep >= 0) ++dim;
1293:   ++dim; /* length in blocks */
1294:   if (ctx->complexVal) ++dim;
1295:   /* get actual number of dimensions in dataset */
1296:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(ctx->dataspace, dims, NULL));
1297:   /* calculate expected dimension indices */
1298:   lenInd = 0;
1299:   if (ctx->timestep >= 0) ++lenInd;
1300:   bsInd = lenInd + 1;
1301:   ctx->dim2 = PETSC_FALSE;
1302:   if (rdim == dim) {
1303:     bs = 1; /* support vectors stored as 1D array */
1304:   } else if (rdim == dim+1) {
1305:     bs = (PetscInt) dims[bsInd];
1306:     if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
1307:   } else {
1308:     SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Number of dimensions %d not %d as expected", rdim, dim);
1309:   }
1310:   len = dims[lenInd];
1311:   if (ctx->horizontal) {
1312:     if (len != 1) SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Cannot have horizontal array with number of rows > 1. In case of MATLAB MAT-file, vectors must be saved as column vectors.");
1313:     len = bs;
1314:     bs = 1;
1315:   }
1316:   N = (PetscInt) len*bs;

1318:   /* Set Vec sizes,blocksize,and type if not already set */
1319:   if (map->bs < 0) {
1320:     PetscLayoutSetBlockSize(map, bs);
1321:   } else if (map->bs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Block size of array in file is %D, not %D as expected",bs,map->bs);
1322:   if (map->N < 0) {
1323:     PetscLayoutSetSize(map, N);
1324:   } else if (map->N != N) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Global size of array in file is %D, not %D as expected",N,map->N);
1325:   return(0);
1326: }

1328: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
1329: {
1330:   hsize_t        count[4], offset[4];
1331:   int            dim;
1332:   PetscInt       bs, n, low;

1336:   /* Compute local size and ownership range */
1337:   PetscLayoutSetUp(map);
1338:   PetscLayoutGetBlockSize(map, &bs);
1339:   PetscLayoutGetLocalSize(map, &n);
1340:   PetscLayoutGetRange(map, &low, NULL);

1342:   /* Each process defines a dataset and reads it from the hyperslab in the file */
1343:   dim  = 0;
1344:   if (ctx->timestep >= 0) {
1345:     count[dim]  = 1;
1346:     offset[dim] = ctx->timestep;
1347:     ++dim;
1348:   }
1349:   if (ctx->horizontal) {
1350:     count[dim]  = 1;
1351:     offset[dim] = 0;
1352:     ++dim;
1353:   }
1354:   {
1355:     PetscHDF5IntCast(n/bs, &count[dim]);
1356:     PetscHDF5IntCast(low/bs, &offset[dim]);
1357:     ++dim;
1358:   }
1359:   if (bs > 1 || ctx->dim2) {
1360:     if (PetscUnlikely(ctx->horizontal)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "cannot have horizontal array with blocksize > 1");
1361:     count[dim]  = bs;
1362:     offset[dim] = 0;
1363:     ++dim;
1364:   }
1365:   if (ctx->complexVal) {
1366:     count[dim]  = 2;
1367:     offset[dim] = 0;
1368:     ++dim;
1369:   }
1370:   PetscStackCallHDF5Return(*memspace,H5Screate_simple,(dim, count, NULL));
1371:   PetscStackCallHDF5(H5Sselect_hyperslab,(ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
1372:   return(0);
1373: }

1375: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
1376: {
1378:   PetscStackCallHDF5(H5Dread,(h->dataset, datatype, memspace, h->dataspace, h->plist, arr));
1379:   return(0);
1380: }

1382: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
1383: {
1384:   HDF5ReadCtx     h=NULL;
1385:   hid_t           memspace=0;
1386:   size_t          unitsize;
1387:   void            *arr;
1388:   PetscErrorCode  ierr;

1391:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
1392: #if defined(PETSC_USE_COMPLEX)
1393:   if (!h->complexVal) {
1394:     H5T_class_t clazz = H5Tget_class(datatype);
1395:     if (clazz == H5T_FLOAT) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"File contains real numbers but PETSc is configured for complex. The conversion is not yet implemented. Configure with --with-scalar-type=real.");
1396:   }
1397: #else
1398:   if (h->complexVal) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"File contains complex numbers but PETSc not configured for them. Configure with --with-scalar-type=complex.");
1399: #endif

1401:   PetscViewerHDF5ReadSizes_Private(viewer, h, &map);
1402:   PetscLayoutSetUp(map);
1403:   PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);

1405:   unitsize = H5Tget_size(datatype);
1406:   if (h->complexVal) unitsize *= 2;
1407:   if (unitsize <= 0 || unitsize > PetscMax(sizeof(PetscInt),sizeof(PetscScalar))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Sanity check failed: HDF5 function H5Tget_size(datatype) returned suspicious value %D",unitsize);
1408:   PetscMalloc(map->n*unitsize, &arr);

1410:   PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
1411:   PetscStackCallHDF5(H5Sclose,(memspace));
1412:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
1413:   *newarr = arr;
1414:   return(0);
1415: }

1417: /*@C
1418:  PetscViewerHDF5ReadSizes - Read block size and global size of a vector (Vec or IS) stored in an HDF5 file.

1420:   Input Parameters:
1421: + viewer - The HDF5 viewer
1422: - name   - The vector name

1424:   Output Parameter:
1425: + bs     - block size
1426: - N      - global size

1428:   Note:
1429:   A vector is stored as an HDF5 dataspace with 1-4 dimensions in this order:
1430:   1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).

1432:   A vectors can be stored as a 2D dataspace even if its blocksize is 1; see PetscViewerHDF5SetBaseDimension2().

1434:   Level: advanced

1436: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
1437: @*/
1438: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
1439: {
1440:   HDF5ReadCtx    h=NULL;
1441:   PetscLayout    map=NULL;

1446:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
1447:   PetscViewerHDF5ReadSizes_Private(viewer, h, &map);
1448:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
1449:   if (bs) *bs = map->bs;
1450:   if (N) *N = map->N;
1451:   PetscLayoutDestroy(&map);
1452:   return(0);
1453: }

1455: /*
1456:   The variable Petsc_Viewer_HDF5_keyval is used to indicate an MPI attribute that
1457:   is attached to a communicator, in this case the attribute is a PetscViewer.
1458: */
1459: PetscMPIInt Petsc_Viewer_HDF5_keyval = MPI_KEYVAL_INVALID;

1461: /*@C
1462:   PETSC_VIEWER_HDF5_ - Creates an HDF5 PetscViewer shared by all processors in a communicator.

1464:   Collective on MPI_Comm

1466:   Input Parameter:
1467: . comm - the MPI communicator to share the HDF5 PetscViewer

1469:   Level: intermediate

1471:   Options Database Keys:
1472: . -viewer_hdf5_filename <name>

1474:   Environmental variables:
1475: . PETSC_VIEWER_HDF5_FILENAME

1477:   Notes:
1478:   Unlike almost all other PETSc routines, PETSC_VIEWER_HDF5_ does not return
1479:   an error code.  The HDF5 PetscViewer is usually used in the form
1480: $       XXXView(XXX object, PETSC_VIEWER_HDF5_(comm));

1482: .seealso: PetscViewerHDF5Open(), PetscViewerCreate(), PetscViewerDestroy()
1483: @*/
1484: PetscViewer PETSC_VIEWER_HDF5_(MPI_Comm comm)
1485: {
1487:   PetscBool      flg;
1488:   PetscViewer    viewer;
1489:   char           fname[PETSC_MAX_PATH_LEN];
1490:   MPI_Comm       ncomm;

1493:   PetscCommDuplicate(comm,&ncomm,NULL);if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1494:   if (Petsc_Viewer_HDF5_keyval == MPI_KEYVAL_INVALID) {
1495:     MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN,MPI_COMM_NULL_DELETE_FN,&Petsc_Viewer_HDF5_keyval,0);
1496:     if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1497:   }
1498:   MPI_Comm_get_attr(ncomm,Petsc_Viewer_HDF5_keyval,(void**)&viewer,(int*)&flg);
1499:   if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1500:   if (!flg) { /* PetscViewer not yet created */
1501:     PetscOptionsGetenv(ncomm,"PETSC_VIEWER_HDF5_FILENAME",fname,PETSC_MAX_PATH_LEN,&flg);
1502:     if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1503:     if (!flg) {
1504:       PetscStrcpy(fname,"output.h5");
1505:       if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1506:     }
1507:     PetscViewerHDF5Open(ncomm,fname,FILE_MODE_WRITE,&viewer);
1508:     if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1509:     PetscObjectRegisterDestroy((PetscObject)viewer);
1510:     if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1511:     MPI_Comm_set_attr(ncomm,Petsc_Viewer_HDF5_keyval,(void*)viewer);
1512:     if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1513:   }
1514:   PetscCommDestroy(&ncomm);
1515:   if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_HDF5_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");return(0);}
1516:   PetscFunctionReturn(viewer);
1517: }