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:   PetscErrorCode  ierr;
     13:   PetscInt        innz = 0, ionz = 0, nlsiz;
     14:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;
     15:   PetscInt        *nnz = 0, *onz = 0;
     16:   char            *tmp = 0;
     17:   Mat             A;
     18:   MatType         mtype;

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

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

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

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

     43:   if (!dmmoab->ltog_map) SETERRQ((((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first.");
     44:   MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);

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

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

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

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

     60:   *J = A;
     61:   return(0);
     62: }

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

     77:   bs = dmmoab->bs;
     78:   nloc = dmmoab->nloc;
     79:   nfields = dmmoab->numFields;
     80:   isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE);
     81:   nlsiz = (isinterlaced ? nloc : nloc * nfields);

     83:   /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */
     84:   for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) {

     86:     vtx = *iter;
     87:     /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
     88:        to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
     89:        non-zero pattern accordingly. */
     90:     adjs.clear();
     91:     if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) {
     92:       merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs); MBERRNM(merr);
     93:     }
     94:     else {
     95:       merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION); MBERRNM(merr);
     96:     }

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

    102:     /* loop over vertices and update the number of connectivity */
    103:     for (unsigned jter = 0; jter < adjs.size(); ++jter) {

    105:       /* Get connectivity information in canonical ordering for the local element */
    106:       const moab::EntityHandle *connect;
    107:       std::vector<moab::EntityHandle> cconnect;
    108:       merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage); MBERRNM(merr);

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

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

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

    148:   for (ivtx = 0; ivtx < nloc; ivtx++) {
    149:     if (!isbaij) {
    150:       for (ibs = 0; ibs < nfields; ibs++) {
    151:         if (dmmoab->dfill) {  /* first address the diagonal block */
    152:           /* just add up the ints -- easier/faster rather than branching based on "1" */
    153:           for (jbs = 0, inbsize = 0; jbs < nfields; jbs++)
    154:             inbsize += dmmoab->dfill[ibs * nfields + jbs];
    155:         }
    156:         else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
    157:         if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize;
    158:         else nnz[ibs * nloc + ivtx] *= inbsize;

    160:         if (onz) {
    161:           if (dmmoab->ofill) {  /* next address the off-diagonal block */
    162:             /* just add up the ints -- easier/faster rather than branching based on "1" */
    163:             for (jbs = 0, iobsize = 0; jbs < nfields; jbs++)
    164:               iobsize += dmmoab->dfill[ibs * nfields + jbs];
    165:           }
    166:           else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
    167:           if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize;
    168:           else onz[ibs * nloc + ivtx] *= iobsize;
    169:         }
    170:       }
    171:     }
    172:     else {
    173:       /* check if we got overzealous in our nnz and onz computations */
    174:       nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]);
    175:       if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]);
    176:     }
    177:   }
    178:   /* update innz and ionz based on local maxima */
    179:   if (innz || ionz) {
    180:     if (innz) *innz = 0;
    181:     if (ionz) *ionz = 0;
    182:     for (i = 0; i < nlsiz; i++) {
    183:       if (innz && (nnz[i] > *innz)) *innz = nnz[i];
    184:       if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i];
    185:     }
    186:   }
    187:   return(0);
    188: }

    191: static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill)
    192: {
    194:   PetscInt       i, j, *ifill;

    197:   if (!fill) return(0);
    198:   PetscMalloc1(w * w, &ifill);
    199:   for (i = 0; i < w; i++) {
    200:     for (j = 0; j < w; j++)
    201:       ifill[i * w + j] = fill[i * w + j];
    202:   }

    204:   *rfill = ifill;
    205:   return(0);
    206: }

    209: /*@C
    210:     DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem
    211:     of the matrix returned by DMCreateMatrix().

    213:     Logically Collective on da

    215:     Input Parameter:
    216: +   dm - the DMMoab object
    217: .   dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
    218: -   ofill - the fill pattern in the off-diagonal blocks

    220:     Level: developer

    222:     Notes:
    223:     This only makes sense when you are doing multicomponent problems but using the
    224:        MPIAIJ matrix format

    226:            The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
    227:        representing coupling and 0 entries for missing coupling. For example
    228: $             dfill[9] = {1, 0, 0,
    229: $                         1, 1, 0,
    230: $                         0, 1, 1}
    231:        means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
    232:        itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
    233:        diagonal block).

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

    238:    Contributed by Glenn Hammond

    240: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()

    242: @*/
    243: PetscErrorCode  DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
    244: {
    245:   DM_Moab       *dmmoab = (DM_Moab*)dm->data;

    250:   DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);
    251:   DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);
    252:   return(0);
    253: }