Actual source code: mhyp.c

      1: /*
      2:     Creates hypre ijmatrix from PETSc matrix
      3: */
      4: #include <petscsys.h>
      5: #include <petsc/private/petschypre.h>
      6: #include <petsc/private/matimpl.h>
      7: #include <petscdmda.h>
      8: #include <../src/dm/impls/da/hypre/mhyp.h>

     10: /* -----------------------------------------------------------------------------------------------------------------*/

     12: /*MC
     13:    MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
     14:           based on the hypre HYPRE_StructMatrix.

     16:    Level: intermediate

     18:    Notes:
     19:     Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
     20:           be defined by a `DMDA`.

     22:           The matrix needs a `DMDA` associated with it by either a call to `MatSetDM()` or if the matrix is obtained from `DMCreateMatrix()`

     24: .seealso: `MatCreate()`, `PCPFMG`, `MatSetDM()`, `DMCreateMatrix()`
     25: M*/

     27: static PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
     28: {
     29:   HYPRE_Int        index[3], entries[9];
     30:   PetscInt         i, j, stencil, row;
     31:   HYPRE_Complex   *values = (HYPRE_Complex *)y;
     32:   Mat_HYPREStruct *ex     = (Mat_HYPREStruct *)mat->data;

     34:   PetscFunctionBegin;
     35:   PetscCheck(ncol <= 9, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncol %" PetscInt_FMT " > 9 too large", ncol);
     36:   for (i = 0; i < nrow; i++) {
     37:     for (j = 0; j < ncol; j++) {
     38:       stencil = icol[j] - irow[i];
     39:       if (!stencil) {
     40:         entries[j] = 3;
     41:       } else if (stencil == -1) {
     42:         entries[j] = 2;
     43:       } else if (stencil == 1) {
     44:         entries[j] = 4;
     45:       } else if (stencil == -ex->gnx) {
     46:         entries[j] = 1;
     47:       } else if (stencil == ex->gnx) {
     48:         entries[j] = 5;
     49:       } else if (stencil == -ex->gnxgny) {
     50:         entries[j] = 0;
     51:       } else if (stencil == ex->gnxgny) {
     52:         entries[j] = 6;
     53:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
     54:     }
     55:     row      = ex->gindices[irow[i]] - ex->rstart;
     56:     index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
     57:     index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
     58:     index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
     59:     if (addv == ADD_VALUES) 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:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: static PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
     67: {
     68:   HYPRE_Int        index[3], entries[7] = {0, 1, 2, 3, 4, 5, 6};
     69:   PetscInt         row, i;
     70:   HYPRE_Complex    values[7];
     71:   Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;

     73:   PetscFunctionBegin;
     74:   PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support");
     75:   PetscCall(PetscArrayzero(values, 7));
     76:   PetscCall(PetscHYPREScalarCast(d, &values[3]));
     77:   for (i = 0; i < nrow; i++) {
     78:     row      = ex->gindices[irow[i]] - ex->rstart;
     79:     index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
     80:     index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
     81:     index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
     82:     PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, 7, entries, values);
     83:   }
     84:   PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: static PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
     89: {
     90:   HYPRE_Int        indices[7] = {0, 1, 2, 3, 4, 5, 6};
     91:   Mat_HYPREStruct *ex         = (Mat_HYPREStruct *)mat->data;

     93:   PetscFunctionBegin;
     94:   /* hypre has no public interface to do this */
     95:   PetscCallExternal(hypre_StructMatrixClearBoxValues, ex->hmat, &ex->hbox, 7, indices, 0, 1);
     96:   PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
    101: {
    102:   Mat_HYPREStruct       *ex = (Mat_HYPREStruct *)mat->data;
    103:   HYPRE_Int              sw[6];
    104:   HYPRE_Int              hlower[3], hupper[3], period[3] = {0, 0, 0};
    105:   PetscInt               dim, dof, psw, Nx, Ny, Nz, nx, ny, nz, ilower[3], iupper[3], ssize, i;
    106:   DMBoundaryType         px, py, pz;
    107:   DMDAStencilType        st;
    108:   ISLocalToGlobalMapping ltog;
    109:   DM                     da;

    111:   PetscFunctionBegin;
    112:   PetscCall(MatGetDM(mat, (DM *)&da));
    113:   ex->da = da;
    114:   PetscCall(PetscObjectReference((PetscObject)da));

    116:   PetscCall(DMDAGetInfo(ex->da, &dim, &Nx, &Ny, &Nz, 0, 0, 0, &dof, &psw, &px, &py, &pz, &st));
    117:   PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));

    119:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
    120:   iupper[0] += ilower[0] - 1;
    121:   iupper[1] += ilower[1] - 1;
    122:   iupper[2] += ilower[2] - 1;
    123:   hlower[0] = (HYPRE_Int)ilower[0];
    124:   hlower[1] = (HYPRE_Int)ilower[1];
    125:   hlower[2] = (HYPRE_Int)ilower[2];
    126:   hupper[0] = (HYPRE_Int)iupper[0];
    127:   hupper[1] = (HYPRE_Int)iupper[1];
    128:   hupper[2] = (HYPRE_Int)iupper[2];

    130:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
    131:   ex->hbox.imin[0] = hlower[0];
    132:   ex->hbox.imin[1] = hlower[1];
    133:   ex->hbox.imin[2] = hlower[2];
    134:   ex->hbox.imax[0] = hupper[0];
    135:   ex->hbox.imax[1] = hupper[1];
    136:   ex->hbox.imax[2] = hupper[2];

    138:   /* create the hypre grid object and set its information */
    139:   PetscCheck(dof <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Currently only support for scalar problems");
    140:   if (px || py || pz) {
    141:     if (px == DM_BOUNDARY_PERIODIC) period[0] = (HYPRE_Int)Nx;
    142:     if (py == DM_BOUNDARY_PERIODIC) period[1] = (HYPRE_Int)Ny;
    143:     if (pz == DM_BOUNDARY_PERIODIC) period[2] = (HYPRE_Int)Nz;
    144:   }
    145:   PetscCallExternal(HYPRE_StructGridCreate, ex->hcomm, dim, &ex->hgrid);
    146:   PetscCallExternal(HYPRE_StructGridSetExtents, ex->hgrid, hlower, hupper);
    147:   PetscCallExternal(HYPRE_StructGridSetPeriodic, ex->hgrid, period);
    148:   PetscCallExternal(HYPRE_StructGridAssemble, ex->hgrid);

    150:   sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw;
    151:   PetscCallExternal(HYPRE_StructGridSetNumGhost, ex->hgrid, sw);

    153:   /* create the hypre stencil object and set its information */
    154:   PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils");
    155:   PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils");
    156:   if (dim == 1) {
    157:     HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
    158:     ssize                   = 3;
    159:     PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
    160:     for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]);
    161:   } else if (dim == 2) {
    162:     HYPRE_Int offsets[5][2] = {
    163:       {0,  -1},
    164:       {-1, 0 },
    165:       {0,  0 },
    166:       {1,  0 },
    167:       {0,  1 }
    168:     };
    169:     ssize = 5;
    170:     PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
    171:     for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]);
    172:   } else if (dim == 3) {
    173:     HYPRE_Int offsets[7][3] = {
    174:       {0,  0,  -1},
    175:       {0,  -1, 0 },
    176:       {-1, 0,  0 },
    177:       {0,  0,  0 },
    178:       {1,  0,  0 },
    179:       {0,  1,  0 },
    180:       {0,  0,  1 }
    181:     };
    182:     ssize = 7;
    183:     PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil);
    184:     for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]);
    185:   }

    187:   /* create the HYPRE vector for rhs and solution */
    188:   PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hb);
    189:   PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hx);
    190:   PetscCallExternal(HYPRE_StructVectorInitialize, ex->hb);
    191:   PetscCallExternal(HYPRE_StructVectorInitialize, ex->hx);
    192:   PetscCallExternal(HYPRE_StructVectorAssemble, ex->hb);
    193:   PetscCallExternal(HYPRE_StructVectorAssemble, ex->hx);

    195:   /* create the hypre matrix object and set its information */
    196:   PetscCallExternal(HYPRE_StructMatrixCreate, ex->hcomm, ex->hgrid, ex->hstencil, &ex->hmat);
    197:   PetscCallExternal(HYPRE_StructGridDestroy, ex->hgrid);
    198:   PetscCallExternal(HYPRE_StructStencilDestroy, ex->hstencil);
    199:   if (ex->needsinitialization) {
    200:     PetscCallExternal(HYPRE_StructMatrixInitialize, ex->hmat);
    201:     ex->needsinitialization = PETSC_FALSE;
    202:   }

    204:   /* set the global and local sizes of the matrix */
    205:   PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz));
    206:   PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE));
    207:   PetscCall(PetscLayoutSetBlockSize(mat->rmap, 1));
    208:   PetscCall(PetscLayoutSetBlockSize(mat->cmap, 1));
    209:   PetscCall(PetscLayoutSetUp(mat->rmap));
    210:   PetscCall(PetscLayoutSetUp(mat->cmap));
    211:   mat->preallocated = PETSC_TRUE;

    213:   if (dim == 3) {
    214:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
    215:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPREStruct_3d;
    216:     mat->ops->zeroentries    = MatZeroEntries_HYPREStruct_3d;

    218:     /*        PetscCall(MatZeroEntries_HYPREStruct_3d(mat)); */
    219:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently");

    221:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
    222:   PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL));
    223:   PetscCall(DMGetLocalToGlobalMapping(ex->da, &ltog));
    224:   PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices));
    225:   PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, 0));
    226:   ex->gnxgny *= ex->gnx;
    227:   PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, 0));
    228:   ex->nxny = ex->nx * ex->ny;
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: static PetscErrorCode MatMult_HYPREStruct(Mat A, Vec x, Vec y)
    233: {
    234:   const PetscScalar *xx;
    235:   PetscScalar       *yy;
    236:   PetscInt           ilower[3], iupper[3];
    237:   HYPRE_Int          hlower[3], hupper[3];
    238:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)A->data;

    240:   PetscFunctionBegin;
    241:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
    242:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
    243:   iupper[0] += ilower[0] - 1;
    244:   iupper[1] += ilower[1] - 1;
    245:   iupper[2] += ilower[2] - 1;
    246:   hlower[0] = (HYPRE_Int)ilower[0];
    247:   hlower[1] = (HYPRE_Int)ilower[1];
    248:   hlower[2] = (HYPRE_Int)ilower[2];
    249:   hupper[0] = (HYPRE_Int)iupper[0];
    250:   hupper[1] = (HYPRE_Int)iupper[1];
    251:   hupper[2] = (HYPRE_Int)iupper[2];

    253:   /* copy x values over to hypre */
    254:   PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
    255:   PetscCall(VecGetArrayRead(x, &xx));
    256:   PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
    257:   PetscCall(VecRestoreArrayRead(x, &xx));
    258:   PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
    259:   PetscCallExternal(HYPRE_StructMatrixMatvec, 1.0, mx->hmat, mx->hb, 0.0, mx->hx);

    261:   /* copy solution values back to PETSc */
    262:   PetscCall(VecGetArray(y, &yy));
    263:   PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
    264:   PetscCall(VecRestoreArray(y, &yy));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: static PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode)
    269: {
    270:   Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;

    272:   PetscFunctionBegin;
    273:   PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat);
    274:   /* PetscCallExternal(HYPRE_StructMatrixPrint,"dummy",ex->hmat,0); */
    275:   PetscFunctionReturn(PETSC_SUCCESS);
    276: }

    278: static PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
    279: {
    280:   PetscFunctionBegin;
    281:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: static PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
    286: {
    287:   Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data;

    289:   PetscFunctionBegin;
    290:   PetscCallExternal(HYPRE_StructMatrixDestroy, ex->hmat);
    291:   PetscCallExternal(HYPRE_StructVectorDestroy, ex->hx);
    292:   PetscCallExternal(HYPRE_StructVectorDestroy, ex->hb);
    293:   PetscCall(PetscObjectDereference((PetscObject)ex->da));
    294:   PetscCallMPI(MPI_Comm_free(&ex->hcomm));
    295:   PetscCall(PetscFree(ex));
    296:   PetscFunctionReturn(PETSC_SUCCESS);
    297: }

    299: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
    300: {
    301:   Mat_HYPREStruct *ex;

    303:   PetscFunctionBegin;
    304:   PetscCall(PetscNew(&ex));
    305:   B->data      = (void *)ex;
    306:   B->rmap->bs  = 1;
    307:   B->assembled = PETSC_FALSE;

    309:   B->insertmode = NOT_SET_VALUES;

    311:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
    312:   B->ops->mult        = MatMult_HYPREStruct;
    313:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
    314:   B->ops->destroy     = MatDestroy_HYPREStruct;
    315:   B->ops->setup       = MatSetUp_HYPREStruct;

    317:   ex->needsinitialization = PETSC_TRUE;

    319:   PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
    320:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT));
    321:   PetscHYPREInitialize();
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: /*MC
    326:    MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
    327:           based on the hypre HYPRE_SStructMatrix.

    329:    Level: intermediate

    331:    Notes:
    332:     Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
    333:           grid objects, we restrict the semi-struct objects to consist of only structured-grid components.

    335:           Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block
    336:           be defined by a `DMDA`.

    338:           The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()

    340: .seealso: `Mat`
    341: M*/

    343: static PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv)
    344: {
    345:   HYPRE_Int         index[3], *entries;
    346:   PetscInt          i, j, stencil;
    347:   HYPRE_Complex    *values = (HYPRE_Complex *)y;
    348:   Mat_HYPRESStruct *ex     = (Mat_HYPRESStruct *)mat->data;

    350:   PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
    351:   PetscInt ordering;
    352:   PetscInt grid_rank, to_grid_rank;
    353:   PetscInt var_type, to_var_type;
    354:   PetscInt to_var_entry = 0;
    355:   PetscInt nvars        = ex->nvars;
    356:   PetscInt row;

    358:   PetscFunctionBegin;
    359:   PetscCall(PetscMalloc1(7 * nvars, &entries));
    360:   ordering = ex->dofs_order; /* ordering= 0   nodal ordering
    361:                                            1   variable ordering */
    362:   /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ...  */
    363:   if (!ordering) {
    364:     for (i = 0; i < nrow; i++) {
    365:       grid_rank = irow[i] / nvars;
    366:       var_type  = (irow[i] % nvars);

    368:       for (j = 0; j < ncol; j++) {
    369:         to_grid_rank = icol[j] / nvars;
    370:         to_var_type  = (icol[j] % nvars);

    372:         to_var_entry = to_var_entry * 7;
    373:         entries[j]   = to_var_entry;

    375:         stencil = to_grid_rank - grid_rank;
    376:         if (!stencil) {
    377:           entries[j] += 3;
    378:         } else if (stencil == -1) {
    379:           entries[j] += 2;
    380:         } else if (stencil == 1) {
    381:           entries[j] += 4;
    382:         } else if (stencil == -ex->gnx) {
    383:           entries[j] += 1;
    384:         } else if (stencil == ex->gnx) {
    385:           entries[j] += 5;
    386:         } else if (stencil == -ex->gnxgny) {
    387:           entries[j] += 0;
    388:         } else if (stencil == ex->gnxgny) {
    389:           entries[j] += 6;
    390:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
    391:       }

    393:       row      = ex->gindices[grid_rank] - ex->rstart;
    394:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
    395:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
    396:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));

    398:       if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
    399:       else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
    400:       values += ncol;
    401:     }
    402:   } else {
    403:     for (i = 0; i < nrow; i++) {
    404:       var_type  = irow[i] / (ex->gnxgnygnz);
    405:       grid_rank = irow[i] - var_type * (ex->gnxgnygnz);

    407:       for (j = 0; j < ncol; j++) {
    408:         to_var_type  = icol[j] / (ex->gnxgnygnz);
    409:         to_grid_rank = icol[j] - to_var_type * (ex->gnxgnygnz);

    411:         to_var_entry = to_var_entry * 7;
    412:         entries[j]   = to_var_entry;

    414:         stencil = to_grid_rank - grid_rank;
    415:         if (!stencil) {
    416:           entries[j] += 3;
    417:         } else if (stencil == -1) {
    418:           entries[j] += 2;
    419:         } else if (stencil == 1) {
    420:           entries[j] += 4;
    421:         } else if (stencil == -ex->gnx) {
    422:           entries[j] += 1;
    423:         } else if (stencil == ex->gnx) {
    424:           entries[j] += 5;
    425:         } else if (stencil == -ex->gnxgny) {
    426:           entries[j] += 0;
    427:         } else if (stencil == ex->gnxgny) {
    428:           entries[j] += 6;
    429:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil);
    430:       }

    432:       row      = ex->gindices[grid_rank] - ex->rstart;
    433:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
    434:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
    435:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));

    437:       if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
    438:       else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values);
    439:       values += ncol;
    440:     }
    441:   }
    442:   PetscCall(PetscFree(entries));
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: static PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b)
    447: {
    448:   HYPRE_Int         index[3], *entries;
    449:   PetscInt          i;
    450:   HYPRE_Complex   **values;
    451:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;

    453:   PetscInt part     = 0; /* Petsc sstruct interface only allows 1 part */
    454:   PetscInt ordering = ex->dofs_order;
    455:   PetscInt grid_rank;
    456:   PetscInt var_type;
    457:   PetscInt nvars = ex->nvars;
    458:   PetscInt row;

    460:   PetscFunctionBegin;
    461:   PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support");
    462:   PetscCall(PetscMalloc1(7 * nvars, &entries));

    464:   PetscCall(PetscMalloc1(nvars, &values));
    465:   PetscCall(PetscMalloc1(7 * nvars * nvars, &values[0]));
    466:   for (i = 1; i < nvars; i++) values[i] = values[i - 1] + nvars * 7;

    468:   for (i = 0; i < nvars; i++) {
    469:     PetscCall(PetscArrayzero(values[i], nvars * 7 * sizeof(HYPRE_Complex)));
    470:     PetscCall(PetscHYPREScalarCast(d, values[i] + 3));
    471:   }

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

    475:   if (!ordering) {
    476:     for (i = 0; i < nrow; i++) {
    477:       grid_rank = irow[i] / nvars;
    478:       var_type  = (irow[i] % nvars);

    480:       row      = ex->gindices[grid_rank] - ex->rstart;
    481:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
    482:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
    483:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
    484:       PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]);
    485:     }
    486:   } else {
    487:     for (i = 0; i < nrow; i++) {
    488:       var_type  = irow[i] / (ex->gnxgnygnz);
    489:       grid_rank = irow[i] - var_type * (ex->gnxgnygnz);

    491:       row      = ex->gindices[grid_rank] - ex->rstart;
    492:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
    493:       index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny));
    494:       index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny)));
    495:       PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]);
    496:     }
    497:   }
    498:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
    499:   PetscCall(PetscFree(values[0]));
    500:   PetscCall(PetscFree(values));
    501:   PetscCall(PetscFree(entries));
    502:   PetscFunctionReturn(PETSC_SUCCESS);
    503: }

    505: static PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
    506: {
    507:   Mat_HYPRESStruct *ex    = (Mat_HYPRESStruct *)mat->data;
    508:   PetscInt          nvars = ex->nvars;
    509:   PetscInt          size;
    510:   PetscInt          part = 0; /* only one part */

    512:   PetscFunctionBegin;
    513:   size = (ex->hbox.imax[0] - ex->hbox.imin[0] + 1) * (ex->hbox.imax[1] - ex->hbox.imin[1] + 1) * (ex->hbox.imax[2] - ex->hbox.imin[2] + 1);
    514:   {
    515:     HYPRE_Int      i, *entries, iupper[3], ilower[3];
    516:     HYPRE_Complex *values;

    518:     for (i = 0; i < 3; i++) {
    519:       ilower[i] = ex->hbox.imin[i];
    520:       iupper[i] = ex->hbox.imax[i];
    521:     }

    523:     PetscCall(PetscMalloc2(nvars * 7, &entries, nvars * 7 * size, &values));
    524:     for (i = 0; i < nvars * 7; i++) entries[i] = i;
    525:     PetscCall(PetscArrayzero(values, nvars * 7 * size));

    527:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructMatrixSetBoxValues, ex->ss_mat, part, ilower, iupper, i, nvars * 7, entries, values);
    528:     PetscCall(PetscFree2(entries, values));
    529:   }
    530:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
    531:   PetscFunctionReturn(PETSC_SUCCESS);
    532: }

    534: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
    535: {
    536:   Mat_HYPRESStruct      *ex = (Mat_HYPRESStruct *)mat->data;
    537:   PetscInt               dim, dof, sw[3], nx, ny, nz;
    538:   PetscInt               ilower[3], iupper[3], ssize, i;
    539:   DMBoundaryType         px, py, pz;
    540:   DMDAStencilType        st;
    541:   PetscInt               nparts = 1; /* assuming only one part */
    542:   PetscInt               part   = 0;
    543:   ISLocalToGlobalMapping ltog;
    544:   DM                     da;

    546:   PetscFunctionBegin;
    547:   PetscCall(MatGetDM(mat, (DM *)&da));
    548:   ex->da = da;
    549:   PetscCall(PetscObjectReference((PetscObject)da));

    551:   PetscCall(DMDAGetInfo(ex->da, &dim, 0, 0, 0, 0, 0, 0, &dof, &sw[0], &px, &py, &pz, &st));
    552:   PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
    553:   iupper[0] += ilower[0] - 1;
    554:   iupper[1] += ilower[1] - 1;
    555:   iupper[2] += ilower[2] - 1;
    556:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
    557:   ex->hbox.imin[0] = ilower[0];
    558:   ex->hbox.imin[1] = ilower[1];
    559:   ex->hbox.imin[2] = ilower[2];
    560:   ex->hbox.imax[0] = iupper[0];
    561:   ex->hbox.imax[1] = iupper[1];
    562:   ex->hbox.imax[2] = iupper[2];

    564:   ex->dofs_order = 0;

    566:   /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
    567:   ex->nvars = (int)dof;

    569:   /* create the hypre grid object and set its information */
    570:   PetscCheck(!px && !py && !pz, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
    571:   PetscCallExternal(HYPRE_SStructGridCreate, ex->hcomm, dim, nparts, &ex->ss_grid);
    572:   PetscCallExternal(HYPRE_SStructGridSetExtents, ex->ss_grid, part, ex->hbox.imin, ex->hbox.imax);
    573:   {
    574:     HYPRE_SStructVariable *vartypes;
    575:     PetscCall(PetscMalloc1(ex->nvars, &vartypes));
    576:     for (i = 0; i < ex->nvars; i++) vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL;
    577:     PetscCallExternal(HYPRE_SStructGridSetVariables, ex->ss_grid, part, ex->nvars, vartypes);
    578:     PetscCall(PetscFree(vartypes));
    579:   }
    580:   PetscCallExternal(HYPRE_SStructGridAssemble, ex->ss_grid);

    582:   sw[1] = sw[0];
    583:   sw[2] = sw[1];
    584:   /* PetscCallExternal(HYPRE_SStructGridSetNumGhost,ex->ss_grid,sw); */

    586:   /* create the hypre stencil object and set its information */
    587:   PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils");
    588:   PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils");

    590:   if (dim == 1) {
    591:     HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}};
    592:     PetscInt  j, cnt;

    594:     ssize = 3 * (ex->nvars);
    595:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
    596:     cnt = 0;
    597:     for (i = 0; i < (ex->nvars); i++) {
    598:       for (j = 0; j < 3; j++) {
    599:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
    600:         cnt++;
    601:       }
    602:     }
    603:   } else if (dim == 2) {
    604:     HYPRE_Int offsets[5][2] = {
    605:       {0,  -1},
    606:       {-1, 0 },
    607:       {0,  0 },
    608:       {1,  0 },
    609:       {0,  1 }
    610:     };
    611:     PetscInt j, cnt;

    613:     ssize = 5 * (ex->nvars);
    614:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
    615:     cnt = 0;
    616:     for (i = 0; i < (ex->nvars); i++) {
    617:       for (j = 0; j < 5; j++) {
    618:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
    619:         cnt++;
    620:       }
    621:     }
    622:   } else if (dim == 3) {
    623:     HYPRE_Int offsets[7][3] = {
    624:       {0,  0,  -1},
    625:       {0,  -1, 0 },
    626:       {-1, 0,  0 },
    627:       {0,  0,  0 },
    628:       {1,  0,  0 },
    629:       {0,  1,  0 },
    630:       {0,  0,  1 }
    631:     };
    632:     PetscInt j, cnt;

    634:     ssize = 7 * (ex->nvars);
    635:     PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil);
    636:     cnt = 0;
    637:     for (i = 0; i < (ex->nvars); i++) {
    638:       for (j = 0; j < 7; j++) {
    639:         PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i);
    640:         cnt++;
    641:       }
    642:     }
    643:   }

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

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

    653:   /* create the HYPRE sstruct vectors for rhs and solution */
    654:   PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_b);
    655:   PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_x);
    656:   PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_b);
    657:   PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_x);
    658:   PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_b);
    659:   PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_x);

    661:   /* create the hypre matrix object and set its information */
    662:   PetscCallExternal(HYPRE_SStructMatrixCreate, ex->hcomm, ex->ss_graph, &ex->ss_mat);
    663:   PetscCallExternal(HYPRE_SStructGridDestroy, ex->ss_grid);
    664:   PetscCallExternal(HYPRE_SStructStencilDestroy, ex->ss_stencil);
    665:   if (ex->needsinitialization) {
    666:     PetscCallExternal(HYPRE_SStructMatrixInitialize, ex->ss_mat);
    667:     ex->needsinitialization = PETSC_FALSE;
    668:   }

    670:   /* set the global and local sizes of the matrix */
    671:   PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz));
    672:   PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE));
    673:   PetscCall(PetscLayoutSetBlockSize(mat->rmap, dof));
    674:   PetscCall(PetscLayoutSetBlockSize(mat->cmap, dof));
    675:   PetscCall(PetscLayoutSetUp(mat->rmap));
    676:   PetscCall(PetscLayoutSetUp(mat->cmap));
    677:   mat->preallocated = PETSC_TRUE;

    679:   if (dim == 3) {
    680:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
    681:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
    682:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

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

    687:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
    688:   PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL));
    689:   PetscCall(DMGetLocalToGlobalMapping(ex->da, &ltog));
    690:   PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices));
    691:   PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, &ex->gnxgnygnz));

    693:   ex->gnxgny *= ex->gnx;
    694:   ex->gnxgnygnz *= ex->gnxgny;

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

    698:   ex->nxny   = ex->nx * ex->ny;
    699:   ex->nxnynz = ex->nz * ex->nxny;
    700:   PetscFunctionReturn(PETSC_SUCCESS);
    701: }

    703: static PetscErrorCode MatMult_HYPRESStruct(Mat A, Vec x, Vec y)
    704: {
    705:   const PetscScalar *xx;
    706:   PetscScalar       *yy;
    707:   HYPRE_Int          hlower[3], hupper[3];
    708:   PetscInt           ilower[3], iupper[3];
    709:   Mat_HYPRESStruct  *mx       = (Mat_HYPRESStruct *)A->data;
    710:   PetscInt           ordering = mx->dofs_order;
    711:   PetscInt           nvars    = mx->nvars;
    712:   PetscInt           part     = 0;
    713:   PetscInt           size;
    714:   PetscInt           i;

    716:   PetscFunctionBegin;
    717:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));

    719:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
    720:   iupper[0] += ilower[0] - 1;
    721:   iupper[1] += ilower[1] - 1;
    722:   iupper[2] += ilower[2] - 1;
    723:   hlower[0] = (HYPRE_Int)ilower[0];
    724:   hlower[1] = (HYPRE_Int)ilower[1];
    725:   hlower[2] = (HYPRE_Int)ilower[2];
    726:   hupper[0] = (HYPRE_Int)iupper[0];
    727:   hupper[1] = (HYPRE_Int)iupper[1];
    728:   hupper[2] = (HYPRE_Int)iupper[2];

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

    733:   /* copy x values over to hypre for variable ordering */
    734:   if (ordering) {
    735:     PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
    736:     PetscCall(VecGetArrayRead(x, &xx));
    737:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i)));
    738:     PetscCall(VecRestoreArrayRead(x, &xx));
    739:     PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
    740:     PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);

    742:     /* copy solution values back to PETSc */
    743:     PetscCall(VecGetArray(y, &yy));
    744:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i)));
    745:     PetscCall(VecRestoreArray(y, &yy));
    746:   } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
    747:     PetscScalar *z;
    748:     PetscInt     j, k;

    750:     PetscCall(PetscMalloc1(nvars * size, &z));
    751:     PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
    752:     PetscCall(VecGetArrayRead(x, &xx));

    754:     /* transform nodal to hypre's variable ordering for sys_pfmg */
    755:     for (i = 0; i < size; i++) {
    756:       k = i * nvars;
    757:       for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
    758:     }
    759:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
    760:     PetscCall(VecRestoreArrayRead(x, &xx));
    761:     PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
    762:     PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);

    764:     /* copy solution values back to PETSc */
    765:     PetscCall(VecGetArray(y, &yy));
    766:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
    767:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
    768:     for (i = 0; i < size; i++) {
    769:       k = i * nvars;
    770:       for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
    771:     }
    772:     PetscCall(VecRestoreArray(y, &yy));
    773:     PetscCall(PetscFree(z));
    774:   }
    775:   PetscFunctionReturn(PETSC_SUCCESS);
    776: }

    778: static PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode)
    779: {
    780:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data;

    782:   PetscFunctionBegin;
    783:   PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat);
    784:   PetscFunctionReturn(PETSC_SUCCESS);
    785: }

    787: static PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
    788: {
    789:   PetscFunctionBegin;
    790:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
    791:   PetscFunctionReturn(PETSC_SUCCESS);
    792: }

    794: static PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
    795: {
    796:   Mat_HYPRESStruct      *ex = (Mat_HYPRESStruct *)mat->data;
    797:   ISLocalToGlobalMapping ltog;

    799:   PetscFunctionBegin;
    800:   PetscCall(DMGetLocalToGlobalMapping(ex->da, &ltog));
    801:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **)&ex->gindices));
    802:   PetscCallExternal(HYPRE_SStructGraphDestroy, ex->ss_graph);
    803:   PetscCallExternal(HYPRE_SStructMatrixDestroy, ex->ss_mat);
    804:   PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_x);
    805:   PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_b);
    806:   PetscCall(PetscObjectDereference((PetscObject)ex->da));
    807:   PetscCallMPI(MPI_Comm_free(&ex->hcomm));
    808:   PetscCall(PetscFree(ex));
    809:   PetscFunctionReturn(PETSC_SUCCESS);
    810: }

    812: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
    813: {
    814:   Mat_HYPRESStruct *ex;

    816:   PetscFunctionBegin;
    817:   PetscCall(PetscNew(&ex));
    818:   B->data      = (void *)ex;
    819:   B->rmap->bs  = 1;
    820:   B->assembled = PETSC_FALSE;

    822:   B->insertmode = NOT_SET_VALUES;

    824:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
    825:   B->ops->mult        = MatMult_HYPRESStruct;
    826:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
    827:   B->ops->destroy     = MatDestroy_HYPRESStruct;
    828:   B->ops->setup       = MatSetUp_HYPRESStruct;

    830:   ex->needsinitialization = PETSC_TRUE;

    832:   PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &ex->hcomm));
    833:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT));
    834:   PetscHYPREInitialize();
    835:   PetscFunctionReturn(PETSC_SUCCESS);
    836: }