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;
26: PetscErrorCode ierr;
28: PetscViewerGetFormat(viewer, &format);
29: switch (format) {
30: case PETSC_VIEWER_HDF5_PETSC:
31: case PETSC_VIEWER_DEFAULT:
32: case PETSC_VIEWER_NATIVE:
33: case PETSC_VIEWER_HDF5_MAT:
34: break;
35: default:
36: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"PetscViewerFormat %s not supported for HDF5 input.",PetscViewerFormats[format]);
37: }
39: PetscObjectGetComm((PetscObject)mat,&comm);
40: MPI_Comm_rank(comm,&rank);
41: MPI_Comm_size(comm,&size);
42: PetscObjectGetName((PetscObject)mat,&mat_name);
43: if (format==PETSC_VIEWER_HDF5_MAT) {
44: PetscStrallocpy("jc",&i_name);
45: PetscStrallocpy("ir",&j_name);
46: PetscStrallocpy("data",&a_name);
47: PetscStrallocpy("MATLAB_sparse",&c_name);
48: } else {
49: /* TODO Once corresponding MatView is implemented, change the names to i,j,a */
50: /* TODO Maybe there could be both namings in the file, using "symbolic link" features of HDF5. */
51: PetscStrallocpy("jc",&i_name);
52: PetscStrallocpy("ir",&j_name);
53: PetscStrallocpy("data",&a_name);
54: PetscStrallocpy("MATLAB_sparse",&c_name);
55: }
57: PetscOptionsBegin(comm,NULL,"Options for loading matrix from HDF5","Mat");
58: PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,&flg);
59: PetscOptionsEnd();
60: if (flg) {
61: MatSetBlockSize(mat, bs);
62: }
64: PetscViewerHDF5PushGroup(viewer,mat_name);
65: PetscViewerHDF5ReadAttribute(viewer,NULL,c_name,PETSC_INT,NULL,&N);
66: PetscViewerHDF5ReadSizes(viewer, i_name, NULL, &M);
67: --M; /* i has size M+1 as there is global number of nonzeros stored at the end */
69: if (format==PETSC_VIEWER_HDF5_MAT && !mat->symmetric) {
70: /* Swap row and columns layout for unallocated matrix. I want to avoid calling MatTranspose() just to transpose sparsity pattern and layout. */
71: if (!mat->preallocated) {
72: PetscLayout tmp;
73: tmp = mat->rmap; mat->rmap = mat->cmap; mat->cmap = tmp;
74: } else SETERRQ(comm,PETSC_ERR_SUP,"Not for preallocated matrix - we would need to transpose it here which we want to avoid");
75: }
77: /* If global sizes are set, check if they are consistent with that given in the file */
81: /* Determine ownership of all (block) rows and columns */
82: mat->rmap->N = M;
83: mat->cmap->N = N;
84: PetscLayoutSetUp(mat->rmap);
85: PetscLayoutSetUp(mat->cmap);
86: m = mat->rmap->n;
88: /* Read array i (array of row indices) */
89: PetscMalloc1(m+1, &i); /* allocate i with one more position for local number of nonzeros on each rank */
90: i[0] = i[m] = 0; /* make the last entry always defined - the code block below overwrites it just on last rank */
91: 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 */
92: M++;
93: ISCreate(comm,&is_i);
94: PetscObjectSetName((PetscObject)is_i,i_name);
95: PetscLayoutSetLocalSize(is_i->map,m);
96: PetscLayoutSetSize(is_i->map,M);
97: ISLoad(is_i,viewer);
98: ISGetIndices(is_i,&i_glob);
99: PetscArraycpy(i,i_glob,m);
101: /* Reset m and M to the matrix sizes */
102: m = mat->rmap->n;
103: M--;
105: /* Create PetscLayout for j and a vectors; construct ranges first */
106: PetscMalloc1(size+1, &range);
107: MPI_Allgather(i, 1, MPIU_INT, range, 1, MPIU_INT, comm);
108: /* Last rank has global number of nonzeros (= length of j and a arrays) in i[m] (last i entry) so broadcast it */
109: range[size] = i[m];
110: MPI_Bcast(&range[size], 1, MPIU_INT, size-1, comm);
111: for (p=size-1; p>0; p--) {
112: if (!range[p]) range[p] = range[p+1]; /* for ranks with 0 rows, take the value from the next processor */
113: }
114: i[m] = range[rank+1]; /* i[m] (last i entry) is equal to next rank's offset */
115: /* Deduce rstart, rend, n and N from the ranges */
116: PetscLayoutCreateFromRanges(comm,range,PETSC_OWN_POINTER,1,&jmap);
118: /* Convert global to local indexing of rows */
119: for (p=1; p<m+1; ++p) i[p] -= i[0];
120: i[0] = 0;
122: /* Read array j (array of column indices) */
123: ISCreate(comm,&is_j);
124: PetscObjectSetName((PetscObject)is_j,j_name);
125: PetscLayoutDuplicate(jmap,&is_j->map);
126: ISLoad(is_j,viewer);
127: ISGetIndices(is_j,&j);
129: /* Read array a (array of values) */
130: VecCreate(comm,&vec_a);
131: PetscObjectSetName((PetscObject)vec_a,a_name);
132: PetscLayoutDuplicate(jmap,&vec_a->map);
133: VecLoad(vec_a,viewer);
134: VecGetArrayRead(vec_a,&a);
136: /* populate matrix */
137: if (!((PetscObject)mat)->type_name) {
138: MatSetType(mat,MATAIJ);
139: }
140: MatSeqAIJSetPreallocationCSR(mat,i,j,a);
141: MatMPIAIJSetPreallocationCSR(mat,i,j,a);
142: /*
143: MatSeqBAIJSetPreallocationCSR(mat,bs,i,j,a);
144: MatMPIBAIJSetPreallocationCSR(mat,bs,i,j,a);
145: */
147: if (format==PETSC_VIEWER_HDF5_MAT && !mat->symmetric) {
148: /* Transpose the input matrix back */
149: MatTranspose(mat,MAT_INPLACE_MATRIX,&mat);
150: }
152: PetscViewerHDF5PopGroup(viewer);
153: PetscFree(i_name);
154: PetscFree(j_name);
155: PetscFree(a_name);
156: PetscFree(c_name);
157: PetscLayoutDestroy(&jmap);
158: PetscFree(i);
159: ISRestoreIndices(is_i,&i_glob);
160: ISRestoreIndices(is_j,&j);
161: VecRestoreArrayRead(vec_a,&a);
162: ISDestroy(&is_i);
163: ISDestroy(&is_j);
164: VecDestroy(&vec_a);
165: return 0;
166: }
167: #endif