1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/ 2 #include <petsc/private/vecimpl.h> 3 4 #include <petscdmmoab.h> 5 #include <MBTagConventions.hpp> 6 #include <moab/NestedRefine.hpp> 7 8 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool); 9 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; 18 19 PetscFunctionBegin; 20 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 21 PetscAssertPointer(J, 3); 22 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); 27 28 /* allocate the nnz, onz arrays based on block size and local nodes */ 29 PetscCall(PetscCalloc2(nlsiz, &nnz, nlsiz, &onz)); 30 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))); 33 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)); 41 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)); 44 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)); 50 51 /* clean up temporary memory */ 52 PetscCall(PetscFree2(nnz, onz)); 53 54 /* set up internal matrix data-structures */ 55 PetscCall(MatSetUp(A)); 56 57 /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */ 58 59 *J = A; 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 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; 73 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); 80 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 } 95 96 /* reset counters */ 97 n_nnz = n_onz = 0; 98 found.clear(); 99 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); 107 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(); 118 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 } 136 137 for (i = 0; i < nlsiz; i++) nnz[i] += 1; /* self count the node */ 138 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; 148 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 } 175 176 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 177 { 178 PetscInt i, j, *ifill; 179 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 } 186 187 *rfill = ifill; 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 /*@C 192 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 193 of the matrix returned by `DMCreateMatrix()`. 194 195 Logically Collective 196 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 201 202 Level: developer 203 204 Notes: 205 This only makes sense when you are doing multicomponent problems but using the 206 MPIAIJ matrix format 207 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). 218 219 `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then 220 can be represented in the dfill, ofill format 221 222 Contributed by Glenn Hammond 223 224 .seealso: `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()` 225 @*/ 226 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 227 { 228 DM_Moab *dmmoab = (DM_Moab *)dm->data; 229 230 PetscFunctionBegin; 231 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 232 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill)); 233 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill)); 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236