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: M*/

333: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
334: {
335:   HYPRE_Int         index[3], *entries;
336:   PetscInt          i, j, stencil;
337:   HYPRE_Complex    *values = (HYPRE_Complex *)y;
338:   Mat_HYPRESStruct *ex     = (Mat_HYPRESStruct *)mat->data;

340:   PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
341:   PetscInt ordering;
342:   PetscInt grid_rank, to_grid_rank;
343:   PetscInt var_type, to_var_type;
344:   PetscInt to_var_entry = 0;
345:   PetscInt nvars        = ex->nvars;
346:   PetscInt row;

348:   PetscMalloc1(7 * nvars, &entries);
349:   ordering = ex->dofs_order; /* ordering= 0   nodal ordering
350:                                            1   variable ordering */
351:   /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ...  */
352:   if (!ordering) {
353:     for (i = 0; i < nrow; i++) {
354:       grid_rank = irow[i] / nvars;
355:       var_type  = (irow[i] % nvars);

357:       for (j = 0; j < ncol; j++) {
358:         to_grid_rank = icol[j] / nvars;
359:         to_var_type  = (icol[j] % nvars);

361:         to_var_entry = to_var_entry * 7;
362:         entries[j]   = to_var_entry;

364:         stencil = to_grid_rank - grid_rank;
365:         if (!stencil) {
366:           entries[j] += 3;
367:         } else if (stencil == -1) {
368:           entries[j] += 2;
369:         } else if (stencil == 1) {
370:           entries[j] += 4;
371:         } else if (stencil == -ex->gnx) {
372:           entries[j] += 1;
373:         } else if (stencil == ex->gnx) {
374:           entries[j] += 5;
375:         } else if (stencil == -ex->gnxgny) {
376:           entries[j] += 0;
377:         } else if (stencil == ex->gnxgny) {
378:           entries[j] += 6;
379:         } 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);
380:       }

382:       row      = ex->gindices[grid_rank] - ex->rstart;
383:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
384:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
385:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));

387:       if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
388:       else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
389:       values += ncol;
390:     }
391:   } else {
392:     for (i = 0; i < nrow; i++) {
393:       var_type  = irow[i] / (ex->gnxgnygnz);
394:       grid_rank = irow[i] - var_type * (ex->gnxgnygnz);

396:       for (j = 0; j < ncol; j++) {
397:         to_var_type  = icol[j] / (ex->gnxgnygnz);
398:         to_grid_rank = icol[j] - to_var_type * (ex->gnxgnygnz);

400:         to_var_entry = to_var_entry * 7;
401:         entries[j]   = to_var_entry;

403:         stencil = to_grid_rank - grid_rank;
404:         if (!stencil) {
405:           entries[j] += 3;
406:         } else if (stencil == -1) {
407:           entries[j] += 2;
408:         } else if (stencil == 1) {
409:           entries[j] += 4;
410:         } else if (stencil == -ex->gnx) {
411:           entries[j] += 1;
412:         } else if (stencil == ex->gnx) {
413:           entries[j] += 5;
414:         } else if (stencil == -ex->gnxgny) {
415:           entries[j] += 0;
416:         } else if (stencil == ex->gnxgny) {
417:           entries[j] += 6;
418:         } 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);
419:       }

421:       row      = ex->gindices[grid_rank] - ex->rstart;
422:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
423:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
424:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));

426:       if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
427:       else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
428:       values += ncol;
429:     }
430:   }
431:   PetscFree(entries);
432:   return 0;
433: }

