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 <petsclayouthdf5.h>

 10: #if defined(PETSC_HAVE_HDF5)
 11: PetscErrorCode MatLoad_Dense_HDF5(Mat mat, PetscViewer viewer)
 12: {
 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:   PetscErrorCode      ierr;

 23:   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:       SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"PetscViewerFormat %s not supported for HDF5 input.",PetscViewerFormats[format]);
 32:   }

 34:   if (!((PetscObject)mat)->name) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Mat name must be set with PetscObjectSetName() before MatLoad()");
 35: #if defined(PETSC_USE_REAL_SINGLE)
 36:   scalartype = H5T_NATIVE_FLOAT;
 37: #elif defined(PETSC_USE_REAL___FLOAT128)
 38: #error "HDF5 output with 128 bit floats not supported."
 39: #elif defined(PETSC_USE_REAL___FP16)
 40: #error "HDF5 output with 16 bit floats not supported."
 41: #else
 42:   scalartype = H5T_NATIVE_DOUBLE;
 43: #endif

 45:   PetscObjectGetComm((PetscObject)mat,&comm);
 46:   MPI_Comm_rank(comm,&rank);
 47:   MPI_Comm_size(comm,&size);
 48:   PetscObjectGetName((PetscObject)mat,&mat_name);

 50:   /* Convert user-defined rmap and cmap to the dataset layout */
 51:   PetscLayoutCreate(PetscObjectComm((PetscObject)mat),&vmap);
 52:   if (mat->rmap->n >= 0 && mat->cmap->N < 0) {
 53:     /* We need to know mat->cmap->N if user specifies custom mat->rmap->n, otherwise the latter would get ignored below */
 54:     PetscViewerHDF5ReadSizes(viewer, mat_name, &mat->cmap->N, NULL);
 55:   }
 56:   vmap->bs = mat->cmap->N;
 57:   vmap->n = (mat->rmap->n < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->n * mat->cmap->N;
 58:   vmap->N = (mat->rmap->N < 0 || mat->cmap->N < 0) ? -1 : mat->rmap->N * mat->cmap->N;

 60:   /* Read the dataset and setup its layout */
 61:   /* Note: PetscViewerHDF5ReadSizes_Private takes into account that the dataset is transposed for MATLAB MAT files */
 62:   PetscViewerHDF5Load(viewer, mat_name, vmap, scalartype, (void**)&a);

 64:   /* Convert the dataset layout back to rmap and cmap */
 65:   mat->cmap->N = vmap->bs;
 66:   mat->rmap->n = vmap->n / mat->cmap->N;
 67:   mat->rmap->N = vmap->N / mat->cmap->N;
 68:   PetscLayoutSetUp(mat->rmap);
 69:   PetscLayoutSetUp(mat->cmap);
 70:   PetscLayoutDestroy(&vmap);

 72:   /* TODO adding PetscCopyMode flag to MatSeqDenseSetPreallocation would make this code cleaner and simpler */
 73:   {
 74:     PetscBool flg;
 75:     Mat_SeqDense *impl;
 76:     PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&flg);
 77:     if (flg) {
 78:       impl = (Mat_SeqDense*)mat->data;
 79:       MatSeqDenseSetPreallocation(mat,a);
 80:     } else {
 81:       Mat_MPIDense *implm = (Mat_MPIDense*)mat->data;
 82:       MatMPIDenseSetPreallocation(mat,a);
 83:       impl = (Mat_SeqDense*)implm->A->data;
 84:     }
 85:     impl->user_alloc = PETSC_FALSE;
 86:   }

 88:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 89:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
 90:   return(0);
 91: }
 92: #endif