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> 9 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; 20 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; 35 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 46 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)); 51 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; 61 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)); 65 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)); 73 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 } 89 90 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 91 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 92 PetscFunctionReturn(PETSC_SUCCESS); 93 } 94