Actual source code: dmmbmat.cxx

      1: #include <petsc/private/dmmbimpl.h>
      2: #include <petsc/private/vecimpl.h>

      4: #include <petscdmmoab.h>
      5: #include <MBTagConventions.hpp>
      6: #include <moab/NestedRefine.hpp>

      8: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);

     10: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J)
     11: {
     12:   PetscInt  innz = 0, ionz = 0, nlsiz;
     13:   DM_Moab  *dmmoab = (DM_Moab *)dm->data;
     14:   PetscInt *nnz = 0, *onz = 0;
     15:   char     *tmp = 0;
     16:   Mat       A;
     17:   MatType   mtype;

     19:   PetscFunctionBegin;
     21:   PetscAssertPointer(J, 3);

     23:   /* next, need to allocate the non-zero arrays to enable pre-allocation */
     24:   mtype = dm->mattype;
     25:   PetscCall(PetscStrstr(mtype, MATBAIJ, &tmp));
     26:   nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields);

     28:   /* allocate the nnz, onz arrays based on block size and local nodes */
     29:   PetscCall(PetscCalloc2(nlsiz, &nnz, nlsiz, &onz));

     31:   /* compute the nonzero pattern based on MOAB connectivity data for local elements */
     32:   PetscCall(DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, tmp ? PETSC_TRUE : PETSC_FALSE));

     34:   /* create the Matrix and set its type as specified by user */
     35:   PetscCall(MatCreate(((PetscObject)dm)->comm, &A));
     36:   PetscCall(MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE));
     37:   PetscCall(MatSetType(A, mtype));
     38:   PetscCall(MatSetBlockSize(A, dmmoab->bs));
     39:   PetscCall(MatSetDM(A, dm)); /* set DM reference */
     40:   PetscCall(MatSetFromOptions(A));

     42:   PetscCheck(dmmoab->ltog_map, ((PetscObject)dm)->comm, PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first.");
     43:   PetscCall(MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map));

     45:   /* set preallocation based on different supported Mat types */
     46:   PetscCall(MatSeqAIJSetPreallocation(A, innz, nnz));
     47:   PetscCall(MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz));
     48:   PetscCall(MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz));
     49:   PetscCall(MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz));

     51:   /* clean up temporary memory */
     52:   PetscCall(PetscFree2(nnz, onz));

     54:   /* set up internal matrix data-structures */
     55:   PetscCall(MatSetUp(A));

     57:   /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */

     59:   *J = A;
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt *innz, PetscInt *nnz, PetscInt *ionz, PetscInt *onz, PetscBool isbaij)
     64: {
     65:   PetscInt                        i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0;
     66:   PetscInt                        ibs, jbs, inbsize, iobsize, nfields, nlsiz;
     67:   DM_Moab                        *dmmoab = (DM_Moab *)dm->data;
     68:   moab::Range                     found;
     69:   std::vector<moab::EntityHandle> adjs, storage;
     70:   PetscBool                       isinterlaced;
     71:   moab::EntityHandle              vtx;
     72:   moab::ErrorCode                 merr;

     74:   PetscFunctionBegin;
     75:   bs           = dmmoab->bs;
     76:   nloc         = dmmoab->nloc;
     77:   nfields      = dmmoab->numFields;
     78:   isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE);
     79:   nlsiz        = (isinterlaced ? nloc : nloc * nfields);

     81:   /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */
     82:   for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) {
     83:     vtx = *iter;
     84:     /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
     85:        to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
     86:        non-zero pattern accordingly. */
     87:     adjs.clear();
     88:     if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) {
     89:       merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs);
     90:       MBERRNM(merr);
     91:     } else {
     92:       merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION);
     93:       MBERRNM(merr);
     94:     }

     96:     /* reset counters */
     97:     n_nnz = n_onz = 0;
     98:     found.clear();

    100:     /* loop over vertices and update the number of connectivity */
    101:     for (unsigned jter = 0; jter < adjs.size(); ++jter) {
    102:       /* Get connectivity information in canonical ordering for the local element */
    103:       const moab::EntityHandle       *connect;
    104:       std::vector<moab::EntityHandle> cconnect;
    105:       merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage);
    106:       MBERRNM(merr);

    108:       /* loop over each element connected to the adjacent vertex and update as needed */
    109:       for (i = 0; i < vpere; ++i) {
    110:         /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
    111:         if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */
    112:         if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++;   /* update out-of-proc onz */
    113:         else n_nnz++;                                                             /* else local vertex */
    114:         found.insert(connect[i]);
    115:       }
    116:     }
    117:     storage.clear();

    119:     if (isbaij) {
    120:       nnz[ivtx] = n_nnz;          /* leave out self to avoid repeats -> node shared by multiple elements */
    121:       if (onz) onz[ivtx] = n_onz; /* add ghost non-owned nodes */
    122:     } else {                      /* AIJ matrices */
    123:       if (!isinterlaced) {
    124:         for (f = 0; f < nfields; f++) {
    125:           nnz[f * nloc + ivtx] = n_nnz;          /* leave out self to avoid repeats -> node shared by multiple elements */
    126:           if (onz) onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */
    127:         }
    128:       } else {
    129:         for (f = 0; f < nfields; f++) {
    130:           nnz[nfields * ivtx + f] = n_nnz;          /* leave out self to avoid repeats -> node shared by multiple elements */
    131:           if (onz) onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */
    132:         }
    133:       }
    134:     }
    135:   }

    137:   for (i = 0; i < nlsiz; i++) nnz[i] += 1; /* self count the node */

    139:   for (ivtx = 0; ivtx < nloc; ivtx++) {
    140:     if (!isbaij) {
    141:       for (ibs = 0; ibs < nfields; ibs++) {
    142:         if (dmmoab->dfill) { /* first address the diagonal block */
    143:           /* just add up the ints -- easier/faster rather than branching based on "1" */
    144:           for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) inbsize += dmmoab->dfill[ibs * nfields + jbs];
    145:         } else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
    146:         if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize;
    147:         else nnz[ibs * nloc + ivtx] *= inbsize;

    149:         if (onz) {
    150:           if (dmmoab->ofill) { /* next address the off-diagonal block */
    151:             /* just add up the ints -- easier/faster rather than branching based on "1" */
    152:             for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) iobsize += dmmoab->dfill[ibs * nfields + jbs];
    153:           } else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
    154:           if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize;
    155:           else onz[ibs * nloc + ivtx] *= iobsize;
    156:         }
    157:       }
    158:     } else {
    159:       /* check if we got overzealous in our nnz and onz computations */
    160:       nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]);
    161:       if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]);
    162:     }
    163:   }
    164:   /* update innz and ionz based on local maxima */
    165:   if (innz || ionz) {
    166:     if (innz) *innz = 0;
    167:     if (ionz) *ionz = 0;
    168:     for (i = 0; i < nlsiz; i++) {
    169:       if (innz && nnz[i] > *innz) *innz = nnz[i];
    170:       if (ionz && onz && onz[i] > *ionz) *ionz = onz[i];
    171:     }
    172:   }
    173:   PetscFunctionReturn(PETSC_SUCCESS);
    174: }

    176: static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill)
    177: {
    178:   PetscInt i, j, *ifill;

    180:   PetscFunctionBegin;
    181:   if (!fill) PetscFunctionReturn(PETSC_SUCCESS);
    182:   PetscCall(PetscMalloc1(w * w, &ifill));
    183:   for (i = 0; i < w; i++) {
    184:     for (j = 0; j < w; j++) ifill[i * w + j] = fill[i * w + j];
    185:   }

    187:   *rfill = ifill;
    188:   PetscFunctionReturn(PETSC_SUCCESS);
    189: }

    191: /*@C
    192:   DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem
    193:   of the matrix returned by `DMCreateMatrix()`.

    195:   Logically Collective

    197:   Input Parameters:
    198: + dm    - the `DMMOAB` object
    199: . dfill - the fill pattern in the diagonal block (may be `NULL`, means use dense block)
    200: - ofill - the fill pattern in the off-diagonal blocks

    202:   Level: developer

    204:   Notes:
    205:   This only makes sense when you are doing multicomponent problems but using the
    206:   MPIAIJ matrix format

    208:   The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
    209:   representing coupling and 0 entries for missing coupling. For example
    210: .vb
    211:              dfill[9] = {1, 0, 0,
    212:                          1, 1, 0,
    213:                          0, 1, 1}
    214: .ve
    215:   means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
    216:   itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
    217:   diagonal block).

    219:   `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then
    220:   can be represented in the dfill, ofill format

    222:   Contributed by\:
    223:   Glenn Hammond

    225: .seealso: `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`
    226: @*/
    227: PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
    228: {
    229:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    231:   PetscFunctionBegin;
    233:   PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill));
    234:   PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill));
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }