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