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