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;
22: PetscErrorCode ierr;
25: PetscViewerGetFormat(viewer, &format);
26: switch (format) {
27: case PETSC_VIEWER_HDF5_PETSC:
28: case PETSC_VIEWER_DEFAULT:
29: case PETSC_VIEWER_NATIVE:
30: case PETSC_VIEWER_HDF5_MAT:
31: break;
32: default:
33: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"PetscViewerFormat %s not supported for HDF5 input.",PetscViewerFormats[format]);
34: }
35: hdf5 = (PetscViewer_HDF5*) viewer->data;
36: /* we store dense matrix columns as blocks, like MATLAB save(filename,variables,'-v7.3') does */
37: hdf5->horizontal = PETSC_TRUE;
39: if (!((PetscObject)mat)->name) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Mat name must be set with PetscObjectSetName() before MatLoad()");
40: #if defined(PETSC_USE_REAL_SINGLE)
41: scalartype = H5T_NATIVE_FLOAT;
42: #elif defined(PETSC_USE_REAL___FLOAT128)
43: #error "HDF5 output with 128 bit floats not supported."
44: #elif defined(PETSC_USE_REAL___FP16)
45: #error "HDF5 output with 16 bit floats not supported."
46: #else
47: scalartype = H5T_NATIVE_DOUBLE;
48: #endif
50: PetscObjectGetComm((PetscObject)mat,&comm);
51: MPI_Comm_rank(comm,&rank);
52: MPI_Comm_size(comm,&size);
53: PetscObjectGetName((PetscObject)mat,&mat_name);
55: /* Convert user-defined rmap and cmap to the dataset layout */
56: PetscLayoutCreate(PetscObjectComm((PetscObject)mat),&vmap);
57: if (mat->rmap->n >= 0 && mat->cmap->N < 0) {
58: /* We need to know mat->cmap->N if user specifies custom mat->rmap->n, otherwise the latter would get ignored below */
59: PetscViewerHDF5ReadSizes(viewer, mat_name, &mat->cmap->N, NULL);
60: }
61: vmap->bs = mat->cmap->N;
62: vmap->n = (mat->rmap->n < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->n * mat->cmap->N;
63: vmap->N = (mat->rmap->N < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->N * mat->cmap->N;
65: /* Read the dataset and setup its layout */
66: /* Note: PetscViewerHDF5ReadSizes_Private takes into account that the dataset is transposed for MATLAB MAT files */
67: PetscViewerHDF5Load(viewer, mat_name, vmap, scalartype, (void**)&a);
69: /* Convert the dataset layout back to rmap and cmap */
70: mat->cmap->N = vmap->bs;
71: mat->rmap->n = vmap->n / mat->cmap->N;
72: mat->rmap->N = vmap->N / mat->cmap->N;
73: PetscLayoutSetUp(mat->rmap);
74: PetscLayoutSetUp(mat->cmap);
75: PetscLayoutDestroy(&vmap);
77: /* TODO adding PetscCopyMode flag to MatSeqDenseSetPreallocation would make this code cleaner and simpler */
78: {
79: PetscBool flg;
80: Mat_SeqDense *impl;
81: PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&flg);
82: if (flg) {
83: impl = (Mat_SeqDense*)mat->data;
84: MatSeqDenseSetPreallocation(mat,a);
85: } else {
86: Mat_MPIDense *implm = (Mat_MPIDense*)mat->data;
87: MatMPIDenseSetPreallocation(mat,a);
88: impl = (Mat_SeqDense*)implm->A->data;
89: }
90: impl->user_alloc = PETSC_FALSE;
91: }
93: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
95: return(0);
96: }
97: #endif