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