435: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
436: {
437:   HYPRE_Int         index[3], *entries;
438:   PetscInt          i;
439:   HYPRE_Complex   **values;
440:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;

442:   PetscInt part     = 0; /* Petsc sstruct interface only allows 1 part */
443:   PetscInt ordering = ex->dofs_order;
444:   PetscInt grid_rank;
445:   PetscInt var_type;
446:   PetscInt nvars = ex->nvars;
447:   PetscInt row;

450:   PetscMalloc1(7 * nvars, &entries);

452:   PetscMalloc1(nvars, &values);
453:   PetscMalloc1(7 * nvars * nvars, &values[0]);
454:   for (i = 1; i < nvars; i++) values[i] = values[i - 1] + nvars * 7;

456:   for (i = 0; i < nvars; i++) {
457:     PetscArrayzero(values[i], nvars * 7 * sizeof(HYPRE_Complex));
458:     PetscHYPREScalarCast(d, values[i] + 3);
459:   }

461:   for (i = 0; i < nvars * 7; i++) entries[i] = i;

463:   if (!ordering) {
464:     for (i = 0; i < nrow; i++) {
465:       grid_rank = irow[i] / nvars;
466:       var_type  = (irow[i] % nvars);

468:       row      = ex->gindices[grid_rank] - ex->rstart;
469:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
470:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
471:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
472:       PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]);
473:     }
474:   } else {
475:     for (i = 0; i < nrow; i++) {
476:       var_type  = irow[i] / (ex->gnxgnygnz);
477:       grid_rank = irow[i] - var_type * (ex->gnxgnygnz);

479:       row      = ex->gindices[grid_rank] - ex->rstart;
480:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
481:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
482:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
483:       PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]);
484:     }
485:   }
486:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
487:   PetscFree(values[0]);
488:   PetscFree(values);
489:   PetscFree(entries);
490:   return 0;
491: }

493: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
494: {
495:   Mat_HYPRESStruct *ex    = (Mat_HYPRESStruct *)mat->data;
496:   PetscInt          nvars = ex->nvars;
497:   PetscInt          size;
498:   PetscInt          part = 0; /* only one part */

500:   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);
501:   {
502:     HYPRE_Int      i, *entries, iupper[3], ilower[3];
503:     HYPRE_Complex *values;

505:     for (i = 0; i < 3; i++) {
506:       ilower[i] = ex->hbox.imin[i];
507:       iupper[i] = ex->hbox.imax[i];
508:     }

510:     PetscMalloc2(nvars * 7, &entries, nvars * 7 * size, &values);
511:     for (i = 0; i < nvars * 7; i++) entries[i] = i;
512:     PetscArrayzero(values, nvars * 7 * size);

514:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructMatrixSetBoxValues, ex->ss_mat, part, ilower, iupper, i, nvars * 7, entries, values);
515:     PetscFree2(entries, values);
516:   }
517:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
518:   return 0;
519: }

521: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
522: {
523:   Mat_HYPRESStruct      *ex = (Mat_HYPRESStruct *)mat->data;
524:   PetscInt               dim, dof, sw[3], nx, ny, nz;
525:   PetscInt               ilower[3], iupper[3], ssize, i;
526:   DMBoundaryType         px, py, pz;
527:   DMDAStencilType        st;
528:   PetscInt               nparts = 1; /* assuming only one part */
529:   PetscInt               part   = 0;
530:   ISLocalToGlobalMapping ltog;
531:   DM                     da;

533:   MatGetDM(mat, (DM *)&da);
534:   ex->da = da;
535:   PetscObjectReference((PetscObject)da);

537:   DMDAGetInfo(ex->da, &dim, 0, 0, 0, 0, 0, 0, &dof, &sw[0], &px, &py, &pz, &st);
538:   DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]);
539:   iupper[0] += ilower[0] - 1;
540:   iupper[1] += ilower[1] - 1;
541:   iupper[2] += ilower[2] - 1;
542:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
543:   ex->hbox.imin[0] = ilower[0];
544:   ex->hbox.imin[1] = ilower[1];
545:   ex->hbox.imin[2] = ilower[2];
546:   ex->hbox.imax[0] = iupper[0];
547:   ex->hbox.imax[1] = iupper[1];
548:   ex->hbox.imax[2] = iupper[2];

550:   ex->dofs_order = 0;

552:   /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
553:   ex->nvars = dof;

555:   /* create the hypre grid object and set its information */
557:   PetscCallExternal(HYPRE_SStructGridCreate, ex->hcomm, dim, nparts, &ex->ss_grid);
558:   PetscCallExternal(HYPRE_SStructGridSetExtents, ex->ss_grid, part, ex->hbox.imin, ex->hbox.imax);
559:   {
560:     HYPRE_SStructVariable *vartypes;
561:     PetscMalloc1(ex->nvars, &vartypes);
562:     for (i = 0; i < ex->nvars; i++) vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL;
563:     PetscCallExternal(HYPRE_SStructGridSetVariables, ex->ss_grid, part, ex->nvars, vartypes);
564:     PetscFree(vartypes);
565:   }
566:   PetscCallExternal(HYPRE_SStructGridAssemble, ex->ss_grid);

