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: #if defined(PETSC_HAVE_HDF5)
11: PetscErrorCode MatLoad_Dense_HDF5(Mat mat, PetscViewer viewer)
12: {
13: PetscViewer_HDF5 *hdf5;
14: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
15: PetscLayout vmap;
16: PetscViewerFormat format;
17: PetscScalar *a = NULL;
18: const char *mat_name = NULL;
19: MPI_Comm comm;
20: PetscMPIInt rank, size;
22: PetscFunctionBegin;
23: PetscCall(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;
37: PetscCheck(((PetscObject)mat)->name, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Mat name must be set with PetscObjectSetName() before MatLoad()");
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: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
49: PetscCallMPI(MPI_Comm_rank(comm, &rank));
50: PetscCallMPI(MPI_Comm_size(comm, &size));
51: PetscCall(PetscObjectGetName((PetscObject)mat, &mat_name));
53: /* Convert user-defined rmap and cmap to the dataset layout */
54: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(PetscLayoutSetUp(mat->rmap));
72: PetscCall(PetscLayoutSetUp(mat->cmap));
73: PetscCall(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: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
80: if (flg) {
81: impl = (Mat_SeqDense *)mat->data;
82: PetscCall(MatSeqDenseSetPreallocation(mat, a));
83: } else {
84: Mat_MPIDense *implm = (Mat_MPIDense *)mat->data;
85: PetscCall(MatMPIDenseSetPreallocation(mat, a));
86: impl = (Mat_SeqDense *)implm->A->data;
87: }
88: impl->user_alloc = PETSC_FALSE;
89: }
91: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
92: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
95: #endif