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;
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) PetscCallExternal(HYPRE_StructMatrixAddToValues, ex->hmat, index, ncol, entries, values);
60: else PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, ncol, entries, values);
61: values += ncol;
62: }
63: return 0;
64: }
66: 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;
74: PetscArrayzero(values, 7);
75: PetscHYPREScalarCast(d, &values[3]);
76: for (i = 0; i < nrow; i++) {
77: row = ex->gindices[irow[i]] - ex->rstart;
78: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
79: index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
80: index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
81: PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, 7, entries, values);
82: }
83: PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
84: return 0;
85: }
87: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
88: {
89: HYPRE_Int indices[7] = {0, 1, 2, 3, 4, 5, 6};
90: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
92: /* hypre has no public interface to do this */
93: PetscCallExternal(hypre_StructMatrixClearBoxValues, ex->hmat, &ex->hbox, 7, indices, 0, 1);
94: PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
95: return 0;
96: }
98: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
99: {
100: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
101: HYPRE_Int sw[6];
102: HYPRE_Int hlower[3], hupper[3], period[3] = {0, 0, 0};
103: PetscInt dim, dof, psw, Nx, Ny, Nz, nx, ny, nz, ilower[3], iupper[3], ssize, i;
104: DMBoundaryType px, py, pz;
105: DMDAStencilType st;
106: ISLocalToGlobalMapping ltog;
107: DM da;
109: MatGetDM(mat, (DM *)&da);
110: ex->da = da;
111: PetscObjectReference((PetscObject)da);
113: DMDAGetInfo(ex->da, &dim, &Nx, &Ny, &Nz, 0, 0, 0, &dof, &psw, &px, &py, &pz, &st);
114: 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 */
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: PetscCallExternal(HYPRE_StructGridCreate, ex->hcomm, dim, &ex->hgrid);
143: PetscCallExternal(HYPRE_StructGridSetExtents, ex->hgrid, hlower, hupper);
144: PetscCallExternal(HYPRE_StructGridSetPeriodic, ex->hgrid, period);
145: PetscCallExternal(HYPRE_StructGridAssemble, ex->hgrid);
147: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw;
148: PetscCallExternal(HYPRE_StructGridSetNumGhost, ex->hgrid, sw);
150: /* create the hypre stencil object and set its information */
153: if (dim == 1) {
154: HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
155: ssize = 3;
156: PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
157: for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, 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: PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
168: for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, 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: PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
181: for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]);
182: }
184: /* create the HYPRE vector for rhs and solution */
185: PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hb);
186: PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hx);
187: PetscCallExternal(HYPRE_StructVectorInitialize, ex->hb);
188: PetscCallExternal(HYPRE_StructVectorInitialize, ex->hx);
189: PetscCallExternal(HYPRE_StructVectorAssemble, ex->hb);
190: PetscCallExternal(HYPRE_StructVectorAssemble, ex->hx);
192: /* create the hypre matrix object and set its information */
193: PetscCallExternal(HYPRE_StructMatrixCreate, ex->hcomm, ex->hgrid, ex->hstencil, &ex->hmat);
194: PetscCallExternal(HYPRE_StructGridDestroy, ex->hgrid);
195: PetscCallExternal(HYPRE_StructStencilDestroy, ex->hstencil);
196: if (ex->needsinitialization) {
197: PetscCallExternal(HYPRE_StructMatrixInitialize, ex->hmat);
198: ex->needsinitialization = PETSC_FALSE;
199: }
201: /* set the global and local sizes of the matrix */
202: DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz);
203: MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE);
204: PetscLayoutSetBlockSize(mat->rmap, 1);
205: PetscLayoutSetBlockSize(mat->cmap, 1);
206: PetscLayoutSetUp(mat->rmap);
207: PetscLayoutSetUp(mat->cmap);
208: mat->preallocated = PETSC_TRUE;
210: if (dim == 3) {
211: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
212: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
213: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
215: /* MatZeroEntries_HYPREStruct_3d(mat); */
216: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");
218: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
219: MatGetOwnershipRange(mat, &ex->rstart, NULL);
220: DMGetLocalToGlobalMapping(ex->da, <og);
221: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices);
222: DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, 0);
223: ex->gnxgny *= ex->gnx;
224: DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, 0);
225: ex->nxny = ex->nx * ex->ny;
226: return 0;
227: }
229: PetscErrorCode MatMult_HYPREStruct(Mat A, Vec x, Vec y)
230: {
231: const PetscScalar *xx;
232: PetscScalar *yy;
233: PetscInt ilower[3], iupper[3];
234: HYPRE_Int hlower[3], hupper[3];
235: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(A->data);
237: 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: PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
251: VecGetArrayRead(x, &xx);
252: PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
253: VecRestoreArrayRead(x, &xx);
254: PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
255: PetscCallExternal(HYPRE_StructMatrixMatvec, 1.0, mx->hmat, mx->hb, 0.0, mx->hx);
257: /* copy solution values back to PETSc */
258: VecGetArray(y, &yy);
259: PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
260: VecRestoreArray(y, &yy);
261: return 0;
262: }
264: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode)
265: {
266: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
268: PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
269: /* PetscCallExternal(HYPRE_StructMatrixPrint,"dummy",ex->hmat,0); */
270: return 0;
271: }
273: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
274: {
275: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
276: return 0;
277: }
279: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
280: {
281: Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;
283: PetscCallExternal(HYPRE_StructMatrixDestroy, ex->hmat);
284: PetscCallExternal(HYPRE_StructVectorDestroy, ex->hx);
285: PetscCallExternal(HYPRE_StructVectorDestroy, ex->hb);
286: PetscObjectDereference((PetscObject)ex->da);
287: MPI_Comm_free(&(ex->hcomm));
288: PetscFree(ex);
289: return 0;
290: }
292: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
293: {
294: Mat_HYPREStruct *ex;
296: PetscNew(&ex);
297: B->data = (void *)ex;
298: B->rmap->bs = 1;
299: B->assembled = PETSC_FALSE;
301: B->insertmode = NOT_SET_VALUES;
303: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
304: B->ops->mult = MatMult_HYPREStruct;
305: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
306: B->ops->destroy = MatDestroy_HYPREStruct;
307: B->ops->setup = MatSetUp_HYPREStruct;
309: ex->needsinitialization = PETSC_TRUE;
311: MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm));
312: PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT);
313: return 0;
314: }
316: /*MC
317: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
318: based on the hypre HYPRE_SStructMatrix.
320: Level: intermediate
322: Notes:
323: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
324: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
326: 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
327: be defined by a `DMDA`.
329: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
331: .seealso: `Mat`
332: M*/
334: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
335: {
336: HYPRE_Int index[3], *entries;
337: PetscInt i, j, stencil;
338: HYPRE_Complex *values = (HYPRE_Complex *)y;
339: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
341: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
342: PetscInt ordering;
343: PetscInt grid_rank, to_grid_rank;
344: PetscInt var_type, to_var_type;
345: PetscInt to_var_entry = 0;
346: PetscInt nvars = ex->nvars;
347: PetscInt row;
349: PetscMalloc1(7 * nvars, &entries);
350: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
351: 1 variable ordering */
352: /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
353: if (!ordering) {
354: for (i = 0; i < nrow; i++) {
355: grid_rank = irow[i] / nvars;
356: var_type = (irow[i] % nvars);
358: for (j = 0; j < ncol; j++) {
359: to_grid_rank = icol[j] / nvars;
360: to_var_type = (icol[j] % nvars);
362: to_var_entry = to_var_entry * 7;
363: entries[j] = to_var_entry;
365: stencil = to_grid_rank - grid_rank;
366: if (!stencil) {
367: entries[j] += 3;
368: } else if (stencil == -1) {
369: entries[j] += 2;
370: } else if (stencil == 1) {
371: entries[j] += 4;
372: } else if (stencil == -ex->gnx) {
373: entries[j] += 1;
374: } else if (stencil == ex->gnx) {
375: entries[j] += 5;
376: } else if (stencil == -ex->gnxgny) {
377: entries[j] += 0;
378: } else if (stencil == ex->gnxgny) {
379: entries[j] += 6;
380: } 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);
381: }
383: row = ex->gindices[grid_rank] - ex->rstart;
384: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
385: index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
386: index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
388: if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
389: else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
390: values += ncol;
391: }
392: } else {
393: for (i = 0; i < nrow; i++) {
394: var_type = irow[i] / (ex->gnxgnygnz);
395: grid_rank = irow[i] - var_type * (ex->gnxgnygnz);
397: for (j = 0; j < ncol; j++) {
398: to_var_type = icol[j] / (ex->gnxgnygnz);
399: to_grid_rank = icol[j] - to_var_type * (ex->gnxgnygnz);
401: to_var_entry = to_var_entry * 7;
402: entries[j] = to_var_entry;
404: stencil = to_grid_rank - grid_rank;
405: if (!stencil) {
406: entries[j] += 3;
407: } else if (stencil == -1) {
408: entries[j] += 2;
409: } else if (stencil == 1) {
410: entries[j] += 4;
411: } else if (stencil == -ex->gnx) {
412: entries[j] += 1;
413: } else if (stencil == ex->gnx) {
414: entries[j] += 5;
415: } else if (stencil == -ex->gnxgny) {
416: entries[j] += 0;
417: } else if (stencil == ex->gnxgny) {
418: entries[j] += 6;
419: } 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);
420: }
422: row = ex->gindices[grid_rank] - ex->rstart;
423: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
424: index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
425: index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
427: if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
428: else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
429: values += ncol;
430: }
431: }
432: PetscFree(entries);
433: return 0;
434: }
436: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
437: {
438: HYPRE_Int index[3], *entries;
439: PetscInt i;
440: HYPRE_Complex **values;
441: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
443: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
444: PetscInt ordering = ex->dofs_order;
445: PetscInt grid_rank;
446: PetscInt var_type;
447: PetscInt nvars = ex->nvars;
448: PetscInt row;
451: PetscMalloc1(7 * nvars, &entries);
453: PetscMalloc1(nvars, &values);
454: PetscMalloc1(7 * nvars * nvars, &values[0]);
455: for (i = 1; i < nvars; i++) values[i] = values[i - 1] + nvars * 7;
457: for (i = 0; i < nvars; i++) {
458: PetscArrayzero(values[i], nvars * 7 * sizeof(HYPRE_Complex));
459: PetscHYPREScalarCast(d, values[i] + 3);
460: }
462: for (i = 0; i < nvars * 7; i++) entries[i] = i;
464: if (!ordering) {
465: for (i = 0; i < nrow; i++) {
466: grid_rank = irow[i] / nvars;
467: var_type = (irow[i] % nvars);
469: row = ex->gindices[grid_rank] - ex->rstart;
470: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
471: index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
472: index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
473: PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]);
474: }
475: } else {
476: for (i = 0; i < nrow; i++) {
477: var_type = irow[i] / (ex->gnxgnygnz);
478: grid_rank = irow[i] - var_type * (ex->gnxgnygnz);
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: }
487: PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
488: PetscFree(values[0]);
489: PetscFree(values);
490: PetscFree(entries);
491: return 0;
492: }
494: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
495: {
496: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
497: PetscInt nvars = ex->nvars;
498: PetscInt size;
499: PetscInt part = 0; /* only one part */
501: 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);
502: {
503: HYPRE_Int i, *entries, iupper[3], ilower[3];
504: HYPRE_Complex *values;
506: for (i = 0; i < 3; i++) {
507: ilower[i] = ex->hbox.imin[i];
508: iupper[i] = ex->hbox.imax[i];
509: }
511: PetscMalloc2(nvars * 7, &entries, nvars * 7 * size, &values);
512: for (i = 0; i < nvars * 7; i++) entries[i] = i;
513: PetscArrayzero(values, nvars * 7 * size);
515: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructMatrixSetBoxValues, ex->ss_mat, part, ilower, iupper, i, nvars * 7, entries, values);
516: PetscFree2(entries, values);
517: }
518: PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
519: return 0;
520: }
522: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
523: {
524: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
525: PetscInt dim, dof, sw[3], nx, ny, nz;
526: PetscInt ilower[3], iupper[3], ssize, i;
527: DMBoundaryType px, py, pz;
528: DMDAStencilType st;
529: PetscInt nparts = 1; /* assuming only one part */
530: PetscInt part = 0;
531: ISLocalToGlobalMapping ltog;
532: DM da;
534: MatGetDM(mat, (DM *)&da);
535: ex->da = da;
536: PetscObjectReference((PetscObject)da);
538: DMDAGetInfo(ex->da, &dim, 0, 0, 0, 0, 0, 0, &dof, &sw[0], &px, &py, &pz, &st);
539: DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]);
540: iupper[0] += ilower[0] - 1;
541: iupper[1] += ilower[1] - 1;
542: iupper[2] += ilower[2] - 1;
543: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
544: ex->hbox.imin[0] = ilower[0];
545: ex->hbox.imin[1] = ilower[1];
546: ex->hbox.imin[2] = ilower[2];
547: ex->hbox.imax[0] = iupper[0];
548: ex->hbox.imax[1] = iupper[1];
549: ex->hbox.imax[2] = iupper[2];
551: ex->dofs_order = 0;
553: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
554: ex->nvars = dof;
556: /* create the hypre grid object and set its information */
558: PetscCallExternal(HYPRE_SStructGridCreate, ex->hcomm, dim, nparts, &ex->ss_grid);
559: PetscCallExternal(HYPRE_SStructGridSetExtents, ex->ss_grid, part, ex->hbox.imin, ex->hbox.imax);
560: {
561: HYPRE_SStructVariable *vartypes;
562: PetscMalloc1(ex->nvars, &vartypes);
563: for (i = 0; i < ex->nvars; i++) vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL;
564: PetscCallExternal(HYPRE_SStructGridSetVariables, ex->ss_grid, part, ex->nvars, vartypes);
565: PetscFree(vartypes);
566: }
567: PetscCallExternal(HYPRE_SStructGridAssemble, ex->ss_grid);
569: sw[1] = sw[0];
570: sw[2] = sw[1];
571: /* PetscCallExternal(HYPRE_SStructGridSetNumGhost,ex->ss_grid,sw); */
573: /* create the hypre stencil object and set its information */
577: if (dim == 1) {
578: HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
579: PetscInt j, cnt;
581: ssize = 3 * (ex->nvars);
582: PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
583: cnt = 0;
584: for (i = 0; i < (ex->nvars); i++) {
585: for (j = 0; j < 3; j++) {
586: PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
587: cnt++;
588: }
589: }
590: } else if (dim == 2) {
591: HYPRE_Int offsets[5][2] = {
592: {0, -1},
593: {-1, 0 },
594: {0, 0 },
595: {1, 0 },
596: {0, 1 }
597: };
598: PetscInt j, cnt;
600: ssize = 5 * (ex->nvars);
601: PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
602: cnt = 0;
603: for (i = 0; i < (ex->nvars); i++) {
604: for (j = 0; j < 5; j++) {
605: PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
606: cnt++;
607: }
608: }
609: } else if (dim == 3) {
610: HYPRE_Int offsets[7][3] = {
611: {0, 0, -1},
612: {0, -1, 0 },
613: {-1, 0, 0 },
614: {0, 0, 0 },
615: {1, 0, 0 },
616: {0, 1, 0 },
617: {0, 0, 1 }
618: };
619: PetscInt j, cnt;
621: ssize = 7 * (ex->nvars);
622: PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
623: cnt = 0;
624: for (i = 0; i < (ex->nvars); i++) {
625: for (j = 0; j < 7; j++) {
626: PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
627: cnt++;
628: }
629: }
630: }
632: /* create the HYPRE graph */
633: PetscCallExternal(HYPRE_SStructGraphCreate, ex->hcomm, ex->ss_grid, &(ex->ss_graph));
635: /* set the stencil graph. Note that each variable has the same graph. This means that each
636: variable couples to all the other variable and with the same stencil pattern. */
637: for (i = 0; i < (ex->nvars); i++) PetscCallExternal(HYPRE_SStructGraphSetStencil, ex->ss_graph, part, i, ex->ss_stencil);
638: PetscCallExternal(HYPRE_SStructGraphAssemble, ex->ss_graph);
640: /* create the HYPRE sstruct vectors for rhs and solution */
641: PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_b);
642: PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_x);
643: PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_b);
644: PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_x);
645: PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_b);
646: PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_x);
648: /* create the hypre matrix object and set its information */
649: PetscCallExternal(HYPRE_SStructMatrixCreate, ex->hcomm, ex->ss_graph, &ex->ss_mat);
650: PetscCallExternal(HYPRE_SStructGridDestroy, ex->ss_grid);
651: PetscCallExternal(HYPRE_SStructStencilDestroy, ex->ss_stencil);
652: if (ex->needsinitialization) {
653: PetscCallExternal(HYPRE_SStructMatrixInitialize, ex->ss_mat);
654: ex->needsinitialization = PETSC_FALSE;
655: }
657: /* set the global and local sizes of the matrix */
658: DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz);
659: MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE);
660: PetscLayoutSetBlockSize(mat->rmap, dof);
661: PetscLayoutSetBlockSize(mat->cmap, dof);
662: PetscLayoutSetUp(mat->rmap);
663: PetscLayoutSetUp(mat->cmap);
664: mat->preallocated = PETSC_TRUE;
666: if (dim == 3) {
667: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
668: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
669: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
671: /* MatZeroEntries_HYPRESStruct_3d(mat); */
672: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");
674: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
675: MatGetOwnershipRange(mat, &ex->rstart, NULL);
676: DMGetLocalToGlobalMapping(ex->da, <og);
677: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices);
678: DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, &ex->gnxgnygnz);
680: ex->gnxgny *= ex->gnx;
681: ex->gnxgnygnz *= ex->gnxgny;
683: DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, &ex->nz);
685: ex->nxny = ex->nx * ex->ny;
686: ex->nxnynz = ex->nz * ex->nxny;
687: return 0;
688: }
690: PetscErrorCode MatMult_HYPRESStruct(Mat A, Vec x, Vec y)
691: {
692: const PetscScalar *xx;
693: PetscScalar *yy;
694: HYPRE_Int hlower[3], hupper[3];
695: PetscInt ilower[3], iupper[3];
696: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(A->data);
697: PetscInt ordering = mx->dofs_order;
698: PetscInt nvars = mx->nvars;
699: PetscInt part = 0;
700: PetscInt size;
701: PetscInt i;
703: DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]);
705: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
706: iupper[0] += ilower[0] - 1;
707: iupper[1] += ilower[1] - 1;
708: iupper[2] += ilower[2] - 1;
709: hlower[0] = (HYPRE_Int)ilower[0];
710: hlower[1] = (HYPRE_Int)ilower[1];
711: hlower[2] = (HYPRE_Int)ilower[2];
712: hupper[0] = (HYPRE_Int)iupper[0];
713: hupper[1] = (HYPRE_Int)iupper[1];
714: hupper[2] = (HYPRE_Int)iupper[2];
716: size = 1;
717: for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);
719: /* copy x values over to hypre for variable ordering */
720: if (ordering) {
721: PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
722: VecGetArrayRead(x, &xx);
723: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i)));
724: VecRestoreArrayRead(x, &xx);
725: PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
726: PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);
728: /* copy solution values back to PETSc */
729: VecGetArray(y, &yy);
730: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i)));
731: VecRestoreArray(y, &yy);
732: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
733: PetscScalar *z;
734: PetscInt j, k;
736: PetscMalloc1(nvars * size, &z);
737: PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
738: VecGetArrayRead(x, &xx);
740: /* transform nodal to hypre's variable ordering for sys_pfmg */
741: for (i = 0; i < size; i++) {
742: k = i * nvars;
743: for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
744: }
745: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
746: VecRestoreArrayRead(x, &xx);
747: PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
748: PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);
750: /* copy solution values back to PETSc */
751: VecGetArray(y, &yy);
752: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
753: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
754: for (i = 0; i < size; i++) {
755: k = i * nvars;
756: for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
757: }
758: VecRestoreArray(y, &yy);
759: PetscFree(z);
760: }
761: return 0;
762: }
764: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode)
765: {
766: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
768: PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
769: return 0;
770: }
772: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
773: {
774: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
775: return 0;
776: }
778: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
779: {
780: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;
781: ISLocalToGlobalMapping ltog;
783: DMGetLocalToGlobalMapping(ex->da, <og);
784: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **)&ex->gindices);
785: PetscCallExternal(HYPRE_SStructGraphDestroy, ex->ss_graph);
786: PetscCallExternal(HYPRE_SStructMatrixDestroy, ex->ss_mat);
787: PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_x);
788: PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_b);
789: PetscObjectDereference((PetscObject)ex->da);
790: MPI_Comm_free(&(ex->hcomm));
791: PetscFree(ex);
792: return 0;
793: }
795: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
796: {
797: Mat_HYPRESStruct *ex;
799: PetscNew(&ex);
800: B->data = (void *)ex;
801: B->rmap->bs = 1;
802: B->assembled = PETSC_FALSE;
804: B->insertmode = NOT_SET_VALUES;
806: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
807: B->ops->mult = MatMult_HYPRESStruct;
808: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
809: B->ops->destroy = MatDestroy_HYPRESStruct;
810: B->ops->setup = MatSetUp_HYPRESStruct;
812: ex->needsinitialization = PETSC_TRUE;
814: MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm));
815: PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT);
816: return 0;
817: }