568:   sw[1] = sw[0];
569:   sw[2] = sw[1];
570:   /* PetscCallExternal(HYPRE_SStructGridSetNumGhost,ex->ss_grid,sw); */

572:   /* create the hypre stencil object and set its information */

576:   if (dim == 1) {
577:     HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
578:     PetscInt  j, cnt;

580:     ssize = 3 * (ex->nvars);
581:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
582:     cnt = 0;
583:     for (i = 0; i < (ex->nvars); i++) {
584:       for (j = 0; j < 3; j++) {
585:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
586:         cnt++;
587:       }
588:     }
589:   } else if (dim == 2) {
590:     HYPRE_Int offsets[5][2] = {
591:       {0,  -1},
592:       {-1, 0 },
593:       {0,  0 },
594:       {1,  0 },
595:       {0,  1 }
596:     };
597:     PetscInt j, cnt;

599:     ssize = 5 * (ex->nvars);
600:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
601:     cnt = 0;
602:     for (i = 0; i < (ex->nvars); i++) {
603:       for (j = 0; j < 5; j++) {
604:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
605:         cnt++;
606:       }
607:     }
608:   } else if (dim == 3) {
609:     HYPRE_Int offsets[7][3] = {
610:       {0,  0,  -1},
611:       {0,  -1, 0 },
612:       {-1, 0,  0 },
613:       {0,  0,  0 },
614:       {1,  0,  0 },
615:       {0,  1,  0 },
616:       {0,  0,  1 }
617:     };
618:     PetscInt j, cnt;

620:     ssize = 7 * (ex->nvars);
621:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
622:     cnt = 0;
623:     for (i = 0; i < (ex->nvars); i++) {
624:       for (j = 0; j < 7; j++) {
625:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
626:         cnt++;
627:       }
628:     }
629:   }

631:   /* create the HYPRE graph */
632:   PetscCallExternal(HYPRE_SStructGraphCreate, ex->hcomm, ex->ss_grid, &(ex->ss_graph));

634:   /* set the stencil graph. Note that each variable has the same graph. This means that each
635:      variable couples to all the other variable and with the same stencil pattern. */
636:   for (i = 0; i < (ex->nvars); i++) PetscCallExternal(HYPRE_SStructGraphSetStencil, ex->ss_graph, part, i, ex->ss_stencil);
637:   PetscCallExternal(HYPRE_SStructGraphAssemble, ex->ss_graph);

639:   /* create the HYPRE sstruct vectors for rhs and solution */
640:   PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_b);
641:   PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_x);
642:   PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_b);
643:   PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_x);
644:   PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_b);
645:   PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_x);

647:   /* create the hypre matrix object and set its information */
648:   PetscCallExternal(HYPRE_SStructMatrixCreate, ex->hcomm, ex->ss_graph, &ex->ss_mat);
649:   PetscCallExternal(HYPRE_SStructGridDestroy, ex->ss_grid);
650:   PetscCallExternal(HYPRE_SStructStencilDestroy, ex->ss_stencil);
651:   if (ex->needsinitialization) {
652:     PetscCallExternal(HYPRE_SStructMatrixInitialize, ex->ss_mat);
653:     ex->needsinitialization = PETSC_FALSE;
654:   }

656:   /* set the global and local sizes of the matrix */
657:   DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz);
658:   MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE);
659:   PetscLayoutSetBlockSize(mat->rmap, dof);
660:   PetscLayoutSetBlockSize(mat->cmap, dof);
661:   PetscLayoutSetUp(mat->rmap);
662:   PetscLayoutSetUp(mat->cmap);
663:   mat->preallocated = PETSC_TRUE;

665:   if (dim == 3) {
666:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
667:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
668:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

670:     /* MatZeroEntries_HYPRESStruct_3d(mat); */
671:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");

673:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
674:   MatGetOwnershipRange(mat, &ex->rstart, NULL);
675:   DMGetLocalToGlobalMapping(ex->da, &ltog);
676:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices);
677:   DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, &ex->gnxgnygnz);

679:   ex->gnxgny *= ex->gnx;
680:   ex->gnxgnygnz *= ex->gnxgny;

