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 static PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM,Mat); 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMCreateMatrix_Moab" 12 PetscErrorCode DMCreateMatrix_Moab(DM dm,Mat *J) 13 { 14 PetscErrorCode ierr; 15 PetscInt innz,ionz,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 = PetscMalloc((nlsiz)*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 33 ierr = PetscMemzero(nnz,sizeof(PetscInt)*(nlsiz));CHKERRQ(ierr); 34 ierr = PetscMalloc(nlsiz*sizeof(PetscInt),&onz);CHKERRQ(ierr); 35 ierr = PetscMemzero(onz,sizeof(PetscInt)*nlsiz);CHKERRQ(ierr); 36 37 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 38 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 39 40 /* create the Matrix and set its type as specified by user */ 41 ierr = MatCreate(dmmoab->pcomm->comm(), &A);CHKERRQ(ierr); 42 ierr = MatSetSizes(A, dmmoab->nloc*dmmoab->numFields, dmmoab->nloc*dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 43 ierr = MatSetType(A, mtype);CHKERRQ(ierr); 44 ierr = MatSetBlockSize(A, dmmoab->bs);CHKERRQ(ierr); 45 ierr = MatSetDM(A, dm);CHKERRQ(ierr); /* set DM reference */ 46 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 47 48 if (!dmmoab->ltog_map) SETERRQ(dmmoab->pcomm->comm(), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 49 ierr = MatSetLocalToGlobalMapping(A,dmmoab->ltog_map,dmmoab->ltog_map);CHKERRQ(ierr); 50 51 /* set preallocation based on different supported Mat types */ 52 ierr = MatSeqAIJSetPreallocation(A, innz, nnz);CHKERRQ(ierr); 53 ierr = MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz);CHKERRQ(ierr); 54 ierr = MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 55 ierr = MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 56 57 /* clean up temporary memory */ 58 ierr = PetscFree(nnz);CHKERRQ(ierr); 59 ierr = PetscFree(onz);CHKERRQ(ierr); 60 61 /* set up internal matrix data-structures */ 62 ierr = MatSetUp(A);CHKERRQ(ierr); 63 64 /* set the correct NNZ pattern by setting matrix entries - make the matrix ready to use */ 65 ierr = DMMoab_MatFillMatrixEntries_Private(dm,A);CHKERRQ(ierr); 66 67 *J = A; 68 PetscFunctionReturn(0); 69 } 70 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 74 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 75 { 76 PetscInt i,f,nloc,vpere,bs,ivtx,n_nnz,n_onz; 77 PetscInt ibs,jbs,inbsize,iobsize,nfields; 78 DM_Moab *dmmoab = (DM_Moab*)dm->data; 79 const moab::EntityHandle *connect; 80 moab::Range adjs,found,allvlocal,allvghost; 81 moab::Range::iterator iter,jter; 82 std::vector<moab::EntityHandle> storage; 83 PetscBool isinterlaced; 84 moab::EntityHandle vtx; 85 moab::ErrorCode merr; 86 87 PetscFunctionBegin; 88 bs = dmmoab->bs; 89 nloc = dmmoab->nloc; 90 nfields = dmmoab->numFields; 91 isinterlaced=(isbaij || bs==nfields ? PETSC_TRUE : PETSC_FALSE); 92 93 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 94 merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,allvlocal,true);MBERRNM(merr); 95 merr = dmmoab->pcomm->filter_pstatus(allvlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&adjs);MBERRNM(merr); 96 allvghost = moab::subtract(allvlocal, adjs); 97 98 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */ 99 for(iter = dmmoab->vowned->begin(),ivtx=0; iter != dmmoab->vowned->end(); iter++,ivtx++) { 100 101 vtx = *iter; 102 adjs.clear(); 103 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 104 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 105 non-zero pattern accordingly. */ 106 merr = dmmoab->mbiface->get_adjacencies(&vtx,1,dmmoab->dim,false,adjs,moab::Interface::INTERSECT); 107 108 /* reset counters */ 109 n_nnz=n_onz=0; 110 found.clear(); 111 112 /* loop over vertices and update the number of connectivity */ 113 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 114 115 /* Get connectivity information in canonical ordering for the local element */ 116 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 117 118 /* loop over each element connected to the adjacent vertex and update as needed */ 119 for (i=0; i<vpere; ++i) { 120 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 121 if (allvghost.find(connect[i]) != allvghost.end()) n_onz++; /* update out-of-proc onz */ 122 else n_nnz++; /* else local vertex */ 123 found.insert(connect[i]); 124 } 125 } 126 127 if (isbaij) { 128 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 129 if (onz) { 130 onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 131 } 132 } 133 else { /* AIJ matrices */ 134 if (!isinterlaced) { 135 for (f=0;f<nfields;f++) { 136 nnz[f*nloc+ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 137 if (onz) 138 onz[f*nloc+ivtx]=n_onz; /* add ghost non-owned nodes */ 139 } 140 } 141 else { 142 for (f=0;f<nfields;f++) { 143 nnz[nfields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 144 if (onz) 145 onz[nfields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 146 } 147 } 148 } 149 } 150 151 for (i=0;i<nloc*(isbaij?1:nfields);i++) 152 nnz[i]+=1; /* self count the node */ 153 154 for (ivtx=0;ivtx<nloc;ivtx++) { 155 if (!isbaij) { 156 for (ibs=0; ibs<nfields; ibs++) { 157 if (dmmoab->dfill) { /* first address the diagonal block */ 158 /* just add up the ints -- easier/faster rather than branching based on "1" */ 159 for (jbs=0,inbsize=0; jbs<nfields; jbs++) 160 inbsize += dmmoab->dfill[ibs*nfields+jbs]; 161 } 162 else inbsize=nfields; /* dense coupling since user didn't specify the component fill explicitly */ 163 if (isinterlaced) nnz[ivtx*nfields+ibs]*=inbsize; 164 else nnz[ibs*nloc+ivtx]*=inbsize; 165 166 if (onz) { 167 if (dmmoab->ofill) { /* next address the off-diagonal block */ 168 /* just add up the ints -- easier/faster rather than branching based on "1" */ 169 for (jbs=0,iobsize=0; jbs<nfields; jbs++) 170 iobsize += dmmoab->dfill[ibs*nfields+jbs]; 171 } 172 else iobsize=nfields; /* dense coupling since user didn't specify the component fill explicitly */ 173 if (isinterlaced) onz[ivtx*nfields+ibs]*=iobsize; 174 else onz[ibs*nloc+ivtx]*=iobsize; 175 } 176 } 177 } 178 else { 179 /* check if we got overzealous in our nnz and onz computations */ 180 nnz[ivtx]=(nnz[ivtx]>dmmoab->nloc?dmmoab->nloc:nnz[ivtx]); 181 if (onz) onz[ivtx]=(onz[ivtx]>dmmoab->nloc?dmmoab->nloc:onz[ivtx]); 182 } 183 } 184 /* update innz and ionz based on local maxima */ 185 if (innz || ionz) { 186 if (innz) *innz=0; 187 if (ionz) *ionz=0; 188 for (i=0;i<nloc*nfields;i++) { 189 if (innz && (nnz[i]>*innz)) *innz=nnz[i]; 190 if ((ionz && onz) && (onz[i]>*ionz)) *ionz=onz[i]; 191 } 192 } 193 PetscFunctionReturn(0); 194 } 195 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "DMMoab_MatFillMatrixEntries_Private" 199 PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM dm, Mat A) 200 { 201 DM_Moab *dmmoab = (DM_Moab*)dm->data; 202 PetscInt nconn = 0,prev_nconn = 0,bs,nloc,nfields; 203 const moab::EntityHandle *connect; 204 PetscScalar *locala=NULL; 205 PetscInt *dof_indices=NULL; 206 PetscBool isinterlaced; 207 char* tmp=0; 208 MatType mtype; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 bs = dmmoab->bs; 213 nloc = dmmoab->nloc; 214 nfields = dmmoab->numFields; 215 216 /* check whether we are updating BAIJ or AIJ matrix */ 217 ierr = MatGetType(A, &mtype);CHKERRQ(ierr); 218 ierr = PetscStrstr(mtype, "baij", &tmp);CHKERRQ(ierr); 219 isinterlaced=(tmp || bs==nfields ? PETSC_TRUE : PETSC_FALSE); 220 221 /* loop over local elements */ 222 for(moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) { 223 const moab::EntityHandle ehandle = *iter; 224 225 /* Get connectivity information: */ 226 ierr = DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);CHKERRQ(ierr); 227 228 /* if we have mixed elements or arrays have not been initialized - Allocate now */ 229 if (prev_nconn != nconn) { 230 if (locala) { 231 ierr = PetscFree(locala);CHKERRQ(ierr); 232 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 233 } 234 ierr = PetscMalloc(sizeof(PetscScalar)*nconn*nconn*nfields*nfields,&locala);CHKERRQ(ierr); 235 ierr = PetscMemzero(locala,sizeof(PetscScalar)*nconn*nconn*nfields*nfields);CHKERRQ(ierr); 236 ierr = PetscMalloc(sizeof(PetscInt)*nconn*(isinterlaced?1:nfields),&dof_indices);CHKERRQ(ierr); 237 prev_nconn=nconn; 238 } 239 240 if (isinterlaced) { 241 /* get the global DOF number to appropriately set the element contribution in the RHS vector */ 242 ierr = DMMoabGetDofsBlockedLocal(dm, nconn, connect, dof_indices);CHKERRQ(ierr); 243 244 /* set the values directly into appropriate locations. Can alternately use VecSetValues */ 245 ierr = MatSetValuesBlockedLocal(A, nconn, dof_indices, nconn, dof_indices, locala, INSERT_VALUES);CHKERRQ(ierr); 246 } 247 else { 248 /* get the global DOF number to appropriately set the element contribution in the RHS vector */ 249 ierr = DMMoabGetDofsLocal(dm, nconn, connect, dof_indices);CHKERRQ(ierr); 250 251 /* set the values directly into appropriate locations. Can alternately use VecSetValues */ 252 ierr = MatSetValuesLocal(A, nconn, dof_indices, nconn, dof_indices, locala, INSERT_VALUES);CHKERRQ(ierr); 253 } 254 } 255 256 /* clean up memory */ 257 ierr = PetscFree(locala);CHKERRQ(ierr); 258 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 259 260 /* finish assembly */ 261 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 262 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "DMMoabSetBlockFills_Private" 269 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w,const PetscInt *fill,PetscInt **rfill) 270 { 271 PetscErrorCode ierr; 272 PetscInt i,j,*ifill; 273 274 PetscFunctionBegin; 275 if (!fill) PetscFunctionReturn(0); 276 ierr = PetscMalloc1(w*w,&ifill);CHKERRQ(ierr); 277 for (i=0;i<w;i++) { 278 for (j=0; j<w; j++) 279 ifill[i*w+j]=fill[i*w+j]; 280 } 281 282 *rfill = ifill; 283 PetscFunctionReturn(0); 284 } 285 286 287 #undef __FUNCT__ 288 #define __FUNCT__ "DMMoabSetBlockFills" 289 /*@ 290 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 291 of the matrix returned by DMCreateMatrix(). 292 293 Logically Collective on DMDA 294 295 Input Parameter: 296 + dm - the DMMoab object 297 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 298 - ofill - the fill pattern in the off-diagonal blocks 299 300 Level: developer 301 302 Notes: This only makes sense when you are doing multicomponent problems but using the 303 MPIAIJ matrix format 304 305 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 306 representing coupling and 0 entries for missing coupling. For example 307 $ dfill[9] = {1, 0, 0, 308 $ 1, 1, 0, 309 $ 0, 1, 1} 310 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 311 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 312 diagonal block). 313 314 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 315 can be represented in the dfill, ofill format 316 317 Contributed by Glenn Hammond 318 319 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 320 321 @*/ 322 PetscErrorCode DMMoabSetBlockFills(DM dm,const PetscInt *dfill,const PetscInt *ofill) 323 { 324 DM_Moab *dmmoab=(DM_Moab*)dm->data; 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 329 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields,dfill,&dmmoab->dfill);CHKERRQ(ierr); 330 ierr = DMMoabSetBlockFills_Private(dmmoab->numFields,ofill,&dmmoab->ofill);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333