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 *J = A; 61 PetscFunctionReturn(0); 62 } 63 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 67 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt* innz, PetscInt* nnz, PetscInt* ionz, PetscInt* onz, PetscBool isbaij) 68 { 69 PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0; 70 PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz; 71 DM_Moab *dmmoab = (DM_Moab*)dm->data; 72 const moab::EntityHandle *connect; 73 moab::Range adjs, found; 74 std::vector<moab::EntityHandle> storage; 75 PetscBool isinterlaced; 76 moab::EntityHandle vtx; 77 moab::ErrorCode merr; 78 79 PetscFunctionBegin; 80 bs = dmmoab->bs; 81 nloc = dmmoab->nloc; 82 nfields = dmmoab->numFields; 83 isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE); 84 nlsiz = (isinterlaced ? nloc : nloc * nfields); 85 86 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 87 for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) { 88 89 vtx = *iter; 90 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 91 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 92 non-zero pattern accordingly. */ 93 adjs.clear(); 94 merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::INTERSECT); MBERRNM(merr); 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 (moab::Range::const_iterator jter = adjs.begin(); jter != adjs.end(); jter++) { 102 103 const moab::EntityHandle jhandle = *jter; 104 105 /* Get connectivity information in canonical ordering for the local element */ 106 merr = dmmoab->mbiface->get_connectivity(jhandle, connect, vpere, false, &storage); 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) { 122 onz[ivtx] = n_onz; /* add ghost non-owned nodes */ 123 } 124 } 125 else { /* AIJ matrices */ 126 if (!isinterlaced) { 127 for (f = 0; f < nfields; f++) { 128 nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 129 if (onz) 130 onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */ 131 } 132 } 133 else { 134 for (f = 0; f < nfields; f++) { 135 nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 136 if (onz) 137 onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */ 138 } 139 } 140 } 141 } 142 143 for (i = 0; i < nlsiz; i++) 144 nnz[i] += 1; /* self count the node */ 145 146 for (ivtx = 0; ivtx < nloc; ivtx++) { 147 if (!isbaij) { 148 for (ibs = 0; ibs < nfields; ibs++) { 149 if (dmmoab->dfill) { /* first address the diagonal block */ 150 /* just add up the ints -- easier/faster rather than branching based on "1" */ 151 for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) 152 inbsize += dmmoab->dfill[ibs * nfields + jbs]; 153 } 154 else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 155 if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize; 156 else nnz[ibs * nloc + ivtx] *= inbsize; 157 158 if (onz) { 159 if (dmmoab->ofill) { /* next address the off-diagonal block */ 160 /* just add up the ints -- easier/faster rather than branching based on "1" */ 161 for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) 162 iobsize += dmmoab->dfill[ibs * nfields + jbs]; 163 } 164 else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 165 if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize; 166 else onz[ibs * nloc + ivtx] *= iobsize; 167 } 168 } 169 } 170 else { 171 /* check if we got overzealous in our nnz and onz computations */ 172 nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]); 173 if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]); 174 } 175 } 176 /* update innz and ionz based on local maxima */ 177 if (innz || ionz) { 178 if (innz) *innz = 0; 179 if (ionz) *ionz = 0; 180 if (isbaij) { 181 for (i = 0; i < nlsiz; i++) { 182 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 183 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 184 } 185 } 186 else { 187 for (i = 0; i < nlsiz * nfields; i++) { 188 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 189 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 190 } 191 } 192 } 193 PetscFunctionReturn(0); 194 } 195 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "DMMoabSetBlockFills_Private" 199 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 200 { 201 PetscErrorCode ierr; 202 PetscInt i, j, *ifill; 203 204 PetscFunctionBegin; 205 if (!fill) PetscFunctionReturn(0); 206 ierr = PetscMalloc1(w * w, &ifill);CHKERRQ(ierr); 207 for (i = 0; i < w; i++) { 208 for (j = 0; j < w; j++) 209 ifill[i * w + j] = fill[i * w + j]; 210 } 211 212 *rfill = ifill; 213 PetscFunctionReturn(0); 214 } 215 216 217 /*@ 218 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 219 of the matrix returned by DMCreateMatrix(). 220 221 Logically Collective on DMDA 222 223 Input Parameter: 224 + dm - the DMMoab object 225 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 226 - ofill - the fill pattern in the off-diagonal blocks 227 228 Level: developer 229 230 Notes: This only makes sense when you are doing multicomponent problems but using the 231 MPIAIJ matrix format 232 233 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 234 representing coupling and 0 entries for missing coupling. For example 235 $ dfill[9] = {1, 0, 0, 236 $ 1, 1, 0, 237 $ 0, 1, 1} 238 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 239 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 240 diagonal block). 241 242 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 243 can be represented in the dfill, ofill format 244 245 Contributed by Glenn Hammond 246 247 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 248 249 @*/ 250 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 251 { 252 DM_Moab *dmmoab = (DM_Moab*)dm->data; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 257 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);CHKERRQ(ierr); 258 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261