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