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