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;


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

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

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

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

 42:   MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);

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

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

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

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

 58:   *J = A;
 59:   return 0;
 60: }

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

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

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

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

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

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

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

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

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

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

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

178:   if (!fill) return 0;
179:   PetscMalloc1(w * w, &ifill);
180:   for (i = 0; i < w; i++) {
181:     for (j = 0; j < w; j++) ifill[i * w + j] = fill[i * w + j];
182:   }

184:   *rfill = ifill;
185:   return 0;
186: }

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

192:     Logically Collective on da

194:     Input Parameters:
195: +   dm - the DMMoab object
196: .   dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
197: -   ofill - the fill pattern in the off-diagonal blocks

199:     Level: developer

201:     Notes:
202:     This only makes sense when you are doing multicomponent problems but using the
203:        MPIAIJ matrix format

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

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

217:    Contributed by Glenn Hammond

219: .seealso `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`

221: @*/
222: PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
223: {
224:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

227:   DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);
228:   DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);
229:   return 0;
230: }