1 /*
2 Creates hypre ijmatrix from PETSc matrix
3 */
4 #include <petscsys.h>
5 #include <petsc/private/petschypre.h>
6 #include <petsc/private/matimpl.h>
7 #include <petscdmda.h> /*I "petscdmda.h" I*/
8 #include <../src/dm/impls/da/hypre/mhyp.h>
9
10 /*MC
11 MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
12 based on the hypre HYPRE_StructMatrix.
13
14 Level: intermediate
15
16 Notes:
17 Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
18 be defined by a `DMDA`.
19
20 The matrix needs a `DMDA` associated with it by either a call to `MatSetDM()` or if the matrix is obtained from `DMCreateMatrix()`
21
22 .seealso: `MatCreate()`, `PCPFMG`, `MatSetDM()`, `DMCreateMatrix()`
23 M*/
24
MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)25 static PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
26 {
27 HYPRE_Int index[3], entries[9];
28 PetscInt i, j, stencil, row;
29 HYPRE_Complex *values = (HYPRE_Complex *)y;
30 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
31
32 PetscFunctionBegin;
33 PetscCheck(ncol <= 9, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncol %" PetscInt_FMT " > 9 too large", ncol);
34 for (i = 0; i < nrow; i++) {
35 for (j = 0; j < ncol; j++) {
36 stencil = icol[j] - irow[i];
37 if (!stencil) {
38 entries[j] = 3;
39 } else if (stencil == -1) {
40 entries[j] = 2;
41 } else if (stencil == 1) {
42 entries[j] = 4;
43 } else if (stencil == -ex->gnx) {
44 entries[j] = 1;
45 } else if (stencil == ex->gnx) {
46 entries[j] = 5;
47 } else if (stencil == -ex->gnxgny) {
48 entries[j] = 0;
49 } else if (stencil == ex->gnxgny) {
50 entries[j] = 6;
51 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
52 }
53 row = ex->gindices[irow[i]] - ex->rstart;
54 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
55 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
56 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
57 if (addv == ADD_VALUES) PetscCallHYPRE(HYPRE_StructMatrixAddToValues(ex->hmat, index, (HYPRE_Int)ncol, entries, values));
58 else PetscCallHYPRE(HYPRE_StructMatrixSetValues(ex->hmat, index, (HYPRE_Int)ncol, entries, values));
59 values += ncol;
60 }
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)64 static PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
65 {
66 HYPRE_Int index[3], entries[7] = {0, 1, 2, 3, 4, 5, 6};
67 PetscInt row, i;
68 HYPRE_Complex values[7];
69 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
70
71 PetscFunctionBegin;
72 PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support");
73 PetscCall(PetscArrayzero(values, 7));
74 PetscCall(PetscHYPREScalarCast(d, &values[3]));
75 for (i = 0; i < nrow; i++) {
76 row = ex->gindices[irow[i]] - ex->rstart;
77 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
78 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
79 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
80 PetscCallHYPRE(HYPRE_StructMatrixSetValues(ex->hmat, index, 7, entries, values));
81 }
82 PetscCallHYPRE(HYPRE_StructMatrixAssemble(ex->hmat));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
MatZeroEntries_HYPREStruct_3d(Mat mat)86 static PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
87 {
88 HYPRE_Int indices[7] = {0, 1, 2, 3, 4, 5, 6};
89 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
90
91 PetscFunctionBegin;
92 /* hypre has no public interface to do this */
93 PetscCallHYPRE(hypre_StructMatrixClearBoxValues(ex->hmat, &ex->hbox, 7, indices, 0, 1));
94 PetscCallHYPRE(HYPRE_StructMatrixAssemble(ex->hmat));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
MatSetUp_HYPREStruct(Mat mat)98 static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
99 {
100 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
101 HYPRE_Int sw[6];
102 HYPRE_Int hlower[3], hupper[3], period[3] = {0, 0, 0};
103 PetscInt dim, dof, psw, Nx, Ny, Nz, nx, ny, nz, ilower[3], iupper[3], ssize, i;
104 DMBoundaryType px, py, pz;
105 DMDAStencilType st;
106 ISLocalToGlobalMapping ltog;
107 DM da;
108
109 PetscFunctionBegin;
110 PetscCall(MatGetDM(mat, (DM *)&da));
111 ex->da = da;
112 PetscCall(PetscObjectReference((PetscObject)da));
113
114 PetscCall(DMDAGetInfo(ex->da, &dim, &Nx, &Ny, &Nz, 0, 0, 0, &dof, &psw, &px, &py, &pz, &st));
115 PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
116
117 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
118 iupper[0] += ilower[0] - 1;
119 iupper[1] += ilower[1] - 1;
120 iupper[2] += ilower[2] - 1;
121 hlower[0] = (HYPRE_Int)ilower[0];
122 hlower[1] = (HYPRE_Int)ilower[1];
123 hlower[2] = (HYPRE_Int)ilower[2];
124 hupper[0] = (HYPRE_Int)iupper[0];
125 hupper[1] = (HYPRE_Int)iupper[1];
126 hupper[2] = (HYPRE_Int)iupper[2];
127
128 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
129 ex->hbox.imin[0] = hlower[0];
130 ex->hbox.imin[1] = hlower[1];
131 ex->hbox.imin[2] = hlower[2];
132 ex->hbox.imax[0] = hupper[0];
133 ex->hbox.imax[1] = hupper[1];
134 ex->hbox.imax[2] = hupper[2];
135
136 /* create the hypre grid object and set its information */
137 PetscCheck(dof <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Currently only support for scalar problems");
138 if (px || py || pz) {
139 if (px == DM_BOUNDARY_PERIODIC) period[0] = (HYPRE_Int)Nx;
140 if (py == DM_BOUNDARY_PERIODIC) period[1] = (HYPRE_Int)Ny;
141 if (pz == DM_BOUNDARY_PERIODIC) period[2] = (HYPRE_Int)Nz;
142 }
143 PetscCallHYPRE(HYPRE_StructGridCreate(ex->hcomm, (HYPRE_Int)dim, &ex->hgrid));
144 PetscCallHYPRE(HYPRE_StructGridSetExtents(ex->hgrid, hlower, hupper));
145 PetscCallHYPRE(HYPRE_StructGridSetPeriodic(ex->hgrid, period));
146 PetscCallHYPRE(HYPRE_StructGridAssemble(ex->hgrid));
147
148 sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = (HYPRE_Int)psw;
149 PetscCallHYPRE(HYPRE_StructGridSetNumGhost(ex->hgrid, sw));
150
151 /* create the hypre stencil object and set its information */
152 PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils");
153 PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils");
154 if (dim == 1) {
155 HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
156 ssize = 3;
157 PetscCallHYPRE(HYPRE_StructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->hstencil));
158 for (i = 0; i < ssize; i++) PetscCallHYPRE(HYPRE_StructStencilSetElement(ex->hstencil, (HYPRE_Int)i, offsets[i]));
159 } else if (dim == 2) {
160 HYPRE_Int offsets[5][2] = {
161 {0, -1},
162 {-1, 0 },
163 {0, 0 },
164 {1, 0 },
165 {0, 1 }
166 };
167 ssize = 5;
168 PetscCallHYPRE(HYPRE_StructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->hstencil));
169 for (i = 0; i < ssize; i++) PetscCallHYPRE(HYPRE_StructStencilSetElement(ex->hstencil, (HYPRE_Int)i, offsets[i]));
170 } else if (dim == 3) {
171 HYPRE_Int offsets[7][3] = {
172 {0, 0, -1},
173 {0, -1, 0 },
174 {-1, 0, 0 },
175 {0, 0, 0 },
176 {1, 0, 0 },
177 {0, 1, 0 },
178 {0, 0, 1 }
179 };
180 ssize = 7;
181 PetscCallHYPRE(HYPRE_StructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->hstencil));
182 for (i = 0; i < ssize; i++) PetscCallHYPRE(HYPRE_StructStencilSetElement(ex->hstencil, (HYPRE_Int)i, offsets[i]));
183 }
184
185 /* create the HYPRE vector for rhs and solution */
186 PetscCallHYPRE(HYPRE_StructVectorCreate(ex->hcomm, ex->hgrid, &ex->hb));
187 PetscCallHYPRE(HYPRE_StructVectorCreate(ex->hcomm, ex->hgrid, &ex->hx));
188 PetscCallHYPRE(HYPRE_StructVectorInitialize(ex->hb));
189 PetscCallHYPRE(HYPRE_StructVectorInitialize(ex->hx));
190 PetscCallHYPRE(HYPRE_StructVectorAssemble(ex->hb));
191 PetscCallHYPRE(HYPRE_StructVectorAssemble(ex->hx));
192
193 /* create the hypre matrix object and set its information */
194 PetscCallHYPRE(HYPRE_StructMatrixCreate(ex->hcomm, ex->hgrid, ex->hstencil, &ex->hmat));
195 PetscCallHYPRE(HYPRE_StructGridDestroy(ex->hgrid));
196 PetscCallHYPRE(HYPRE_StructStencilDestroy(ex->hstencil));
197 if (ex->needsinitialization) {
198 PetscCallHYPRE(HYPRE_StructMatrixInitialize(ex->hmat));
199 ex->needsinitialization = PETSC_FALSE;
200 }
201
202 /* set the global and local sizes of the matrix */
203 PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz));
204 PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE));
205 PetscCall(PetscLayoutSetBlockSize(mat->rmap, 1));
206 PetscCall(PetscLayoutSetBlockSize(mat->cmap, 1));
207 PetscCall(PetscLayoutSetUp(mat->rmap));
208 PetscCall(PetscLayoutSetUp(mat->cmap));
209 mat->preallocated = PETSC_TRUE;
210
211 PetscCheck(dim == 3, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");
212 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
213 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
214 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
215
216 /* PetscCall(MatZeroEntries_HYPREStruct_3d(mat)); */
217
218 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
219 PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL));
220 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og));
221 PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices));
222 PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, 0));
223 ex->gnxgny *= ex->gnx;
224 PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, 0));
225 ex->nxny = ex->nx * ex->ny;
226 PetscFunctionReturn(PETSC_SUCCESS);
227 }
228
MatMult_HYPREStruct(Mat A,Vec x,Vec y)229 static PetscErrorCode MatMult_HYPREStruct(Mat A, Vec x, Vec y)
230 {
231 const PetscScalar *xx;
232 PetscScalar *yy;
233 PetscInt ilower[3], iupper[3];
234 HYPRE_Int hlower[3], hupper[3];
235 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)A->data;
236
237 PetscFunctionBegin;
238 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
239 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
240 iupper[0] += ilower[0] - 1;
241 iupper[1] += ilower[1] - 1;
242 iupper[2] += ilower[2] - 1;
243 hlower[0] = (HYPRE_Int)ilower[0];
244 hlower[1] = (HYPRE_Int)ilower[1];
245 hlower[2] = (HYPRE_Int)ilower[2];
246 hupper[0] = (HYPRE_Int)iupper[0];
247 hupper[1] = (HYPRE_Int)iupper[1];
248 hupper[2] = (HYPRE_Int)iupper[2];
249
250 /* copy x values over to hypre */
251 PetscCallHYPRE(HYPRE_StructVectorSetConstantValues(mx->hb, 0.0));
252 PetscCall(VecGetArrayRead(x, &xx));
253 PetscCallHYPRE(HYPRE_StructVectorSetBoxValues(mx->hb, hlower, hupper, (HYPRE_Complex *)xx));
254 PetscCall(VecRestoreArrayRead(x, &xx));
255 PetscCallHYPRE(HYPRE_StructVectorAssemble(mx->hb));
256 PetscCallHYPRE(HYPRE_StructMatrixMatvec(1.0, mx->hmat, mx->hb, 0.0, mx->hx));
257
258 /* copy solution values back to PETSc */
259 PetscCall(VecGetArray(y, &yy));
260 PetscCallHYPRE(HYPRE_StructVectorGetBoxValues(mx->hx, hlower, hupper, (HYPRE_Complex *)yy));
261 PetscCall(VecRestoreArray(y, &yy));
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)265 static PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode)
266 {
267 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
268
269 PetscFunctionBegin;
270 PetscCallHYPRE(HYPRE_StructMatrixAssemble(ex->hmat));
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
MatZeroEntries_HYPREStruct(Mat mat)274 static PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
275 {
276 PetscFunctionBegin;
277 /* before the DMDA is set to the matrix the zero doesn't need to do anything */
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
MatDestroy_HYPREStruct(Mat mat)281 static PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
282 {
283 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
284
285 PetscFunctionBegin;
286 PetscCallHYPRE(HYPRE_StructMatrixDestroy(ex->hmat));
287 PetscCallHYPRE(HYPRE_StructVectorDestroy(ex->hx));
288 PetscCallHYPRE(HYPRE_StructVectorDestroy(ex->hb));
289 PetscCall(PetscObjectDereference((PetscObject)ex->da));
290 PetscCallMPI(MPI_Comm_free(&ex->hcomm));
291 PetscCall(PetscFree(ex));
292 PetscFunctionReturn(PETSC_SUCCESS);
293 }
294
MatCreate_HYPREStruct(Mat B)295 PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
296 {
297 Mat_HYPREStruct *ex;
298
299 PetscFunctionBegin;
300 PetscCall(PetscNew(&ex));
301 B->data = (void *)ex;
302 B->rmap->bs = 1;
303 B->assembled = PETSC_FALSE;
304
305 B->insertmode = NOT_SET_VALUES;
306
307 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
308 B->ops->mult = MatMult_HYPREStruct;
309 B->ops->zeroentries = MatZeroEntries_HYPREStruct;
310 B->ops->destroy = MatDestroy_HYPREStruct;
311 B->ops->setup = MatSetUp_HYPREStruct;
312
313 ex->needsinitialization = PETSC_TRUE;
314
315 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
316 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT));
317 PetscCall(PetscHYPREInitialize());
318 PetscFunctionReturn(PETSC_SUCCESS);
319 }
320
321 /*MC
322 MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
323 based on the hypre HYPRE_SStructMatrix.
324
325 Level: intermediate
326
327 Notes:
328 Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
329 grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
330
331 Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block
332 be defined by a `DMDA`.
333
334 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
335
336 .seealso: `Mat`
337 M*/
338
MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)339 static PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
340 {
341 HYPRE_Int index[3], *entries;
342 PetscInt i, j, stencil;
343 HYPRE_Complex *values = (HYPRE_Complex *)y;
344 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
345
346 HYPRE_Int part = 0; /* PETSc sstruct interface only allows 1 part */
347 PetscInt ordering;
348 PetscInt grid_rank, to_grid_rank;
349 PetscInt var_type, to_var_type;
350 PetscInt to_var_entry = 0;
351 PetscInt nvars = ex->nvars;
352 PetscInt row;
353
354 PetscFunctionBegin;
355 PetscCall(PetscMalloc1(7 * nvars, &entries));
356 ordering = ex->dofs_order; /* ordering= 0 nodal ordering
357 1 variable ordering */
358 /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
359 if (!ordering) {
360 for (i = 0; i < nrow; i++) {
361 grid_rank = irow[i] / nvars;
362 var_type = (irow[i] % nvars);
363
364 for (j = 0; j < ncol; j++) {
365 to_grid_rank = icol[j] / nvars;
366 to_var_type = (icol[j] % nvars);
367
368 to_var_entry = to_var_entry * 7;
369 entries[j] = (HYPRE_Int)to_var_entry;
370
371 stencil = to_grid_rank - grid_rank;
372 if (!stencil) {
373 entries[j] += 3;
374 } else if (stencil == -1) {
375 entries[j] += 2;
376 } else if (stencil == 1) {
377 entries[j] += 4;
378 } else if (stencil == -ex->gnx) {
379 entries[j] += 1;
380 } else if (stencil == ex->gnx) {
381 entries[j] += 5;
382 } else if (stencil == -ex->gnxgny) {
383 entries[j] += 0;
384 } else if (stencil == ex->gnxgny) {
385 entries[j] += 6;
386 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
387 }
388
389 row = ex->gindices[grid_rank] - ex->rstart;
390 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
391 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
392 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
393
394 if (addv == ADD_VALUES) PetscCallHYPRE(HYPRE_SStructMatrixAddToValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, (HYPRE_Int)ncol, entries, values));
395 else PetscCallHYPRE(HYPRE_SStructMatrixSetValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, (HYPRE_Int)ncol, entries, values));
396 values += ncol;
397 }
398 } else {
399 for (i = 0; i < nrow; i++) {
400 var_type = irow[i] / (ex->gnxgnygnz);
401 grid_rank = irow[i] - var_type * (ex->gnxgnygnz);
402
403 for (j = 0; j < ncol; j++) {
404 to_var_type = icol[j] / (ex->gnxgnygnz);
405 to_grid_rank = icol[j] - to_var_type * (ex->gnxgnygnz);
406
407 to_var_entry = to_var_entry * 7;
408 entries[j] = (HYPRE_Int)to_var_entry;
409
410 stencil = to_grid_rank - grid_rank;
411 if (!stencil) {
412 entries[j] += 3;
413 } else if (stencil == -1) {
414 entries[j] += 2;
415 } else if (stencil == 1) {
416 entries[j] += 4;
417 } else if (stencil == -ex->gnx) {
418 entries[j] += 1;
419 } else if (stencil == ex->gnx) {
420 entries[j] += 5;
421 } else if (stencil == -ex->gnxgny) {
422 entries[j] += 0;
423 } else if (stencil == ex->gnxgny) {
424 entries[j] += 6;
425 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
426 }
427
428 row = ex->gindices[grid_rank] - ex->rstart;
429 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
430 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
431 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
432
433 if (addv == ADD_VALUES) PetscCallHYPRE(HYPRE_SStructMatrixAddToValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, (HYPRE_Int)ncol, entries, values));
434 else PetscCallHYPRE(HYPRE_SStructMatrixSetValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, (HYPRE_Int)ncol, entries, values));
435 values += ncol;
436 }
437 }
438 PetscCall(PetscFree(entries));
439 PetscFunctionReturn(PETSC_SUCCESS);
440 }
441
MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)442 static PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
443 {
444 HYPRE_Int index[3], *entries;
445 PetscInt i;
446 HYPRE_Complex **values;
447 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
448
449 HYPRE_Int part = 0; /* PETSc sstruct interface only allows 1 part */
450 PetscInt ordering = ex->dofs_order;
451 PetscInt grid_rank;
452 PetscInt var_type;
453 PetscInt nvars = ex->nvars;
454 PetscInt row;
455
456 PetscFunctionBegin;
457 PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support");
458 PetscCall(PetscMalloc1(7 * nvars, &entries));
459
460 PetscCall(PetscMalloc1(nvars, &values));
461 PetscCall(PetscMalloc1(7 * nvars * nvars, &values[0]));
462 for (i = 1; i < nvars; i++) values[i] = values[i - 1] + nvars * 7;
463
464 for (i = 0; i < nvars; i++) {
465 PetscCall(PetscArrayzero(values[i], nvars * 7 * sizeof(HYPRE_Complex)));
466 PetscCall(PetscHYPREScalarCast(d, values[i] + 3));
467 }
468
469 for (i = 0; i < nvars * 7; i++) entries[i] = (HYPRE_Int)i;
470
471 if (!ordering) {
472 for (i = 0; i < nrow; i++) {
473 grid_rank = irow[i] / nvars;
474 var_type = (irow[i] % nvars);
475
476 row = ex->gindices[grid_rank] - ex->rstart;
477 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
478 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
479 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
480 PetscCallHYPRE(HYPRE_SStructMatrixSetValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, 7 * (HYPRE_Int)nvars, entries, values[var_type]));
481 }
482 } else {
483 for (i = 0; i < nrow; i++) {
484 var_type = irow[i] / (ex->gnxgnygnz);
485 grid_rank = irow[i] - var_type * (ex->gnxgnygnz);
486
487 row = ex->gindices[grid_rank] - ex->rstart;
488 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
489 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
490 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
491 PetscCallHYPRE(HYPRE_SStructMatrixSetValues(ex->ss_mat, part, index, (HYPRE_Int)var_type, 7 * (HYPRE_Int)nvars, entries, values[var_type]));
492 }
493 }
494 PetscCallHYPRE(HYPRE_SStructMatrixAssemble(ex->ss_mat));
495 PetscCall(PetscFree(values[0]));
496 PetscCall(PetscFree(values));
497 PetscCall(PetscFree(entries));
498 PetscFunctionReturn(PETSC_SUCCESS);
499 }
500
MatZeroEntries_HYPRESStruct_3d(Mat mat)501 static PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
502 {
503 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
504 PetscInt nvars = ex->nvars;
505 PetscInt size;
506 HYPRE_Int part = 0; /* only one part */
507
508 PetscFunctionBegin;
509 size = (ex->hbox.imax[0] - ex->hbox.imin[0] + 1) * (ex->hbox.imax[1] - ex->hbox.imin[1] + 1) * (ex->hbox.imax[2] - ex->hbox.imin[2] + 1);
510 {
511 HYPRE_Int i, *entries, iupper[3], ilower[3];
512 HYPRE_Complex *values;
513
514 for (i = 0; i < 3; i++) {
515 ilower[i] = (HYPRE_Int)ex->hbox.imin[i];
516 iupper[i] = (HYPRE_Int)ex->hbox.imax[i];
517 }
518
519 PetscCall(PetscMalloc2(nvars * 7, &entries, nvars * 7 * size, &values));
520 for (i = 0; i < nvars * 7; i++) entries[i] = i;
521 PetscCall(PetscArrayzero(values, nvars * 7 * size));
522
523 for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructMatrixSetBoxValues(ex->ss_mat, part, ilower, iupper, (HYPRE_Int)i, (HYPRE_Int)nvars * 7, entries, values));
524 PetscCall(PetscFree2(entries, values));
525 }
526 PetscCallHYPRE(HYPRE_SStructMatrixAssemble(ex->ss_mat));
527 PetscFunctionReturn(PETSC_SUCCESS);
528 }
529
MatSetUp_HYPRESStruct(Mat mat)530 static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
531 {
532 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
533 PetscInt dim, dof, sw[3], nx, ny, nz;
534 PetscInt ilower[3], iupper[3], ssize, i;
535 DMBoundaryType px, py, pz;
536 DMDAStencilType st;
537 PetscInt nparts = 1; /* assuming only one part */
538 HYPRE_Int part = 0;
539 ISLocalToGlobalMapping ltog;
540 DM da;
541
542 PetscFunctionBegin;
543 PetscCall(MatGetDM(mat, (DM *)&da));
544 ex->da = da;
545 PetscCall(PetscObjectReference((PetscObject)da));
546
547 PetscCall(DMDAGetInfo(ex->da, &dim, 0, 0, 0, 0, 0, 0, &dof, &sw[0], &px, &py, &pz, &st));
548 PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
549 iupper[0] += ilower[0] - 1;
550 iupper[1] += ilower[1] - 1;
551 iupper[2] += ilower[2] - 1;
552 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
553 ex->hbox.imin[0] = (HYPRE_Int)ilower[0];
554 ex->hbox.imin[1] = (HYPRE_Int)ilower[1];
555 ex->hbox.imin[2] = (HYPRE_Int)ilower[2];
556 ex->hbox.imax[0] = (HYPRE_Int)iupper[0];
557 ex->hbox.imax[1] = (HYPRE_Int)iupper[1];
558 ex->hbox.imax[2] = (HYPRE_Int)iupper[2];
559
560 ex->dofs_order = 0;
561
562 /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
563 ex->nvars = (int)dof;
564
565 /* create the hypre grid object and set its information */
566 PetscCheck(!px && !py && !pz, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
567 PetscCallHYPRE(HYPRE_SStructGridCreate(ex->hcomm, (HYPRE_Int)dim, (HYPRE_Int)nparts, &ex->ss_grid));
568 PetscCallHYPRE(HYPRE_SStructGridSetExtents(ex->ss_grid, part, ex->hbox.imin, ex->hbox.imax));
569 {
570 HYPRE_SStructVariable *vartypes;
571 PetscCall(PetscMalloc1(ex->nvars, &vartypes));
572 for (i = 0; i < ex->nvars; i++) vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL;
573 PetscCallHYPRE(HYPRE_SStructGridSetVariables(ex->ss_grid, part, (HYPRE_Int)ex->nvars, vartypes));
574 PetscCall(PetscFree(vartypes));
575 }
576 PetscCallHYPRE(HYPRE_SStructGridAssemble(ex->ss_grid));
577
578 sw[1] = sw[0];
579 sw[2] = sw[1];
580 /* PetscCallHYPRE(HYPRE_SStructGridSetNumGhost(ex->ss_grid,sw)); */
581
582 /* create the hypre stencil object and set its information */
583 PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils");
584 PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils");
585
586 if (dim == 1) {
587 HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
588 PetscInt j, cnt;
589
590 ssize = 3 * (ex->nvars);
591 PetscCallHYPRE(HYPRE_SStructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->ss_stencil));
592 cnt = 0;
593 for (i = 0; i < (ex->nvars); i++) {
594 for (j = 0; j < 3; j++) {
595 PetscCallHYPRE(HYPRE_SStructStencilSetEntry(ex->ss_stencil, (HYPRE_Int)cnt, offsets[j], (HYPRE_Int)i));
596 cnt++;
597 }
598 }
599 } else if (dim == 2) {
600 HYPRE_Int offsets[5][2] = {
601 {0, -1},
602 {-1, 0 },
603 {0, 0 },
604 {1, 0 },
605 {0, 1 }
606 };
607 PetscInt j, cnt;
608
609 ssize = 5 * (ex->nvars);
610 PetscCallHYPRE(HYPRE_SStructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->ss_stencil));
611 cnt = 0;
612 for (i = 0; i < (ex->nvars); i++) {
613 for (j = 0; j < 5; j++) {
614 PetscCallHYPRE(HYPRE_SStructStencilSetEntry(ex->ss_stencil, (HYPRE_Int)cnt, offsets[j], (HYPRE_Int)i));
615 cnt++;
616 }
617 }
618 } else if (dim == 3) {
619 HYPRE_Int offsets[7][3] = {
620 {0, 0, -1},
621 {0, -1, 0 },
622 {-1, 0, 0 },
623 {0, 0, 0 },
624 {1, 0, 0 },
625 {0, 1, 0 },
626 {0, 0, 1 }
627 };
628 PetscInt j, cnt;
629
630 ssize = 7 * (ex->nvars);
631 PetscCallHYPRE(HYPRE_SStructStencilCreate((HYPRE_Int)dim, (HYPRE_Int)ssize, &ex->ss_stencil));
632 cnt = 0;
633 for (i = 0; i < (ex->nvars); i++) {
634 for (j = 0; j < 7; j++) {
635 PetscCallHYPRE(HYPRE_SStructStencilSetEntry(ex->ss_stencil, (HYPRE_Int)cnt, offsets[j], (HYPRE_Int)i));
636 cnt++;
637 }
638 }
639 }
640
641 /* create the HYPRE graph */
642 PetscCallHYPRE(HYPRE_SStructGraphCreate(ex->hcomm, ex->ss_grid, &ex->ss_graph));
643
644 /* set the stencil graph. Note that each variable has the same graph. This means that each
645 variable couples to all the other variable and with the same stencil pattern. */
646 for (i = 0; i < (ex->nvars); i++) PetscCallHYPRE(HYPRE_SStructGraphSetStencil(ex->ss_graph, part, (HYPRE_Int)i, ex->ss_stencil));
647 PetscCallHYPRE(HYPRE_SStructGraphAssemble(ex->ss_graph));
648
649 /* create the HYPRE sstruct vectors for rhs and solution */
650 PetscCallHYPRE(HYPRE_SStructVectorCreate(ex->hcomm, ex->ss_grid, &ex->ss_b));
651 PetscCallHYPRE(HYPRE_SStructVectorCreate(ex->hcomm, ex->ss_grid, &ex->ss_x));
652 PetscCallHYPRE(HYPRE_SStructVectorInitialize(ex->ss_b));
653 PetscCallHYPRE(HYPRE_SStructVectorInitialize(ex->ss_x));
654 PetscCallHYPRE(HYPRE_SStructVectorAssemble(ex->ss_b));
655 PetscCallHYPRE(HYPRE_SStructVectorAssemble(ex->ss_x));
656
657 /* create the hypre matrix object and set its information */
658 PetscCallHYPRE(HYPRE_SStructMatrixCreate(ex->hcomm, ex->ss_graph, &ex->ss_mat));
659 PetscCallHYPRE(HYPRE_SStructGridDestroy(ex->ss_grid));
660 PetscCallHYPRE(HYPRE_SStructStencilDestroy(ex->ss_stencil));
661 if (ex->needsinitialization) {
662 PetscCallHYPRE(HYPRE_SStructMatrixInitialize(ex->ss_mat));
663 ex->needsinitialization = PETSC_FALSE;
664 }
665
666 /* set the global and local sizes of the matrix */
667 PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz));
668 PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE));
669 PetscCall(PetscLayoutSetBlockSize(mat->rmap, dof));
670 PetscCall(PetscLayoutSetBlockSize(mat->cmap, dof));
671 PetscCall(PetscLayoutSetUp(mat->rmap));
672 PetscCall(PetscLayoutSetUp(mat->cmap));
673 mat->preallocated = PETSC_TRUE;
674
675 PetscCheck(dim == 3, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");
676 mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
677 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
678 mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
679
680 /* PetscCall(MatZeroEntries_HYPRESStruct_3d(mat)); */
681
682 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
683 PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL));
684 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og));
685 PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices));
686 PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, &ex->gnxgnygnz));
687
688 ex->gnxgny *= ex->gnx;
689 ex->gnxgnygnz *= ex->gnxgny;
690
691 PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, &ex->nz));
692
693 ex->nxny = ex->nx * ex->ny;
694 ex->nxnynz = ex->nz * ex->nxny;
695 PetscFunctionReturn(PETSC_SUCCESS);
696 }
697
MatMult_HYPRESStruct(Mat A,Vec x,Vec y)698 static PetscErrorCode MatMult_HYPRESStruct(Mat A, Vec x, Vec y)
699 {
700 const PetscScalar *xx;
701 PetscScalar *yy;
702 HYPRE_Int hlower[3], hupper[3];
703 PetscInt ilower[3], iupper[3];
704 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)A->data;
705 PetscInt ordering = mx->dofs_order;
706 PetscInt nvars = mx->nvars;
707 HYPRE_Int part = 0;
708 PetscInt size;
709 PetscInt i;
710
711 PetscFunctionBegin;
712 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
713
714 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
715 iupper[0] += ilower[0] - 1;
716 iupper[1] += ilower[1] - 1;
717 iupper[2] += ilower[2] - 1;
718 hlower[0] = (HYPRE_Int)ilower[0];
719 hlower[1] = (HYPRE_Int)ilower[1];
720 hlower[2] = (HYPRE_Int)ilower[2];
721 hupper[0] = (HYPRE_Int)iupper[0];
722 hupper[1] = (HYPRE_Int)iupper[1];
723 hupper[2] = (HYPRE_Int)iupper[2];
724
725 size = 1;
726 for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);
727
728 /* copy x values over to hypre for variable ordering */
729 if (ordering) {
730 PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0));
731 PetscCall(VecGetArrayRead(x, &xx));
732 for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(xx + (size * i))));
733 PetscCall(VecRestoreArrayRead(x, &xx));
734 PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b));
735 PetscCallHYPRE(HYPRE_SStructMatrixMatvec(1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x));
736
737 /* copy solution values back to PETSc */
738 PetscCall(VecGetArray(y, &yy));
739 for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(yy + (size * i))));
740 PetscCall(VecRestoreArray(y, &yy));
741 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
742 PetscScalar *z;
743 PetscInt j, k;
744
745 PetscCall(PetscMalloc1(nvars * size, &z));
746 PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0));
747 PetscCall(VecGetArrayRead(x, &xx));
748
749 /* transform nodal to hypre's variable ordering for sys_pfmg */
750 for (i = 0; i < size; i++) {
751 k = i * nvars;
752 for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
753 }
754 for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i))));
755 PetscCall(VecRestoreArrayRead(x, &xx));
756 PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b));
757 PetscCallHYPRE(HYPRE_SStructMatrixMatvec(1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x));
758
759 /* copy solution values back to PETSc */
760 PetscCall(VecGetArray(y, &yy));
761 for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i))));
762 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
763 for (i = 0; i < size; i++) {
764 k = i * nvars;
765 for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
766 }
767 PetscCall(VecRestoreArray(y, &yy));
768 PetscCall(PetscFree(z));
769 }
770 PetscFunctionReturn(PETSC_SUCCESS);
771 }
772
MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)773 static PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode)
774 {
775 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
776
777 PetscFunctionBegin;
778 PetscCallHYPRE(HYPRE_SStructMatrixAssemble(ex->ss_mat));
779 PetscFunctionReturn(PETSC_SUCCESS);
780 }
781
MatZeroEntries_HYPRESStruct(Mat mat)782 static PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
783 {
784 PetscFunctionBegin;
785 /* before the DMDA is set to the matrix the zero doesn't need to do anything */
786 PetscFunctionReturn(PETSC_SUCCESS);
787 }
788
MatDestroy_HYPRESStruct(Mat mat)789 static PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
790 {
791 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
792 ISLocalToGlobalMapping ltog;
793
794 PetscFunctionBegin;
795 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og));
796 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **)&ex->gindices));
797 PetscCallHYPRE(HYPRE_SStructGraphDestroy(ex->ss_graph));
798 PetscCallHYPRE(HYPRE_SStructMatrixDestroy(ex->ss_mat));
799 PetscCallHYPRE(HYPRE_SStructVectorDestroy(ex->ss_x));
800 PetscCallHYPRE(HYPRE_SStructVectorDestroy(ex->ss_b));
801 PetscCall(PetscObjectDereference((PetscObject)ex->da));
802 PetscCallMPI(MPI_Comm_free(&ex->hcomm));
803 PetscCall(PetscFree(ex));
804 PetscFunctionReturn(PETSC_SUCCESS);
805 }
806
MatCreate_HYPRESStruct(Mat B)807 PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
808 {
809 Mat_HYPRESStruct *ex;
810
811 PetscFunctionBegin;
812 PetscCall(PetscNew(&ex));
813 B->data = (void *)ex;
814 B->rmap->bs = 1;
815 B->assembled = PETSC_FALSE;
816
817 B->insertmode = NOT_SET_VALUES;
818
819 B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
820 B->ops->mult = MatMult_HYPRESStruct;
821 B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
822 B->ops->destroy = MatDestroy_HYPRESStruct;
823 B->ops->setup = MatSetUp_HYPRESStruct;
824
825 ex->needsinitialization = PETSC_TRUE;
826
827 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
828 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT));
829 PetscCall(PetscHYPREInitialize());
830 PetscFunctionReturn(PETSC_SUCCESS);
831 }
832