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 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; 19 20 PetscFunctionBegin; 21 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 22 PetscValidPointer(J, 3); 23 24 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 25 mtype = dm->mattype; 26 ierr = PetscStrstr(mtype, "baij", &tmp);CHKERRQ(ierr); 27 nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields); 28 29 /* allocate the nnz, onz arrays based on block size and local nodes */ 30 ierr = PetscCalloc2(nlsiz, &nnz, nlsiz, &onz);CHKERRQ(ierr); 31 32 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 33 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, (tmp ? PETSC_TRUE : PETSC_FALSE));CHKERRQ(ierr); 34 35 /* create the Matrix and set its type as specified by user */ 36 ierr = MatCreate((((PetscObject)dm)->comm), &A);CHKERRQ(ierr); 37 ierr = MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 38 ierr = MatSetType(A, mtype);CHKERRQ(ierr); 39 ierr = MatSetBlockSize(A, dmmoab->bs);CHKERRQ(ierr); 40 ierr = MatSetDM(A, dm);CHKERRQ(ierr); /* set DM reference */ 41 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 42 43 if (!dmmoab->ltog_map) SETERRQ((((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 44 ierr = MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);CHKERRQ(ierr); 45 46 /* set preallocation based on different supported Mat types */ 47 ierr = MatSeqAIJSetPreallocation(A, innz, nnz);CHKERRQ(ierr); 48 ierr = MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz);CHKERRQ(ierr); 49 ierr = MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 50 ierr = MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 51 52 /* clean up temporary memory */ 53 ierr = PetscFree2(nnz, onz);CHKERRQ(ierr); 54 55 /* set up internal matrix data-structures */ 56 ierr = MatSetUp(A);CHKERRQ(ierr); 57 58 /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */ 59 60 *J = A; 61 PetscFunctionReturn(0); 62 } 63 64 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; 75 76 PetscFunctionBegin; 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); 82 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++) { 85 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 } 97 98 /* reset counters */ 99 n_nnz = n_onz = 0; 100 found.clear(); 101 102 /* loop over vertices and update the number of connectivity */ 103 for (unsigned jter = 0; jter < adjs.size(); ++jter) { 104 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); 109 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(); 120 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 } 144 145 for (i = 0; i < nlsiz; i++) 146 nnz[i] += 1; /* self count the node */ 147 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; 159 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 PetscFunctionReturn(0); 188 } 189 190 191 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 192 { 193 PetscErrorCode ierr; 194 PetscInt i, j, *ifill; 195 196 PetscFunctionBegin; 197 if (!fill) PetscFunctionReturn(0); 198 ierr = PetscMalloc1(w * w, &ifill);CHKERRQ(ierr); 199 for (i = 0; i < w; i++) { 200 for (j = 0; j < w; j++) 201 ifill[i * w + j] = fill[i * w + j]; 202 } 203 204 *rfill = ifill; 205 PetscFunctionReturn(0); 206 } 207 208 209 /*@ 210 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 211 of the matrix returned by DMCreateMatrix(). 212 213 Logically Collective on DMDA 214 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 219 220 Level: developer 221 222 Notes: This only makes sense when you are doing multicomponent problems but using the 223 MPIAIJ matrix format 224 225 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 226 representing coupling and 0 entries for missing coupling. For example 227 $ dfill[9] = {1, 0, 0, 228 $ 1, 1, 0, 229 $ 0, 1, 1} 230 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 231 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 232 diagonal block). 233 234 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 235 can be represented in the dfill, ofill format 236 237 Contributed by Glenn Hammond 238 239 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 240 241 @*/ 242 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 243 { 244 DM_Moab *dmmoab = (DM_Moab*)dm->data; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 249 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);CHKERRQ(ierr); 250 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253