Actual source code: densehdf5.c

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

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

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

     36:   PetscCheck(((PetscObject)mat)->name, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Mat name must be set with PetscObjectSetName() before MatLoad()");
     37: #if defined(PETSC_USE_REAL_SINGLE)
     38:   scalartype = H5T_NATIVE_FLOAT;
     39: #elif defined(PETSC_USE_REAL___FLOAT128)
     40:   #error "HDF5 output with 128 bit floats not supported."
     41: #elif defined(PETSC_USE_REAL___FP16)
     42:   #error "HDF5 output with 16 bit floats not supported."
     43: #else
     44:   scalartype = H5T_NATIVE_DOUBLE;
     45: #endif

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

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

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

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

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

     90:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
     91:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }