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