Actual source code: aijhdf5.c
petsc-3.11.4 2019-09-28
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: const char *a_name = NULL, *i_name = NULL, *j_name = NULL, *mat_name = NULL, *c_name = NULL;
16: PetscInt p, m, M, N;
17: PetscInt bs = mat->rmap->bs;
18: PetscBool flg;
19: IS is_i = NULL, is_j = NULL;
20: Vec vec_a = NULL;
21: PetscLayout jmap = NULL;
22: MPI_Comm comm;
23: PetscMPIInt rank, size;
24: PetscErrorCode ierr;
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: SETERRQ1(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: PetscViewerHDF5GetAIJNames(viewer,&i_name,&j_name,&a_name,&c_name);
44: PetscOptionsBegin(comm,NULL,"Options for loading matrix from HDF5","Mat");
45: PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,&flg);
46: PetscOptionsEnd();
47: if (flg) {
48: MatSetBlockSize(mat, bs);
49: }
51: PetscViewerHDF5PushGroup(viewer,mat_name);
52: PetscViewerHDF5ReadAttribute(viewer,NULL,c_name,PETSC_INT,&N);
53: PetscViewerHDF5ReadSizes(viewer, i_name, NULL, &M);
54: --M; /* i has size M+1 as there is global number of nonzeros stored at the end */
56: if (format==PETSC_VIEWER_HDF5_MAT && !mat->symmetric) {
57: /* Swap row and columns layout for unallocated matrix. I want to avoid calling MatTranspose() just to transpose sparsity pattern and layout. */
58: if (!mat->preallocated) {
59: PetscLayout tmp;
60: tmp = mat->rmap; mat->rmap = mat->cmap; mat->cmap = tmp;
61: } else SETERRQ(comm,PETSC_ERR_SUP,"Not for preallocated matrix - we would need to transpose it here which we want to avoid");
62: }
64: /* If global sizes are set, check if they are consistent with that given in the file */
65: if (mat->rmap->N >= 0 && mat->rmap->N != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent # of rows: Matrix in file has (%D) and input matrix has (%D)",mat->rmap->N,M);
66: if (mat->cmap->N >= 0 && mat->cmap->N != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent # of cols: Matrix in file has (%D) and input matrix has (%D)",mat->cmap->N,N);
68: /* Determine ownership of all (block) rows and columns */
69: mat->rmap->N = M;
70: mat->cmap->N = N;
71: PetscLayoutSetUp(mat->rmap);
72: PetscLayoutSetUp(mat->cmap);
73: m = mat->rmap->n;
75: /* Read array i (array of row indices) */
76: PetscMalloc1(m+1, &i); /* allocate i with one more position for local number of nonzeros on each rank */
77: i[0] = i[m] = 0; /* make the last entry always defined - the code block below overwrites it just on last rank */
78: 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 */
79: M++;
80: ISCreate(comm,&is_i);
81: PetscObjectSetName((PetscObject)is_i,i_name);
82: PetscLayoutSetLocalSize(is_i->map,m);
83: PetscLayoutSetSize(is_i->map,M);
84: ISLoad(is_i,viewer);
85: ISGetIndices(is_i,&i_glob);
86: PetscMemcpy(i,i_glob,m*sizeof(PetscInt));
88: /* Reset m and M to the matrix sizes */
89: m = mat->rmap->n;
90: M--;
92: /* Create PetscLayout for j and a vectors; construct ranges first */
93: PetscLayoutCreate(comm,&jmap);
94: PetscCalloc1(size+1, &jmap->range);
95: MPI_Allgather(&i[0], 1, MPIU_INT, jmap->range, 1, MPIU_INT, comm);
96: /* Last rank has global number of nonzeros (= length of j and a arrays) in i[m] (last i entry) so broadcast it */
97: jmap->range[size] = i[m];
98: MPI_Bcast(&jmap->range[size], 1, MPIU_INT, size-1, comm);
99: for (p=size-1; p>0; p--) {
100: if (!jmap->range[p]) jmap->range[p] = jmap->range[p+1]; /* for ranks with 0 rows, take the value from the next processor */
101: }
102: i[m] = jmap->range[rank+1]; /* i[m] (last i entry) is equal to next rank's offset */
103: /* Deduce rstart, rend, n and N from the ranges */
104: PetscLayoutSetUp(jmap);
106: /* Convert global to local indexing of rows */
107: for (p=1; p<m+1; ++p) i[p] -= i[0];
108: i[0] = 0;
110: /* Read array j (array of column indices) */
111: ISCreate(comm,&is_j);
112: PetscObjectSetName((PetscObject)is_j,j_name);
113: PetscLayoutDuplicate(jmap,&is_j->map);
114: ISLoad(is_j,viewer);
115: ISGetIndices(is_j,&j);
117: /* Read array a (array of values) */
118: VecCreate(comm,&vec_a);
119: PetscObjectSetName((PetscObject)vec_a,a_name);
120: PetscLayoutDuplicate(jmap,&vec_a->map);
121: VecLoad(vec_a,viewer);
122: VecGetArrayRead(vec_a,&a);
124: /* populate matrix */
125: if (!((PetscObject)mat)->type_name) {
126: MatSetType(mat,MATAIJ);
127: }
128: MatSeqAIJSetPreallocationCSR(mat,i,j,a);
129: MatMPIAIJSetPreallocationCSR(mat,i,j,a);
130: /*
131: MatSeqBAIJSetPreallocationCSR(mat,bs,i,j,a);
132: MatMPIBAIJSetPreallocationCSR(mat,bs,i,j,a);
133: */
135: if (format==PETSC_VIEWER_HDF5_MAT && !mat->symmetric) {
136: /* Transpose the input matrix back */
137: MatTranspose(mat,MAT_INPLACE_MATRIX,&mat);
138: }
140: PetscViewerHDF5PopGroup(viewer);
141: PetscLayoutDestroy(&jmap);
142: PetscFree(i);
143: ISRestoreIndices(is_i,&i_glob);
144: ISRestoreIndices(is_j,&j);
145: VecRestoreArrayRead(vec_a,&a);
146: ISDestroy(&is_i);
147: ISDestroy(&is_j);
148: VecDestroy(&vec_a);
149: return(0);
150: }
151: #endif