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