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 #undef __FUNCT__ 11 #define __FUNCT__ "DMCreateMatrix_Moab" 12 PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J) 13 { 14 PetscErrorCode ierr; 15 PetscInt innz = 0, ionz = 0, nlsiz; 16 DM_Moab *dmmoab = (DM_Moab*)dm->data; 17 PetscInt *nnz = 0, *onz = 0; 18 char *tmp = 0; 19 Mat A; 20 MatType mtype; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 24 PetscValidPointer(J, 3); 25 26 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 27 mtype = dm->mattype; 28 ierr = PetscStrstr(mtype, "baij", &tmp);CHKERRQ(ierr); 29 nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields); 30 31 /* allocate the nnz, onz arrays based on block size and local nodes */ 32 ierr = PetscCalloc2(nlsiz, &nnz, nlsiz, &onz);CHKERRQ(ierr); 33 34 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 35 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, (tmp ? PETSC_TRUE : PETSC_FALSE));CHKERRQ(ierr); 36 37 /* create the Matrix and set its type as specified by user */ 38 ierr = MatCreate((((PetscObject)dm)->comm), &A);CHKERRQ(ierr); 39 ierr = MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 40 ierr = MatSetType(A, mtype);CHKERRQ(ierr); 41 ierr = MatSetBlockSize(A, dmmoab->bs);CHKERRQ(ierr); 42 ierr = MatSetDM(A, dm);CHKERRQ(ierr); /* set DM reference */ 43 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 44 45 if (!dmmoab->ltog_map) SETERRQ((((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 46 ierr = MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);CHKERRQ(ierr); 47 48 /* set preallocation based on different supported Mat types */ 49 ierr = MatSeqAIJSetPreallocation(A, innz, nnz);CHKERRQ(ierr); 50 ierr = MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz);CHKERRQ(ierr); 51 ierr = MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 52 ierr = MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 53 54 /* clean up temporary memory */ 55 ierr = PetscFree2(nnz, onz);CHKERRQ(ierr); 56 57 /* set up internal matrix data-structures */ 58 ierr = MatSetUp(A);CHKERRQ(ierr); 59 60 /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */ 61 62 *J = A; 63 PetscFunctionReturn(0); 64 } 65 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 69 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt* innz, PetscInt* nnz, PetscInt* ionz, PetscInt* onz, PetscBool isbaij) 70 { 71 PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0; 72 PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz; 73 DM_Moab *dmmoab = (DM_Moab*)dm->data; 74 moab::Range found; 75 std::vector<moab::EntityHandle> adjs, storage; 76 PetscBool isinterlaced; 77 moab::EntityHandle vtx; 78 moab::ErrorCode merr; 79 80 PetscFunctionBegin; 81 bs = dmmoab->bs; 82 nloc = dmmoab->nloc; 83 nfields = dmmoab->numFields; 84 isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE); 85 nlsiz = (isinterlaced ? nloc : nloc * nfields); 86 87 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 88 for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) { 89 90 vtx = *iter; 91 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 92 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 93 non-zero pattern accordingly. */ 94 adjs.clear(); 95 if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) { 96 merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs); MBERRNM(merr); 97 } 98 else { 99 merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION); MBERRNM(merr); 100 } 101 102 /* reset counters */ 103 n_nnz = n_onz = 0; 104 found.clear(); 105 106 /* loop over vertices and update the number of connectivity */ 107 for (unsigned jter = 0; jter < adjs.size(); ++jter) { 108 109 /* Get connectivity information in canonical ordering for the local element */ 110 const moab::EntityHandle *connect; 111 std::vector<moab::EntityHandle> cconnect; 112 merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage); MBERRNM(merr); 113 114 /* loop over each element connected to the adjacent vertex and update as needed */ 115 for (i = 0; i < vpere; ++i) { 116 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 117 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 118 if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */ 119 else n_nnz++; /* else local vertex */ 120 found.insert(connect[i]); 121 } 122 } 123 storage.clear(); 124 125 if (isbaij) { 126 nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 127 if (onz) { 128 onz[ivtx] = n_onz; /* add ghost non-owned nodes */ 129 } 130 } 131 else { /* AIJ matrices */ 132 if (!isinterlaced) { 133 for (f = 0; f < nfields; f++) { 134 nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 135 if (onz) 136 onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */ 137 } 138 } 139 else { 140 for (f = 0; f < nfields; f++) { 141 nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 142 if (onz) 143 onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */ 144 } 145 } 146 } 147 } 148 149 for (i = 0; i < nlsiz; i++) 150 nnz[i] += 1; /* self count the node */ 151 152 for (ivtx = 0; ivtx < nloc; ivtx++) { 153 if (!isbaij) { 154 for (ibs = 0; ibs < nfields; ibs++) { 155 if (dmmoab->dfill) { /* first address the diagonal block */ 156 /* just add up the ints -- easier/faster rather than branching based on "1" */ 157 for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) 158 inbsize += dmmoab->dfill[ibs * nfields + jbs]; 159 } 160 else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 161 if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize; 162 else nnz[ibs * nloc + ivtx] *= inbsize; 163 164 if (onz) { 165 if (dmmoab->ofill) { /* next address the off-diagonal block */ 166 /* just add up the ints -- easier/faster rather than branching based on "1" */ 167 for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) 168 iobsize += dmmoab->dfill[ibs * nfields + jbs]; 169 } 170 else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 171 if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize; 172 else onz[ibs * nloc + ivtx] *= iobsize; 173 } 174 } 175 } 176 else { 177 /* check if we got overzealous in our nnz and onz computations */ 178 nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]); 179 if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]); 180 } 181 } 182 /* update innz and ionz based on local maxima */ 183 if (innz || ionz) { 184 if (innz) *innz = 0; 185 if (ionz) *ionz = 0; 186 if (isbaij) { 187 for (i = 0; i < nlsiz; i++) { 188 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 189 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 190 } 191 } 192 else { 193 for (i = 0; i < nlsiz * nfields; i++) { 194 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 195 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 196 } 197 } 198 } 199 PetscFunctionReturn(0); 200 } 201 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "DMMoabSetBlockFills_Private" 205 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 206 { 207 PetscErrorCode ierr; 208 PetscInt i, j, *ifill; 209 210 PetscFunctionBegin; 211 if (!fill) PetscFunctionReturn(0); 212 ierr = PetscMalloc1(w * w, &ifill);CHKERRQ(ierr); 213 for (i = 0; i < w; i++) { 214 for (j = 0; j < w; j++) 215 ifill[i * w + j] = fill[i * w + j]; 216 } 217 218 *rfill = ifill; 219 PetscFunctionReturn(0); 220 } 221 222 223 /*@ 224 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 225 of the matrix returned by DMCreateMatrix(). 226 227 Logically Collective on DMDA 228 229 Input Parameter: 230 + dm - the DMMoab object 231 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 232 - ofill - the fill pattern in the off-diagonal blocks 233 234 Level: developer 235 236 Notes: This only makes sense when you are doing multicomponent problems but using the 237 MPIAIJ matrix format 238 239 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 240 representing coupling and 0 entries for missing coupling. For example 241 $ dfill[9] = {1, 0, 0, 242 $ 1, 1, 0, 243 $ 0, 1, 1} 244 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 245 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 246 diagonal block). 247 248 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 249 can be represented in the dfill, ofill format 250 251 Contributed by Glenn Hammond 252 253 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 254 255 @*/ 256 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 257 { 258 DM_Moab *dmmoab = (DM_Moab*)dm->data; 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 263 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);CHKERRQ(ierr); 264 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267