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 PetscInt innz = 0, ionz = 0, nlsiz; 12 DM_Moab *dmmoab = (DM_Moab *)dm->data; 13 PetscInt *nnz = 0, *onz = 0; 14 char *tmp = 0; 15 Mat A; 16 MatType mtype; 17 18 PetscFunctionBegin; 19 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 20 PetscValidPointer(J, 3); 21 22 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 23 mtype = dm->mattype; 24 PetscCall(PetscStrstr(mtype, MATBAIJ, &tmp)); 25 nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields); 26 27 /* allocate the nnz, onz arrays based on block size and local nodes */ 28 PetscCall(PetscCalloc2(nlsiz, &nnz, nlsiz, &onz)); 29 30 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 31 PetscCall(DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, (tmp ? PETSC_TRUE : PETSC_FALSE))); 32 33 /* create the Matrix and set its type as specified by user */ 34 PetscCall(MatCreate((((PetscObject)dm)->comm), &A)); 35 PetscCall(MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE)); 36 PetscCall(MatSetType(A, mtype)); 37 PetscCall(MatSetBlockSize(A, dmmoab->bs)); 38 PetscCall(MatSetDM(A, dm)); /* set DM reference */ 39 PetscCall(MatSetFromOptions(A)); 40 41 PetscCheck(dmmoab->ltog_map, (((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 42 PetscCall(MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map)); 43 44 /* set preallocation based on different supported Mat types */ 45 PetscCall(MatSeqAIJSetPreallocation(A, innz, nnz)); 46 PetscCall(MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz)); 47 PetscCall(MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz)); 48 PetscCall(MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz)); 49 50 /* clean up temporary memory */ 51 PetscCall(PetscFree2(nnz, onz)); 52 53 /* set up internal matrix data-structures */ 54 PetscCall(MatSetUp(A)); 55 56 /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */ 57 58 *J = A; 59 PetscFunctionReturn(0); 60 } 61 62 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt *innz, PetscInt *nnz, PetscInt *ionz, PetscInt *onz, PetscBool isbaij) { 63 PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0; 64 PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz; 65 DM_Moab *dmmoab = (DM_Moab *)dm->data; 66 moab::Range found; 67 std::vector<moab::EntityHandle> adjs, storage; 68 PetscBool isinterlaced; 69 moab::EntityHandle vtx; 70 moab::ErrorCode merr; 71 72 PetscFunctionBegin; 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); 78 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 } 93 94 /* reset counters */ 95 n_nnz = n_onz = 0; 96 found.clear(); 97 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); 105 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(); 116 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 } 134 135 for (i = 0; i < nlsiz; i++) nnz[i] += 1; /* self count the node */ 136 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; 146 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 PetscFunctionReturn(0); 172 } 173 174 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) { 175 PetscInt i, j, *ifill; 176 177 PetscFunctionBegin; 178 if (!fill) PetscFunctionReturn(0); 179 PetscCall(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 } 183 184 *rfill = ifill; 185 PetscFunctionReturn(0); 186 } 187 188 /*@C 189 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 190 of the matrix returned by DMCreateMatrix(). 191 192 Logically Collective on da 193 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 198 199 Level: developer 200 201 Notes: 202 This only makes sense when you are doing multicomponent problems but using the 203 MPIAIJ matrix format 204 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). 213 214 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 215 can be represented in the dfill, ofill format 216 217 Contributed by Glenn Hammond 218 219 .seealso `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()` 220 221 @*/ 222 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) { 223 DM_Moab *dmmoab = (DM_Moab *)dm->data; 224 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 227 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill)); 228 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill)); 229 PetscFunctionReturn(0); 230 } 231