Actual source code: aijhdf5.c


  2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <petsc/private/isimpl.h>
  4: #include <petsc/private/vecimpl.h>
  5: #include <petsclayouthdf5.h>

  7: #if defined(PETSC_HAVE_HDF5)
  8: PetscErrorCode MatLoad_AIJ_HDF5(Mat mat, PetscViewer viewer)
  9: {
 10:   PetscViewerFormat  format;
 11:   const PetscInt    *i_glob = NULL;
 12:   PetscInt          *i      = NULL;
 13:   const PetscInt    *j      = NULL;
 14:   const PetscScalar *a      = NULL;
 15:   char              *a_name = NULL, *i_name = NULL, *j_name = NULL, *c_name = NULL;
 16:   const char        *mat_name = NULL;
 17:   PetscInt           p, m, M, N;
 18:   PetscInt           bs = mat->rmap->bs;
 19:   PetscInt          *range;
 20:   PetscBool          flg;
 21:   IS                 is_i = NULL, is_j = NULL;
 22:   Vec                vec_a = NULL;
 23:   PetscLayout        jmap  = NULL;
 24:   MPI_Comm           comm;
 25:   PetscMPIInt        rank, size;

 27:   PetscViewerGetFormat(viewer, &format);
 28:   switch (format) {
 29:   case PETSC_VIEWER_HDF5_PETSC:
 30:   case PETSC_VIEWER_DEFAULT:
 31:   case PETSC_VIEWER_NATIVE:
 32:   case PETSC_VIEWER_HDF5_MAT:
 33:     break;
 34:   default:
 35:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 input.", PetscViewerFormats[format]);
 36:   }

 38:   PetscObjectGetComm((PetscObject)mat, &comm);
 39:   MPI_Comm_rank(comm, &rank);
 40:   MPI_Comm_size(comm, &size);
 41:   PetscObjectGetName((PetscObject)mat, &mat_name);
 42:   if (format == PETSC_VIEWER_HDF5_MAT) {
 43:     PetscStrallocpy("jc", &i_name);
 44:     PetscStrallocpy("ir", &j_name);
 45:     PetscStrallocpy("data", &a_name);
 46:     PetscStrallocpy("MATLAB_sparse", &c_name);
 47:   } else {
 48:     /* TODO Once corresponding MatView is implemented, change the names to i,j,a */
 49:     /* TODO Maybe there could be both namings in the file, using "symbolic link" features of HDF5. */
 50:     PetscStrallocpy("jc", &i_name);
 51:     PetscStrallocpy("ir", &j_name);
 52:     PetscStrallocpy("data", &a_name);
 53:     PetscStrallocpy("MATLAB_sparse", &c_name);
 54:   }

 56:   PetscOptionsBegin(comm, NULL, "Options for loading matrix from HDF5", "Mat");
 57:   PetscOptionsInt("-matload_block_size", "Set the blocksize used to store the matrix", "MatLoad", bs, &bs, &flg);
 58:   PetscOptionsEnd();
 59:   if (flg) MatSetBlockSize(mat, bs);

 61:   PetscViewerHDF5PushGroup(viewer, mat_name);
 62:   PetscViewerHDF5ReadAttribute(viewer, NULL, c_name, PETSC_INT, NULL, &N);
 63:   PetscViewerHDF5ReadSizes(viewer, i_name, NULL, &M);
 64:   --M; /* i has size M+1 as there is global number of nonzeros stored at the end */

 66:   if (format == PETSC_VIEWER_HDF5_MAT && mat->symmetric != PETSC_BOOL3_TRUE) {
 67:     /* Swap row and columns layout for unallocated matrix. I want to avoid calling MatTranspose() just to transpose sparsity pattern and layout. */
 68:     if (!mat->preallocated) {
 69:       PetscLayout tmp;
 70:       tmp       = mat->rmap;
 71:       mat->rmap = mat->cmap;
 72:       mat->cmap = tmp;
 73:     } else SETERRQ(comm, PETSC_ERR_SUP, "Not for preallocated matrix - we would need to transpose it here which we want to avoid");
 74:   }

 76:   /* If global sizes are set, check if they are consistent with that given in the file */

 80:   /* Determine ownership of all (block) rows and columns */
 81:   mat->rmap->N = M;
 82:   mat->cmap->N = N;
 83:   PetscLayoutSetUp(mat->rmap);
 84:   PetscLayoutSetUp(mat->cmap);
 85:   m = mat->rmap->n;

 87:   /* Read array i (array of row indices) */
 88:   PetscMalloc1(m + 1, &i); /* allocate i with one more position for local number of nonzeros on each rank */
 89:   i[0] = i[m] = 0;                    /* make the last entry always defined - the code block below overwrites it just on last rank */
 90:   if (rank == size - 1) m++;          /* in the loaded array i_glob, only the last rank has one more position with the global number of nonzeros */
 91:   M++;
 92:   ISCreate(comm, &is_i);
 93:   PetscObjectSetName((PetscObject)is_i, i_name);
 94:   PetscLayoutSetLocalSize(is_i->map, m);
 95:   PetscLayoutSetSize(is_i->map, M);
 96:   ISLoad(is_i, viewer);
 97:   ISGetIndices(is_i, &i_glob);
 98:   PetscArraycpy(i, i_glob, m);

100:   /* Reset m and M to the matrix sizes */
101:   m = mat->rmap->n;
102:   M--;

104:   /* Create PetscLayout for j and a vectors; construct ranges first */
105:   PetscMalloc1(size + 1, &range);
106:   MPI_Allgather(i, 1, MPIU_INT, range, 1, MPIU_INT, comm);
107:   /* Last rank has global number of nonzeros (= length of j and a arrays) in i[m] (last i entry) so broadcast it */
108:   range[size] = i[m];
109:   MPI_Bcast(&range[size], 1, MPIU_INT, size - 1, comm);
110:   for (p = size - 1; p > 0; p--) {
111:     if (!range[p]) range[p] = range[p + 1]; /* for ranks with 0 rows, take the value from the next processor */
112:   }
113:   i[m] = range[rank + 1]; /* i[m] (last i entry) is equal to next rank's offset */
114:   /* Deduce rstart, rend, n and N from the ranges */
115:   PetscLayoutCreateFromRanges(comm, range, PETSC_OWN_POINTER, 1, &jmap);

117:   /* Convert global to local indexing of rows */
118:   for (p = 1; p < m + 1; ++p) i[p] -= i[0];
119:   i[0] = 0;

121:   /* Read array j (array of column indices) */
122:   ISCreate(comm, &is_j);
123:   PetscObjectSetName((PetscObject)is_j, j_name);
124:   PetscLayoutDuplicate(jmap, &is_j->map);
125:   ISLoad(is_j, viewer);
126:   ISGetIndices(is_j, &j);

128:   /* Read array a (array of values) */
129:   VecCreate(comm, &vec_a);
130:   PetscObjectSetName((PetscObject)vec_a, a_name);
131:   PetscLayoutDuplicate(jmap, &vec_a->map);
132:   VecLoad(vec_a, viewer);
133:   VecGetArrayRead(vec_a, &a);

135:   /* populate matrix */
136:   if (!((PetscObject)mat)->type_name) MatSetType(mat, MATAIJ);
137:   MatSeqAIJSetPreallocationCSR(mat, i, j, a);
138:   MatMPIAIJSetPreallocationCSR(mat, i, j, a);
139:   /*
140:   MatSeqBAIJSetPreallocationCSR(mat,bs,i,j,a);
141:   MatMPIBAIJSetPreallocationCSR(mat,bs,i,j,a);
142:   */

144:   if (format == PETSC_VIEWER_HDF5_MAT && mat->symmetric != PETSC_BOOL3_TRUE) {
145:     /* Transpose the input matrix back */
146:     MatTranspose(mat, MAT_INPLACE_MATRIX, &mat);
147:   }

149:   PetscViewerHDF5PopGroup(viewer);
150:   PetscFree(i_name);
151:   PetscFree(j_name);
152:   PetscFree(a_name);
153:   PetscFree(c_name);
154:   PetscLayoutDestroy(&jmap);
155:   PetscFree(i);
156:   ISRestoreIndices(is_i, &i_glob);
157:   ISRestoreIndices(is_j, &j);
158:   VecRestoreArrayRead(vec_a, &a);
159:   ISDestroy(&is_i);
160:   ISDestroy(&is_j);
161:   VecDestroy(&vec_a);
162:   return 0;
163: }
164: #endif