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