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