1 #include <petscdmsliced.h> /*I "petscdmsliced.h" I*/ 2 #include <petscmat.h> /*I "petscmat.h" I*/ 3 #include <private/dmimpl.h> /*I "petscdm.h" I*/ 4 5 /* CSR storage of the nonzero structure of a bs*bs matrix */ 6 typedef struct { 7 PetscInt bs,nz,*i,*j; 8 } DMSlicedBlockFills; 9 10 typedef struct { 11 PetscInt bs,n,N,Nghosts,*ghosts; 12 PetscInt d_nz,o_nz,*d_nnz,*o_nnz; 13 DMSlicedBlockFills *dfill,*ofill; 14 } DM_Sliced; 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "DMCreateMatrix_Sliced" 18 PetscErrorCode DMCreateMatrix_Sliced(DM dm, const MatType mtype,Mat *J) 19 { 20 PetscErrorCode ierr; 21 PetscInt *globals,*sd_nnz,*so_nnz,rstart,bs,i; 22 ISLocalToGlobalMapping lmap,blmap; 23 void (*aij)(void) = PETSC_NULL; 24 DM_Sliced *slice = (DM_Sliced*)dm->data; 25 26 PetscFunctionBegin; 27 bs = slice->bs; 28 ierr = MatCreate(((PetscObject)dm)->comm,J);CHKERRQ(ierr); 29 ierr = MatSetSizes(*J,slice->n*bs,slice->n*bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 30 ierr = MatSetType(*J,mtype);CHKERRQ(ierr); 31 ierr = MatSeqBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz);CHKERRQ(ierr); 32 ierr = MatMPIBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);CHKERRQ(ierr); 33 /* In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any 34 * good before going on with it. */ 35 ierr = PetscObjectQueryFunction((PetscObject)*J,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 36 if (!aij) { 37 ierr = PetscObjectQueryFunction((PetscObject)*J,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 38 } 39 if (aij) { 40 if (bs == 1) { 41 ierr = MatSeqAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz);CHKERRQ(ierr); 42 ierr = MatMPIAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);CHKERRQ(ierr); 43 } else if (!slice->d_nnz) { 44 ierr = MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,PETSC_NULL);CHKERRQ(ierr); 45 ierr = MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,PETSC_NULL,slice->o_nz*bs,PETSC_NULL);CHKERRQ(ierr); 46 } else { 47 /* The user has provided preallocation per block-row, convert it to per scalar-row respecting DMSlicedSetBlockFills() if applicable */ 48 ierr = PetscMalloc2(slice->n*bs,PetscInt,&sd_nnz,(!!slice->o_nnz)*slice->n*bs,PetscInt,&so_nnz);CHKERRQ(ierr); 49 for (i=0; i<slice->n*bs; i++) { 50 sd_nnz[i] = (slice->d_nnz[i/bs]-1) * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs) 51 + (slice->dfill ? slice->dfill->i[i%bs+1]-slice->dfill->i[i%bs] : bs); 52 if (so_nnz) { 53 so_nnz[i] = slice->o_nnz[i/bs] * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs); 54 } 55 } 56 ierr = MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz);CHKERRQ(ierr); 57 ierr = MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz,slice->o_nz*bs,so_nnz);CHKERRQ(ierr); 58 ierr = PetscFree2(sd_nnz,so_nnz);CHKERRQ(ierr); 59 } 60 } 61 62 ierr = MatSetBlockSize(*J,bs);CHKERRQ(ierr); 63 64 /* Set up the local to global map. For the scalar map, we have to translate to entry-wise indexing instead of block-wise. */ 65 ierr = PetscMalloc((slice->n+slice->Nghosts)*bs*sizeof(PetscInt),&globals);CHKERRQ(ierr); 66 ierr = MatGetOwnershipRange(*J,&rstart,PETSC_NULL);CHKERRQ(ierr); 67 for (i=0; i<slice->n*bs; i++) { 68 globals[i] = rstart + i; 69 } 70 for (i=0; i<slice->Nghosts*bs; i++) { 71 globals[slice->n*bs+i] = slice->ghosts[i/bs]*bs + i%bs; 72 } 73 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,(slice->n+slice->Nghosts)*bs,globals,PETSC_OWN_POINTER,&lmap);CHKERRQ(ierr); 74 ierr = ISLocalToGlobalMappingBlock(lmap,bs,&blmap);CHKERRQ(ierr); 75 ierr = MatSetLocalToGlobalMapping(*J,lmap,lmap);CHKERRQ(ierr); 76 ierr = MatSetLocalToGlobalMappingBlock(*J,blmap,blmap);CHKERRQ(ierr); 77 ierr = ISLocalToGlobalMappingDestroy(&lmap);CHKERRQ(ierr); 78 ierr = ISLocalToGlobalMappingDestroy(&blmap);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 #undef __FUNCT__ 83 #define __FUNCT__ "DMSlicedSetGhosts" 84 /*@C 85 DMSlicedSetGhosts - Sets the global indices of other processes elements that will 86 be ghosts on this process 87 88 Not Collective 89 90 Input Parameters: 91 + slice - the DM object 92 . bs - block size 93 . nlocal - number of local (owned, non-ghost) blocks 94 . Nghosts - number of ghost blocks on this process 95 - ghosts - global indices of each ghost block 96 97 Level: advanced 98 99 .seealso DMDestroy(), DMCreateGlobalVector() 100 101 @*/ 102 PetscErrorCode DMSlicedSetGhosts(DM dm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[]) 103 { 104 PetscErrorCode ierr; 105 DM_Sliced *slice = (DM_Sliced*)dm->data; 106 107 PetscFunctionBegin; 108 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 109 ierr = PetscFree(slice->ghosts);CHKERRQ(ierr); 110 ierr = PetscMalloc(Nghosts*sizeof(PetscInt),&slice->ghosts);CHKERRQ(ierr); 111 ierr = PetscMemcpy(slice->ghosts,ghosts,Nghosts*sizeof(PetscInt));CHKERRQ(ierr); 112 slice->bs = bs; 113 slice->n = nlocal; 114 slice->Nghosts = Nghosts; 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "DMSlicedSetPreallocation" 120 /*@C 121 DMSlicedSetPreallocation - sets the matrix memory preallocation for matrices computed by DMSliced 122 123 Not Collective 124 125 Input Parameters: 126 + slice - the DM object 127 . d_nz - number of block nonzeros per block row in diagonal portion of local 128 submatrix (same for all local rows) 129 . d_nnz - array containing the number of block nonzeros in the various block rows 130 of the in diagonal portion of the local (possibly different for each block 131 row) or PETSC_NULL. 132 . o_nz - number of block nonzeros per block row in the off-diagonal portion of local 133 submatrix (same for all local rows). 134 - o_nnz - array containing the number of nonzeros in the various block rows of the 135 off-diagonal portion of the local submatrix (possibly different for 136 each block row) or PETSC_NULL. 137 138 Notes: 139 See MatMPIBAIJSetPreallocation() for more details on preallocation. If a scalar matrix (AIJ) is 140 obtained with DMSlicedGetMatrix(), the correct preallocation will be set, respecting DMSlicedSetBlockFills(). 141 142 Level: advanced 143 144 .seealso DMDestroy(), DMCreateGlobalVector(), MatMPIAIJSetPreallocation(), 145 MatMPIBAIJSetPreallocation(), DMSlicedGetMatrix(), DMSlicedSetBlockFills() 146 147 @*/ 148 PetscErrorCode DMSlicedSetPreallocation(DM dm,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 149 { 150 DM_Sliced *slice = (DM_Sliced*)dm->data; 151 152 PetscFunctionBegin; 153 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 154 slice->d_nz = d_nz; 155 slice->d_nnz = (PetscInt*)d_nnz; 156 slice->o_nz = o_nz; 157 slice->o_nnz = (PetscInt*)o_nnz; 158 PetscFunctionReturn(0); 159 } 160 161 #undef __FUNCT__ 162 #define __FUNCT__ "DMSlicedSetBlockFills_Private" 163 static PetscErrorCode DMSlicedSetBlockFills_Private(PetscInt bs,const PetscInt *fill,DMSlicedBlockFills **inf) 164 { 165 PetscErrorCode ierr; 166 PetscInt i,j,nz,*fi,*fj; 167 DMSlicedBlockFills *f; 168 169 PetscFunctionBegin; 170 PetscValidPointer(inf,3); 171 if (*inf) {ierr = PetscFree3((*inf)->i,(*inf)->j,*inf);CHKERRQ(ierr);} 172 if (!fill) PetscFunctionReturn(0); 173 for (i=0,nz=0; i<bs*bs; i++) if (fill[i]) nz++; 174 ierr = PetscMalloc3(1,DMSlicedBlockFills,&f,bs+1,PetscInt,&fi,nz,PetscInt,&fj);CHKERRQ(ierr); 175 f->bs = bs; 176 f->nz = nz; 177 f->i = fi; 178 f->j = fj; 179 for (i=0,nz=0; i<bs; i++) { 180 fi[i] = nz; 181 for (j=0; j<bs; j++) if (fill[i*bs+j]) fj[nz++] = j; 182 } 183 fi[i] = nz; 184 *inf = f; 185 PetscFunctionReturn(0); 186 } 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "DMSlicedSetBlockFills" 190 /*@C 191 DMSlicedSetBlockFills - Sets the fill pattern in each block for a multi-component problem 192 of the matrix returned by DMSlicedGetMatrix(). 193 194 Logically Collective on DM 195 196 Input Parameter: 197 + sliced - the DM object 198 . dfill - the fill pattern in the diagonal block (may be PETSC_NULL, means use dense block) 199 - ofill - the fill pattern in the off-diagonal blocks 200 201 Notes: 202 This only makes sense for multicomponent problems using scalar matrix formats (AIJ). 203 See DMDASetBlockFills() for example usage. 204 205 Level: advanced 206 207 .seealso DMSlicedGetMatrix(), DMDASetBlockFills() 208 @*/ 209 PetscErrorCode DMSlicedSetBlockFills(DM dm,const PetscInt *dfill,const PetscInt *ofill) 210 { 211 DM_Sliced *slice = (DM_Sliced*)dm->data; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 216 ierr = DMSlicedSetBlockFills_Private(slice->bs,dfill,&slice->dfill);CHKERRQ(ierr); 217 ierr = DMSlicedSetBlockFills_Private(slice->bs,ofill,&slice->ofill);CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "DMDestroy_Sliced" 223 PetscErrorCode DMDestroy_Sliced(DM dm) 224 { 225 PetscErrorCode ierr; 226 DM_Sliced *slice = (DM_Sliced*)dm->data; 227 228 PetscFunctionBegin; 229 ierr = PetscFree(slice->ghosts);CHKERRQ(ierr); 230 if (slice->dfill) {ierr = PetscFree3(slice->dfill,slice->dfill->i,slice->dfill->j);CHKERRQ(ierr);} 231 if (slice->ofill) {ierr = PetscFree3(slice->ofill,slice->ofill->i,slice->ofill->j);CHKERRQ(ierr);} 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "DMCreateGlobalVector_Sliced" 237 PetscErrorCode DMCreateGlobalVector_Sliced(DM dm,Vec *gvec) 238 { 239 PetscErrorCode ierr; 240 DM_Sliced *slice = (DM_Sliced*)dm->data; 241 242 PetscFunctionBegin; 243 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 244 PetscValidPointer(gvec,2); 245 *gvec = 0; 246 ierr = VecCreateGhostBlock(((PetscObject)dm)->comm,slice->bs,slice->n*slice->bs,PETSC_DETERMINE,slice->Nghosts,slice->ghosts,gvec);CHKERRQ(ierr); 247 ierr = PetscObjectCompose((PetscObject)*gvec,"DM",(PetscObject)dm);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 EXTERN_C_BEGIN 252 #undef __FUNCT__ 253 #define __FUNCT__ "DMCreate_Sliced" 254 PetscErrorCode DMCreate_Sliced(DM p) 255 { 256 PetscErrorCode ierr; 257 DM_Sliced *slice; 258 259 PetscFunctionBegin; 260 ierr = PetscNewLog(p,DM_Sliced,&slice);CHKERRQ(ierr); 261 p->data = slice; 262 263 ierr = PetscObjectChangeTypeName((PetscObject)p,DMSLICED);CHKERRQ(ierr); 264 p->ops->createglobalvector = DMCreateGlobalVector_Sliced; 265 p->ops->creatematrix = DMCreateMatrix_Sliced; 266 p->ops->destroy = DMDestroy_Sliced; 267 PetscFunctionReturn(0); 268 } 269 EXTERN_C_END 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "DMSlicedCreate" 273 /*@C 274 DMSlicedCreate - Creates a DM object, used to manage data for a unstructured problem 275 276 Collective on MPI_Comm 277 278 Input Parameter: 279 . comm - the processors that will share the global vector 280 281 Output Parameters: 282 . slice - the slice object 283 284 Level: advanced 285 286 .seealso DMDestroy(), DMCreateGlobalVector() 287 288 @*/ 289 PetscErrorCode DMSlicedCreate(MPI_Comm comm,DM *dm) 290 { 291 PetscErrorCode ierr; 292 293 PetscFunctionBegin; 294 PetscValidPointer(dm,2); 295 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 296 ierr = DMSetType(*dm,DMSLICED);CHKERRQ(ierr); 297 PetscFunctionReturn(0); 298 } 299 300 /* Explanation of the missing functions for DMDA-style handling of the local vector: 301 302 DMSlicedCreateLocalVector() 303 DMSlicedGlobalToLocalBegin() 304 DMSlicedGlobalToLocalEnd() 305 306 There is no way to get the global form from a local form, so DMSlicedCreateLocalVector() is a memory leak without 307 external accounting for the global vector. Also, DMSliced intends the user to work with the VecGhost interface since the 308 ghosts are already ordered after the owned entries. Contrast this to a DMDA where the local vector has a special 309 ordering described by the structured grid, hence it cannot share memory with the global form. For this reason, users 310 of DMSliced should work with the global vector and use 311 312 VecGhostGetLocalForm(), VecGhostRestoreLocalForm() 313 VecGhostUpdateBegin(), VecGhostUpdateEnd() 314 315 rather than the missing DMDA-style functions. This is conceptually simpler and offers better performance than is 316 possible with the DMDA-style interface. 317 */ 318