Actual source code: mhyp.c
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>
8: #include <../src/dm/impls/da/hypre/mhyp.h>
10: /*MC
11: MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
12: based on the hypre HYPRE_StructMatrix.
14: Level: intermediate
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`.
20: The matrix needs a `DMDA` associated with it by either a call to `MatSetDM()` or if the matrix is obtained from `DMCreateMatrix()`
22: .seealso: `MatCreate()`, `PCPFMG`, `MatSetDM()`, `DMCreateMatrix()`
23: M*/
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;
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: }
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;
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: }
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;
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: }
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;
109: PetscFunctionBegin;
110: PetscCall(MatGetDM(mat, (DM *)&da));
111: ex->da = da;
112: PetscCall(PetscObjectReference((PetscObject)da));
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]));
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];
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];
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));
148: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = (HYPRE_Int)psw;
149: PetscCallHYPRE(HYPRE_StructGridSetNumGhost(ex->hgrid, sw));
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: }
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));
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: }
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;
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;
216: /* PetscCall(MatZeroEntries_HYPREStruct_3d(mat)); */
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: }
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;
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];
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));
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: }
265: static PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode)
266: {
267: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
269: PetscFunctionBegin;
270: PetscCallHYPRE(HYPRE_StructMatrixAssemble(ex->hmat));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
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: }
281: static PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
282: {
283: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
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: }
295: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
296: {
297: Mat_HYPREStruct *ex;
299: PetscFunctionBegin;
300: PetscCall(PetscNew(&ex));
301: B->data = (void *)ex;
302: B->rmap->bs = 1;
303: B->assembled = PETSC_FALSE;
305: B->insertmode = NOT_SET_VALUES;
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;
313: ex->needsinitialization = PETSC_TRUE;
315: PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
316: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT));
317: PetscCall(PetscHYPREInitialize());
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*MC
322: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
323: based on the hypre HYPRE_SStructMatrix.
325: Level: intermediate
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.
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`.
334: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
336: .seealso: `Mat`
337: M*/
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;
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;
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);
364: for (j = 0; j < ncol; j++) {
365: to_grid_rank = icol[j] / nvars;
366: to_var_type = (icol[j] % nvars);
368: to_var_entry = to_var_entry * 7;
369: entries[j] = (HYPRE_Int)to_var_entry;
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: }
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)));
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);
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);
407: to_var_entry = to_var_entry * 7;
408: entries[j] = (HYPRE_Int)to_var_entry;
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: }
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)));
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: }
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;
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;
456: PetscFunctionBegin;
457: PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support");
458: PetscCall(PetscMalloc1(7 * nvars, &entries));
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;
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: }
469: for (i = 0; i < nvars * 7; i++) entries[i] = (HYPRE_Int)i;
471: if (!ordering) {
472: for (i = 0; i < nrow; i++) {
473: grid_rank = irow[i] / nvars;
474: var_type = (irow[i] % nvars);
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);
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: }
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 */
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;
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: }
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));
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: }
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;
542: PetscFunctionBegin;
543: PetscCall(MatGetDM(mat, (DM *)&da));
544: ex->da = da;
545: PetscCall(PetscObjectReference((PetscObject)da));
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];
560: ex->dofs_order = 0;
562: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
563: ex->nvars = (int)dof;
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));
578: sw[1] = sw[0];
579: sw[2] = sw[1];
580: /* PetscCallHYPRE(HYPRE_SStructGridSetNumGhost(ex->ss_grid,sw)); */
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");
586: if (dim == 1) {
587: HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
588: PetscInt j, cnt;
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;
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;
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: }
641: /* create the HYPRE graph */
642: PetscCallHYPRE(HYPRE_SStructGraphCreate(ex->hcomm, ex->ss_grid, &ex->ss_graph));
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));
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));
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: }
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;
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;
680: /* PetscCall(MatZeroEntries_HYPRESStruct_3d(mat)); */
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));
688: ex->gnxgny *= ex->gnx;
689: ex->gnxgnygnz *= ex->gnxgny;
691: PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, &ex->nz));
693: ex->nxny = ex->nx * ex->ny;
694: ex->nxnynz = ex->nz * ex->nxny;
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
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;
711: PetscFunctionBegin;
712: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
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];
725: size = 1;
726: for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);
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));
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;
745: PetscCall(PetscMalloc1(nvars * size, &z));
746: PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0));
747: PetscCall(VecGetArrayRead(x, &xx));
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));
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: }
773: static PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode)
774: {
775: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
777: PetscFunctionBegin;
778: PetscCallHYPRE(HYPRE_SStructMatrixAssemble(ex->ss_mat));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
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: }
789: static PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
790: {
791: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
792: ISLocalToGlobalMapping ltog;
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: }
807: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
808: {
809: Mat_HYPRESStruct *ex;
811: PetscFunctionBegin;
812: PetscCall(PetscNew(&ex));
813: B->data = (void *)ex;
814: B->rmap->bs = 1;
815: B->assembled = PETSC_FALSE;
817: B->insertmode = NOT_SET_VALUES;
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;
825: ex->needsinitialization = PETSC_TRUE;
827: PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
828: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT));
829: PetscCall(PetscHYPREInitialize());
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }