1*2e1d0745SJose E. Roman /* TODO change to
2*2e1d0745SJose E. Roman #include <../src/mat/impls/dense/seq/dense.h>
3*2e1d0745SJose E. Roman */
4*2e1d0745SJose E. Roman #include <../src/mat/impls/dense/mpi/mpidense.h>
5*2e1d0745SJose E. Roman #include <petsc/private/isimpl.h>
6*2e1d0745SJose E. Roman #include <petsc/private/vecimpl.h>
7*2e1d0745SJose E. Roman #include <petsc/private/viewerhdf5impl.h>
8*2e1d0745SJose E. Roman #include <petsclayouthdf5.h>
9*2e1d0745SJose E. Roman
MatLoad_Dense_HDF5(Mat mat,PetscViewer viewer)10*2e1d0745SJose E. Roman PetscErrorCode MatLoad_Dense_HDF5(Mat mat, PetscViewer viewer)
11*2e1d0745SJose E. Roman {
12*2e1d0745SJose E. Roman PetscViewer_HDF5 *hdf5;
13*2e1d0745SJose E. Roman hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
14*2e1d0745SJose E. Roman PetscLayout vmap;
15*2e1d0745SJose E. Roman PetscViewerFormat format;
16*2e1d0745SJose E. Roman PetscScalar *a = NULL;
17*2e1d0745SJose E. Roman const char *mat_name = NULL;
18*2e1d0745SJose E. Roman MPI_Comm comm;
19*2e1d0745SJose E. Roman PetscMPIInt rank, size;
20*2e1d0745SJose E. Roman
21*2e1d0745SJose E. Roman PetscFunctionBegin;
22*2e1d0745SJose E. Roman PetscCall(PetscViewerGetFormat(viewer, &format));
23*2e1d0745SJose E. Roman switch (format) {
24*2e1d0745SJose E. Roman case PETSC_VIEWER_HDF5_PETSC:
25*2e1d0745SJose E. Roman case PETSC_VIEWER_DEFAULT:
26*2e1d0745SJose E. Roman case PETSC_VIEWER_NATIVE:
27*2e1d0745SJose E. Roman case PETSC_VIEWER_HDF5_MAT:
28*2e1d0745SJose E. Roman break;
29*2e1d0745SJose E. Roman default:
30*2e1d0745SJose E. Roman SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 input.", PetscViewerFormats[format]);
31*2e1d0745SJose E. Roman }
32*2e1d0745SJose E. Roman hdf5 = (PetscViewer_HDF5 *)viewer->data;
33*2e1d0745SJose E. Roman /* we store dense matrix columns as blocks, like MATLAB save(filename,variables,'-v7.3') does */
34*2e1d0745SJose E. Roman hdf5->horizontal = PETSC_TRUE;
35*2e1d0745SJose E. Roman
36*2e1d0745SJose E. Roman PetscCheck(((PetscObject)mat)->name, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Mat name must be set with PetscObjectSetName() before MatLoad()");
37*2e1d0745SJose E. Roman #if defined(PETSC_USE_REAL_SINGLE)
38*2e1d0745SJose E. Roman scalartype = H5T_NATIVE_FLOAT;
39*2e1d0745SJose E. Roman #elif defined(PETSC_USE_REAL___FLOAT128)
40*2e1d0745SJose E. Roman #error "HDF5 output with 128 bit floats not supported."
41*2e1d0745SJose E. Roman #elif defined(PETSC_USE_REAL___FP16)
42*2e1d0745SJose E. Roman #error "HDF5 output with 16 bit floats not supported."
43*2e1d0745SJose E. Roman #else
44*2e1d0745SJose E. Roman scalartype = H5T_NATIVE_DOUBLE;
45*2e1d0745SJose E. Roman #endif
46*2e1d0745SJose E. Roman
47*2e1d0745SJose E. Roman PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
48*2e1d0745SJose E. Roman PetscCallMPI(MPI_Comm_rank(comm, &rank));
49*2e1d0745SJose E. Roman PetscCallMPI(MPI_Comm_size(comm, &size));
50*2e1d0745SJose E. Roman PetscCall(PetscObjectGetName((PetscObject)mat, &mat_name));
51*2e1d0745SJose E. Roman
52*2e1d0745SJose E. Roman /* Convert user-defined rmap and cmap to the dataset layout */
53*2e1d0745SJose E. Roman PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)mat), &vmap));
54*2e1d0745SJose E. Roman if (mat->rmap->n >= 0 && mat->cmap->N < 0) {
55*2e1d0745SJose E. Roman /* We need to know mat->cmap->N if user specifies custom mat->rmap->n, otherwise the latter would get ignored below */
56*2e1d0745SJose E. Roman PetscCall(PetscViewerHDF5ReadSizes(viewer, mat_name, &mat->cmap->N, NULL));
57*2e1d0745SJose E. Roman }
58*2e1d0745SJose E. Roman vmap->bs = mat->cmap->N;
59*2e1d0745SJose E. Roman vmap->n = (mat->rmap->n < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->n * mat->cmap->N;
60*2e1d0745SJose E. Roman vmap->N = (mat->rmap->N < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->N * mat->cmap->N;
61*2e1d0745SJose E. Roman
62*2e1d0745SJose E. Roman /* Read the dataset and setup its layout */
63*2e1d0745SJose E. Roman /* Note: PetscViewerHDF5ReadSizes_Private takes into account that the dataset is transposed for MATLAB MAT files */
64*2e1d0745SJose E. Roman PetscCall(PetscViewerHDF5Load(viewer, mat_name, vmap, scalartype, (void **)&a));
65*2e1d0745SJose E. Roman
66*2e1d0745SJose E. Roman /* Convert the dataset layout back to rmap and cmap */
67*2e1d0745SJose E. Roman mat->cmap->N = vmap->bs;
68*2e1d0745SJose E. Roman mat->rmap->n = vmap->n / mat->cmap->N;
69*2e1d0745SJose E. Roman mat->rmap->N = vmap->N / mat->cmap->N;
70*2e1d0745SJose E. Roman PetscCall(PetscLayoutSetUp(mat->rmap));
71*2e1d0745SJose E. Roman PetscCall(PetscLayoutSetUp(mat->cmap));
72*2e1d0745SJose E. Roman PetscCall(PetscLayoutDestroy(&vmap));
73*2e1d0745SJose E. Roman
74*2e1d0745SJose E. Roman /* TODO adding PetscCopyMode flag to MatSeqDenseSetPreallocation would make this code cleaner and simpler */
75*2e1d0745SJose E. Roman {
76*2e1d0745SJose E. Roman PetscBool flg;
77*2e1d0745SJose E. Roman Mat_SeqDense *impl;
78*2e1d0745SJose E. Roman PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
79*2e1d0745SJose E. Roman if (flg) {
80*2e1d0745SJose E. Roman impl = (Mat_SeqDense *)mat->data;
81*2e1d0745SJose E. Roman PetscCall(MatSeqDenseSetPreallocation(mat, a));
82*2e1d0745SJose E. Roman } else {
83*2e1d0745SJose E. Roman Mat_MPIDense *implm = (Mat_MPIDense *)mat->data;
84*2e1d0745SJose E. Roman PetscCall(MatMPIDenseSetPreallocation(mat, a));
85*2e1d0745SJose E. Roman impl = (Mat_SeqDense *)implm->A->data;
86*2e1d0745SJose E. Roman }
87*2e1d0745SJose E. Roman impl->user_alloc = PETSC_FALSE;
88*2e1d0745SJose E. Roman }
89*2e1d0745SJose E. Roman
90*2e1d0745SJose E. Roman PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
91*2e1d0745SJose E. Roman PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
92*2e1d0745SJose E. Roman PetscFunctionReturn(PETSC_SUCCESS);
93*2e1d0745SJose E. Roman }
94