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