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