Actual source code: densehdf5.c


  2: /* TODO change to
  3: #include <../src/mat/impls/dense/seq/dense.h>
  4: */
  5: #include <../src/mat/impls/dense/mpi/mpidense.h>
  6: #include <petsc/private/isimpl.h>
  7: #include <petsc/private/vecimpl.h>
  8: #include <petsc/private/viewerhdf5impl.h>
  9: #include <petsclayouthdf5.h>

 11: #if defined(PETSC_HAVE_HDF5)
 12: PetscErrorCode MatLoad_Dense_HDF5(Mat mat, PetscViewer viewer)
 13: {
 14:   PetscViewer_HDF5 *hdf5;
 15:   hid_t             scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
 16:   PetscLayout       vmap;
 17:   PetscViewerFormat format;
 18:   PetscScalar      *a        = NULL;
 19:   const char       *mat_name = NULL;
 20:   MPI_Comm          comm;
 21:   PetscMPIInt       rank, size;

 23:   PetscViewerGetFormat(viewer, &format);
 24:   switch (format) {
 25:   case PETSC_VIEWER_HDF5_PETSC:
 26:   case PETSC_VIEWER_DEFAULT:
 27:   case PETSC_VIEWER_NATIVE:
 28:   case PETSC_VIEWER_HDF5_MAT:
 29:     break;
 30:   default:
 31:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 input.", PetscViewerFormats[format]);
 32:   }
 33:   hdf5 = (PetscViewer_HDF5 *)viewer->data;
 34:   /* we store dense matrix columns as blocks, like MATLAB save(filename,variables,'-v7.3') does */
 35:   hdf5->horizontal = PETSC_TRUE;

 38:   #if defined(PETSC_USE_REAL_SINGLE)
 39:   scalartype = H5T_NATIVE_FLOAT;
 40:   #elif defined(PETSC_USE_REAL___FLOAT128)
 41:     #error "HDF5 output with 128 bit floats not supported."
 42:   #elif defined(PETSC_USE_REAL___FP16)
 43:     #error "HDF5 output with 16 bit floats not supported."
 44:   #else
 45:   scalartype = H5T_NATIVE_DOUBLE;
 46:   #endif

 48:   PetscObjectGetComm((PetscObject)mat, &comm);
 49:   MPI_Comm_rank(comm, &rank);
 50:   MPI_Comm_size(comm, &size);
 51:   PetscObjectGetName((PetscObject)mat, &mat_name);

 53:   /* Convert user-defined rmap and cmap to the dataset layout */
 54:   PetscLayoutCreate(PetscObjectComm((PetscObject)mat), &vmap);
 55:   if (mat->rmap->n >= 0 && mat->cmap->N < 0) {
 56:     /* We need to know mat->cmap->N if user specifies custom mat->rmap->n, otherwise the latter would get ignored below */
 57:     PetscViewerHDF5ReadSizes(viewer, mat_name, &mat->cmap->N, NULL);
 58:   }
 59:   vmap->bs = mat->cmap->N;
 60:   vmap->n  = (mat->rmap->n < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->n * mat->cmap->N;
 61:   vmap->N  = (mat->rmap->N < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->N * mat->cmap->N;

 63:   /* Read the dataset and setup its layout */
 64:   /* Note: PetscViewerHDF5ReadSizes_Private takes into account that the dataset is transposed for MATLAB MAT files */
 65:   PetscViewerHDF5Load(viewer, mat_name, vmap, scalartype, (void **)&a);

 67:   /* Convert the dataset layout back to rmap and cmap */
 68:   mat->cmap->N = vmap->bs;
 69:   mat->rmap->n = vmap->n / mat->cmap->N;
 70:   mat->rmap->N = vmap->N / mat->cmap->N;
 71:   PetscLayoutSetUp(mat->rmap);
 72:   PetscLayoutSetUp(mat->cmap);
 73:   PetscLayoutDestroy(&vmap);

 75:   /* TODO adding PetscCopyMode flag to MatSeqDenseSetPreallocation would make this code cleaner and simpler */
 76:   {
 77:     PetscBool     flg;
 78:     Mat_SeqDense *impl;
 79:     PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
 80:     if (flg) {
 81:       impl = (Mat_SeqDense *)mat->data;
 82:       MatSeqDenseSetPreallocation(mat, a);
 83:     } else {
 84:       Mat_MPIDense *implm = (Mat_MPIDense *)mat->data;
 85:       MatMPIDenseSetPreallocation(mat, a);
 86:       impl = (Mat_SeqDense *)implm->A->data;
 87:     }
 88:     impl->user_alloc = PETSC_FALSE;
 89:   }

 91:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
 92:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
 93:   return 0;
 94: }
 95: #endif