682:   DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, &ex->nz);

684:   ex->nxny   = ex->nx * ex->ny;
685:   ex->nxnynz = ex->nz * ex->nxny;
686:   return 0;
687: }

689: PetscErrorCode MatMult_HYPRESStruct(Mat A, Vec x, Vec y)
690: {
691:   const PetscScalar *xx;
692:   PetscScalar       *yy;
693:   HYPRE_Int          hlower[3], hupper[3];
694:   PetscInt           ilower[3], iupper[3];
695:   Mat_HYPRESStruct  *mx       = (Mat_HYPRESStruct *)(A->data);
696:   PetscInt           ordering = mx->dofs_order;
697:   PetscInt           nvars    = mx->nvars;
698:   PetscInt           part     = 0;
699:   PetscInt           size;
700:   PetscInt           i;

702:   DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]);

704:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
705:   iupper[0] += ilower[0] - 1;
706:   iupper[1] += ilower[1] - 1;
707:   iupper[2] += ilower[2] - 1;
708:   hlower[0] = (HYPRE_Int)ilower[0];
709:   hlower[1] = (HYPRE_Int)ilower[1];
710:   hlower[2] = (HYPRE_Int)ilower[2];
711:   hupper[0] = (HYPRE_Int)iupper[0];
712:   hupper[1] = (HYPRE_Int)iupper[1];
713:   hupper[2] = (HYPRE_Int)iupper[2];

715:   size = 1;
716:   for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);

718:   /* copy x values over to hypre for variable ordering */
719:   if (ordering) {
720:     PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
721:     VecGetArrayRead(x, &xx);
722:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i)));
723:     VecRestoreArrayRead(x, &xx);
724:     PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
725:     PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);

727:     /* copy solution values back to PETSc */
728:     VecGetArray(y, &yy);
729:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i)));
730:     VecRestoreArray(y, &yy);
731:   } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
732:     PetscScalar *z;
733:     PetscInt     j, k;

735:     PetscMalloc1(nvars * size, &z);
736:     PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
737:     VecGetArrayRead(x, &xx);

739:     /* transform nodal to hypre's variable ordering for sys_pfmg */
740:     for (i = 0; i < size; i++) {
741:       k = i * nvars;
742:       for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
743:     }
744:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
745:     VecRestoreArrayRead(x, &xx);
746:     PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
747:     PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);

749:     /* copy solution values back to PETSc */
750:     VecGetArray(y, &yy);
751:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
752:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
753:     for (i = 0; i < size; i++) {
754:       k = i * nvars;
755:       for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
756:     }
757:     VecRestoreArray(y, &yy);
758:     PetscFree(z);
759:   }
760:   return 0;
761: }

763: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode)
764: {
765:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;

767:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
768:   return 0;
769: }

771: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
772: {
773:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
774:   return 0;
775: }

777: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
778: {
779:   Mat_HYPRESStruct      *ex = (Mat_HYPRESStruct *)mat->data;
780:   ISLocalToGlobalMapping ltog;

782:   DMGetLocalToGlobalMapping(ex->da, &ltog);
783:   ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **)&ex->gindices);
784:   PetscCallExternal(HYPRE_SStructGraphDestroy, ex->ss_graph);
785:   PetscCallExternal(HYPRE_SStructMatrixDestroy, ex->ss_mat);
786:   PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_x);
787:   PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_b);
788:   PetscObjectDereference((PetscObject)ex->da);
789:   MPI_Comm_free(&(ex->hcomm));
790:   PetscFree(ex);
791:   return 0;
792: }

794: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
795: {
796:   Mat_HYPRESStruct *ex;

798:   PetscNew(&ex);
799:   B->data      = (void *)ex;
800:   B->rmap->bs  = 1;
801:   B->assembled = PETSC_FALSE;

803:   B->insertmode = NOT_SET_VALUES;

805:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
806:   B->ops->mult        = MatMult_HYPRESStruct;
807:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
808:   B->ops->destroy     = MatDestroy_HYPRESStruct;
809:   B->ops->setup       = MatSetUp_HYPRESStruct;

811:   ex->needsinitialization = PETSC_TRUE;

813:   MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm));
814:   PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT);
815:   return 0;
816: }