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
MatLoad_Dense_HDF5(Mat mat,PetscViewer viewer)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