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, &ltog);
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, &ltog);
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, &ltog);
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: }