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