1 2 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 3 #include <petscmat.h> /*I "petscmat.h" I*/ 4 #include <petsc-private/matimpl.h> 5 6 extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring *); 7 extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring *); 8 extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring *); 9 extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring *); 10 11 /* 12 For ghost i that may be negative or greater than the upper bound this 13 maps it into the 0:m-1 range using periodicity 14 */ 15 #define SetInRange(i,m) ((i < 0) ? m+i:((i >= m) ? i-m:i)) 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "DMDASetBlockFills_Private" 19 static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill) 20 { 21 PetscErrorCode ierr; 22 PetscInt i,j,nz,*fill; 23 24 PetscFunctionBegin; 25 if (!dfill) PetscFunctionReturn(0); 26 27 /* count number nonzeros */ 28 nz = 0; 29 for (i=0; i<w; i++) { 30 for (j=0; j<w; j++) { 31 if (dfill[w*i+j]) nz++; 32 } 33 } 34 ierr = PetscMalloc((nz + w + 1)*sizeof(PetscInt),&fill);CHKERRQ(ierr); 35 /* construct modified CSR storage of nonzero structure */ 36 /* fill[0 -- w] marks starts of each row of column indices (and end of last row) 37 so fill[1] - fill[0] gives number of nonzeros in first row etc */ 38 nz = w + 1; 39 for (i=0; i<w; i++) { 40 fill[i] = nz; 41 for (j=0; j<w; j++) { 42 if (dfill[w*i+j]) { 43 fill[nz] = j; 44 nz++; 45 } 46 } 47 } 48 fill[w] = nz; 49 50 *rfill = fill; 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "DMDASetBlockFills" 56 /*@ 57 DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem 58 of the matrix returned by DMCreateMatrix(). 59 60 Logically Collective on DMDA 61 62 Input Parameter: 63 + da - the distributed array 64 . dfill - the fill pattern in the diagonal block (may be PETSC_NULL, means use dense block) 65 - ofill - the fill pattern in the off-diagonal blocks 66 67 68 Level: developer 69 70 Notes: This only makes sense when you are doing multicomponent problems but using the 71 MPIAIJ matrix format 72 73 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 74 representing coupling and 0 entries for missing coupling. For example 75 $ dfill[9] = {1, 0, 0, 76 $ 1, 1, 0, 77 $ 0, 1, 1} 78 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 79 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 80 diagonal block). 81 82 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 83 can be represented in the dfill, ofill format 84 85 Contributed by Glenn Hammond 86 87 .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 88 89 @*/ 90 PetscErrorCode DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill) 91 { 92 DM_DA *dd = (DM_DA*)da->data; 93 PetscErrorCode ierr; 94 95 PetscFunctionBegin; 96 ierr = DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);CHKERRQ(ierr); 97 ierr = DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "DMCreateColoring_DA" 104 PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,MatType mtype,ISColoring *coloring) 105 { 106 PetscErrorCode ierr; 107 PetscInt dim,m,n,p,nc; 108 DMDABoundaryType bx,by,bz; 109 MPI_Comm comm; 110 PetscMPIInt size; 111 PetscBool isBAIJ; 112 DM_DA *dd = (DM_DA*)da->data; 113 114 PetscFunctionBegin; 115 /* 116 m 117 ------------------------------------------------------ 118 | | 119 | | 120 | ---------------------- | 121 | | | | 122 n | yn | | | 123 | | | | 124 | .--------------------- | 125 | (xs,ys) xn | 126 | . | 127 | (gxs,gys) | 128 | | 129 ----------------------------------------------------- 130 */ 131 132 /* 133 nc - number of components per grid point 134 col - number of colors needed in one direction for single component problem 135 136 */ 137 ierr = DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);CHKERRQ(ierr); 138 139 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 140 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 141 if (ctype == IS_COLORING_GHOSTED){ 142 if (size == 1) { 143 ctype = IS_COLORING_GLOBAL; 144 } else if (dim > 1){ 145 if ((m==1 && bx == DMDA_BOUNDARY_PERIODIC) || (n==1 && by == DMDA_BOUNDARY_PERIODIC) || (p==1 && bz == DMDA_BOUNDARY_PERIODIC)){ 146 SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"IS_COLORING_GHOSTED cannot be used for periodic boundary condition having both ends of the domain on the same process"); 147 } 148 } 149 } 150 151 /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ 152 matrices is for the blocks, not the individual matrix elements */ 153 ierr = PetscStrcmp(mtype,MATBAIJ,&isBAIJ);CHKERRQ(ierr); 154 if (!isBAIJ) {ierr = PetscStrcmp(mtype,MATMPIBAIJ,&isBAIJ);CHKERRQ(ierr);} 155 if (!isBAIJ) {ierr = PetscStrcmp(mtype,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr);} 156 if (isBAIJ) { 157 dd->w = 1; 158 dd->xs = dd->xs/nc; 159 dd->xe = dd->xe/nc; 160 dd->Xs = dd->Xs/nc; 161 dd->Xe = dd->Xe/nc; 162 } 163 164 /* 165 We do not provide a getcoloring function in the DMDA operations because 166 the basic DMDA does not know about matrices. We think of DMDA as being more 167 more low-level then matrices. 168 */ 169 if (dim == 1) { 170 ierr = DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 171 } else if (dim == 2) { 172 ierr = DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 173 } else if (dim == 3) { 174 ierr = DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 175 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim); 176 if (isBAIJ) { 177 dd->w = nc; 178 dd->xs = dd->xs*nc; 179 dd->xe = dd->xe*nc; 180 dd->Xs = dd->Xs*nc; 181 dd->Xe = dd->Xe*nc; 182 } 183 PetscFunctionReturn(0); 184 } 185 186 /* ---------------------------------------------------------------------------------*/ 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "DMCreateColoring_DA_2d_MPIAIJ" 190 PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 191 { 192 PetscErrorCode ierr; 193 PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col; 194 PetscInt ncolors; 195 MPI_Comm comm; 196 DMDABoundaryType bx,by; 197 DMDAStencilType st; 198 ISColoringValue *colors; 199 DM_DA *dd = (DM_DA*)da->data; 200 201 PetscFunctionBegin; 202 /* 203 nc - number of components per grid point 204 col - number of colors needed in one direction for single component problem 205 206 */ 207 ierr = DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 208 col = 2*s + 1; 209 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 210 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 211 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 212 213 /* special case as taught to us by Paul Hovland */ 214 if (st == DMDA_STENCIL_STAR && s == 1) { 215 ierr = DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 216 } else { 217 218 if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\ 219 by 2*stencil_width + 1 (%d)\n", m, col); 220 if (by == DMDA_BOUNDARY_PERIODIC && (n % col)) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\ 221 by 2*stencil_width + 1 (%d)\n", n, col); 222 if (ctype == IS_COLORING_GLOBAL) { 223 if (!dd->localcoloring) { 224 ierr = PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 225 ii = 0; 226 for (j=ys; j<ys+ny; j++) { 227 for (i=xs; i<xs+nx; i++) { 228 for (k=0; k<nc; k++) { 229 colors[ii++] = k + nc*((i % col) + col*(j % col)); 230 } 231 } 232 } 233 ncolors = nc + nc*(col-1 + col*(col-1)); 234 ierr = ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);CHKERRQ(ierr); 235 } 236 *coloring = dd->localcoloring; 237 } else if (ctype == IS_COLORING_GHOSTED) { 238 if (!dd->ghostedcoloring) { 239 ierr = PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 240 ii = 0; 241 for (j=gys; j<gys+gny; j++) { 242 for (i=gxs; i<gxs+gnx; i++) { 243 for (k=0; k<nc; k++) { 244 /* the complicated stuff is to handle periodic boundaries */ 245 colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col)); 246 } 247 } 248 } 249 ncolors = nc + nc*(col - 1 + col*(col-1)); 250 ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);CHKERRQ(ierr); 251 /* PetscIntView(ncolors,(PetscInt *)colors,0); */ 252 253 ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);CHKERRQ(ierr); 254 } 255 *coloring = dd->ghostedcoloring; 256 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 257 } 258 ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 /* ---------------------------------------------------------------------------------*/ 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "DMCreateColoring_DA_3d_MPIAIJ" 266 PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 267 { 268 PetscErrorCode ierr; 269 PetscInt xs,ys,nx,ny,i,j,gxs,gys,gnx,gny,m,n,p,dim,s,k,nc,col,zs,gzs,ii,l,nz,gnz,M,N,P; 270 PetscInt ncolors; 271 MPI_Comm comm; 272 DMDABoundaryType bx,by,bz; 273 DMDAStencilType st; 274 ISColoringValue *colors; 275 DM_DA *dd = (DM_DA*)da->data; 276 277 PetscFunctionBegin; 278 /* 279 nc - number of components per grid point 280 col - number of colors needed in one direction for single component problem 281 282 */ 283 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 284 col = 2*s + 1; 285 if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\ 286 by 2*stencil_width + 1\n"); 287 if (by == DMDA_BOUNDARY_PERIODIC && (n % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\ 288 by 2*stencil_width + 1\n"); 289 if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\ 290 by 2*stencil_width + 1\n"); 291 292 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 293 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 294 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 295 296 /* create the coloring */ 297 if (ctype == IS_COLORING_GLOBAL) { 298 if (!dd->localcoloring) { 299 ierr = PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 300 ii = 0; 301 for (k=zs; k<zs+nz; k++) { 302 for (j=ys; j<ys+ny; j++) { 303 for (i=xs; i<xs+nx; i++) { 304 for (l=0; l<nc; l++) { 305 colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col)); 306 } 307 } 308 } 309 } 310 ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1)); 311 ierr = ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,&dd->localcoloring);CHKERRQ(ierr); 312 } 313 *coloring = dd->localcoloring; 314 } else if (ctype == IS_COLORING_GHOSTED) { 315 if (!dd->ghostedcoloring) { 316 ierr = PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 317 ii = 0; 318 for (k=gzs; k<gzs+gnz; k++) { 319 for (j=gys; j<gys+gny; j++) { 320 for (i=gxs; i<gxs+gnx; i++) { 321 for (l=0; l<nc; l++) { 322 /* the complicated stuff is to handle periodic boundaries */ 323 colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col)); 324 } 325 } 326 } 327 } 328 ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1)); 329 ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,&dd->ghostedcoloring);CHKERRQ(ierr); 330 ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);CHKERRQ(ierr); 331 } 332 *coloring = dd->ghostedcoloring; 333 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 334 ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 /* ---------------------------------------------------------------------------------*/ 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "DMCreateColoring_DA_1d_MPIAIJ" 342 PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 343 { 344 PetscErrorCode ierr; 345 PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col; 346 PetscInt ncolors; 347 MPI_Comm comm; 348 DMDABoundaryType bx; 349 ISColoringValue *colors; 350 DM_DA *dd = (DM_DA*)da->data; 351 352 PetscFunctionBegin; 353 /* 354 nc - number of components per grid point 355 col - number of colors needed in one direction for single component problem 356 357 */ 358 ierr = DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 359 col = 2*s + 1; 360 361 if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\ 362 by 2*stencil_width + 1 %d\n",(int)m,(int)col); 363 364 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 365 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 366 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 367 368 /* create the coloring */ 369 if (ctype == IS_COLORING_GLOBAL) { 370 if (!dd->localcoloring) { 371 ierr = PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 372 i1 = 0; 373 for (i=xs; i<xs+nx; i++) { 374 for (l=0; l<nc; l++) { 375 colors[i1++] = l + nc*(i % col); 376 } 377 } 378 ncolors = nc + nc*(col-1); 379 ierr = ISColoringCreate(comm,ncolors,nc*nx,colors,&dd->localcoloring);CHKERRQ(ierr); 380 } 381 *coloring = dd->localcoloring; 382 } else if (ctype == IS_COLORING_GHOSTED) { 383 if (!dd->ghostedcoloring) { 384 ierr = PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 385 i1 = 0; 386 for (i=gxs; i<gxs+gnx; i++) { 387 for (l=0; l<nc; l++) { 388 /* the complicated stuff is to handle periodic boundaries */ 389 colors[i1++] = l + nc*(SetInRange(i,m) % col); 390 } 391 } 392 ncolors = nc + nc*(col-1); 393 ierr = ISColoringCreate(comm,ncolors,nc*gnx,colors,&dd->ghostedcoloring);CHKERRQ(ierr); 394 ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);CHKERRQ(ierr); 395 } 396 *coloring = dd->ghostedcoloring; 397 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 398 ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "DMCreateColoring_DA_2d_5pt_MPIAIJ" 404 PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 405 { 406 PetscErrorCode ierr; 407 PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc; 408 PetscInt ncolors; 409 MPI_Comm comm; 410 DMDABoundaryType bx,by; 411 ISColoringValue *colors; 412 DM_DA *dd = (DM_DA*)da->data; 413 414 PetscFunctionBegin; 415 /* 416 nc - number of components per grid point 417 col - number of colors needed in one direction for single component problem 418 419 */ 420 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);CHKERRQ(ierr); 421 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 422 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 423 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 424 425 if (bx == DMDA_BOUNDARY_PERIODIC && (m % 5)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible by 5\n"); 426 if (by == DMDA_BOUNDARY_PERIODIC && (n % 5)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible by 5\n"); 427 428 /* create the coloring */ 429 if (ctype == IS_COLORING_GLOBAL) { 430 if (!dd->localcoloring) { 431 ierr = PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 432 ii = 0; 433 for (j=ys; j<ys+ny; j++) { 434 for (i=xs; i<xs+nx; i++) { 435 for (k=0; k<nc; k++) { 436 colors[ii++] = k + nc*((3*j+i) % 5); 437 } 438 } 439 } 440 ncolors = 5*nc; 441 ierr = ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);CHKERRQ(ierr); 442 } 443 *coloring = dd->localcoloring; 444 } else if (ctype == IS_COLORING_GHOSTED) { 445 if (!dd->ghostedcoloring) { 446 ierr = PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 447 ii = 0; 448 for (j=gys; j<gys+gny; j++) { 449 for (i=gxs; i<gxs+gnx; i++) { 450 for (k=0; k<nc; k++) { 451 colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5); 452 } 453 } 454 } 455 ncolors = 5*nc; 456 ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);CHKERRQ(ierr); 457 ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);CHKERRQ(ierr); 458 } 459 *coloring = dd->ghostedcoloring; 460 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 461 PetscFunctionReturn(0); 462 } 463 464 /* =========================================================================== */ 465 extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat); 466 extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat); 467 extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat); 468 extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat); 469 extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat); 470 extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat); 471 extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat); 472 extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat); 473 extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat); 474 extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat); 475 476 #undef __FUNCT__ 477 #define __FUNCT__ "MatSetupDM" 478 /*@C 479 MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix 480 481 Logically Collective on Mat 482 483 Input Parameters: 484 + mat - the matrix 485 - da - the da 486 487 Level: intermediate 488 489 @*/ 490 PetscErrorCode MatSetupDM(Mat mat,DM da) 491 { 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 496 PetscValidHeaderSpecific(da,DM_CLASSID,1); 497 ierr = PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 EXTERN_C_BEGIN 502 #undef __FUNCT__ 503 #define __FUNCT__ "MatView_MPI_DA" 504 PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer) 505 { 506 DM da; 507 PetscErrorCode ierr; 508 const char *prefix; 509 Mat Anatural; 510 AO ao; 511 PetscInt rstart,rend,*petsc,i; 512 IS is; 513 MPI_Comm comm; 514 PetscViewerFormat format; 515 516 PetscFunctionBegin; 517 /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */ 518 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 519 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 520 521 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 522 ierr = MatGetDM(A, &da);CHKERRQ(ierr); 523 if (!da) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 524 525 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 526 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 527 ierr = PetscMalloc((rend-rstart)*sizeof(PetscInt),&petsc);CHKERRQ(ierr); 528 for (i=rstart; i<rend; i++) petsc[i-rstart] = i; 529 ierr = AOApplicationToPetsc(ao,rend-rstart,petsc);CHKERRQ(ierr); 530 ierr = ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 531 532 /* call viewer on natural ordering */ 533 ierr = MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);CHKERRQ(ierr); 534 ierr = ISDestroy(&is);CHKERRQ(ierr); 535 ierr = PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);CHKERRQ(ierr); 536 ierr = PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);CHKERRQ(ierr); 537 ierr = PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);CHKERRQ(ierr); 538 ierr = MatView(Anatural,viewer);CHKERRQ(ierr); 539 ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 540 PetscFunctionReturn(0); 541 } 542 EXTERN_C_END 543 544 EXTERN_C_BEGIN 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatLoad_MPI_DA" 547 PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer) 548 { 549 DM da; 550 PetscErrorCode ierr; 551 Mat Anatural,Aapp; 552 AO ao; 553 PetscInt rstart,rend,*app,i; 554 IS is; 555 MPI_Comm comm; 556 557 PetscFunctionBegin; 558 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 559 ierr = MatGetDM(A, &da);CHKERRQ(ierr); 560 if (!da) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 561 562 /* Load the matrix in natural ordering */ 563 ierr = MatCreate(((PetscObject)A)->comm,&Anatural);CHKERRQ(ierr); 564 ierr = MatSetType(Anatural,((PetscObject)A)->type_name);CHKERRQ(ierr); 565 ierr = MatSetSizes(Anatural,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 566 ierr = MatLoad(Anatural,viewer);CHKERRQ(ierr); 567 568 /* Map natural ordering to application ordering and create IS */ 569 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 570 ierr = MatGetOwnershipRange(Anatural,&rstart,&rend);CHKERRQ(ierr); 571 ierr = PetscMalloc((rend-rstart)*sizeof(PetscInt),&app);CHKERRQ(ierr); 572 for (i=rstart; i<rend; i++) app[i-rstart] = i; 573 ierr = AOPetscToApplication(ao,rend-rstart,app);CHKERRQ(ierr); 574 ierr = ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 575 576 /* Do permutation and replace header */ 577 ierr = MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);CHKERRQ(ierr); 578 ierr = MatHeaderReplace(A,Aapp);CHKERRQ(ierr); 579 ierr = ISDestroy(&is);CHKERRQ(ierr); 580 ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 EXTERN_C_END 584 585 #undef __FUNCT__ 586 #define __FUNCT__ "DMCreateMatrix_DA" 587 PetscErrorCode DMCreateMatrix_DA(DM da, MatType mtype,Mat *J) 588 { 589 PetscErrorCode ierr; 590 PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P; 591 Mat A; 592 MPI_Comm comm; 593 MatType Atype; 594 PetscSection section, sectionGlobal; 595 void (*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL; 596 MatType ttype[256]; 597 PetscBool flg; 598 PetscMPIInt size; 599 DM_DA *dd = (DM_DA*)da->data; 600 601 PetscFunctionBegin; 602 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 603 ierr = MatInitializePackage(PETSC_NULL);CHKERRQ(ierr); 604 #endif 605 if (!mtype) mtype = MATAIJ; 606 ierr = PetscStrcpy((char*)ttype,mtype);CHKERRQ(ierr); 607 ierr = PetscOptionsBegin(((PetscObject)da)->comm,((PetscObject)da)->prefix,"DMDA options","Mat");CHKERRQ(ierr); 608 ierr = PetscOptionsList("-dm_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)ttype,256,&flg);CHKERRQ(ierr); 609 ierr = PetscOptionsEnd(); 610 611 ierr = DMGetDefaultSection(da, §ion);CHKERRQ(ierr); 612 if (section) { 613 PetscInt bs = -1; 614 PetscInt localSize; 615 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric; 616 617 ierr = DMGetDefaultGlobalSection(da, §ionGlobal);CHKERRQ(ierr); 618 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);CHKERRQ(ierr); 619 ierr = MatCreate(((PetscObject) da)->comm, J);CHKERRQ(ierr); 620 ierr = MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 621 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 622 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 623 ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr); 624 ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr); 625 ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr); 626 ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr); 627 ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr); 628 ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr); 629 ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr); 630 /* Check for symmetric storage */ 631 isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock); 632 if (isSymmetric) { 633 ierr = MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr); 634 } 635 if (!isShell) { 636 /* PetscBool fillMatrix = (PetscBool) !da->prealloc_only; */ 637 PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal; 638 639 if (bs < 0) { 640 if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) { 641 PetscInt pStart, pEnd, p, dof; 642 643 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 644 for (p = pStart; p < pEnd; ++p) { 645 ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr); 646 if (dof) { 647 bs = dof; 648 break; 649 } 650 } 651 } else { 652 bs = 1; 653 } 654 /* Must have same blocksize on all procs (some might have no points) */ 655 bsLocal = bs; 656 ierr = MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, ((PetscObject) da)->comm);CHKERRQ(ierr); 657 } 658 ierr = PetscMalloc4(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz, localSize/bs, PetscInt, &dnzu, localSize/bs, PetscInt, &onzu);CHKERRQ(ierr); 659 ierr = PetscMemzero(dnz, localSize/bs * sizeof(PetscInt));CHKERRQ(ierr); 660 ierr = PetscMemzero(onz, localSize/bs * sizeof(PetscInt));CHKERRQ(ierr); 661 ierr = PetscMemzero(dnzu, localSize/bs * sizeof(PetscInt));CHKERRQ(ierr); 662 ierr = PetscMemzero(onzu, localSize/bs * sizeof(PetscInt));CHKERRQ(ierr); 663 /* ierr = DMComplexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);CHKERRQ(ierr); */ 664 ierr = PetscFree4(dnz, onz, dnzu, onzu);CHKERRQ(ierr); 665 } 666 } 667 /* 668 m 669 ------------------------------------------------------ 670 | | 671 | | 672 | ---------------------- | 673 | | | | 674 n | ny | | | 675 | | | | 676 | .--------------------- | 677 | (xs,ys) nx | 678 | . | 679 | (gxs,gys) | 680 | | 681 ----------------------------------------------------- 682 */ 683 684 /* 685 nc - number of components per grid point 686 col - number of colors needed in one direction for single component problem 687 688 */ 689 M = dd->M; 690 N = dd->N; 691 P = dd->P; 692 dim = dd->dim; 693 dof = dd->w; 694 /* ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); */ 695 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 696 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 697 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 698 ierr = MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);CHKERRQ(ierr); 699 ierr = MatSetType(A,(MatType)ttype);CHKERRQ(ierr); 700 ierr = MatSetDM(A,da);CHKERRQ(ierr); 701 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 702 ierr = MatGetType(A,&Atype);CHKERRQ(ierr); 703 /* 704 We do not provide a getmatrix function in the DMDA operations because 705 the basic DMDA does not know about matrices. We think of DMDA as being more 706 more low-level than matrices. This is kind of cheating but, cause sometimes 707 we think of DMDA has higher level than matrices. 708 709 We could switch based on Atype (or mtype), but we do not since the 710 specialized setting routines depend only the particular preallocation 711 details of the matrix, not the type itself. 712 */ 713 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 714 if (!aij) { 715 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 716 } 717 if (!aij) { 718 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 719 if (!baij) { 720 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 721 } 722 if (!baij){ 723 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 724 if (!sbaij) { 725 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 726 } 727 } 728 } 729 if (aij) { 730 if (dim == 1) { 731 if (dd->ofill) { 732 ierr = DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 733 } else { 734 ierr = DMCreateMatrix_DA_1d_MPIAIJ(da,A);CHKERRQ(ierr); 735 } 736 } else if (dim == 2) { 737 if (dd->ofill) { 738 ierr = DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 739 } else { 740 ierr = DMCreateMatrix_DA_2d_MPIAIJ(da,A);CHKERRQ(ierr); 741 } 742 } else if (dim == 3) { 743 if (dd->ofill) { 744 ierr = DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 745 } else { 746 ierr = DMCreateMatrix_DA_3d_MPIAIJ(da,A);CHKERRQ(ierr); 747 } 748 } 749 } else if (baij) { 750 if (dim == 2) { 751 ierr = DMCreateMatrix_DA_2d_MPIBAIJ(da,A);CHKERRQ(ierr); 752 } else if (dim == 3) { 753 ierr = DMCreateMatrix_DA_3d_MPIBAIJ(da,A);CHKERRQ(ierr); 754 } else SETERRQ3(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim); 755 } else if (sbaij) { 756 if (dim == 2) { 757 ierr = DMCreateMatrix_DA_2d_MPISBAIJ(da,A);CHKERRQ(ierr); 758 } else if (dim == 3) { 759 ierr = DMCreateMatrix_DA_3d_MPISBAIJ(da,A);CHKERRQ(ierr); 760 } else SETERRQ3(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim); 761 } else { 762 ISLocalToGlobalMapping ltog,ltogb; 763 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 764 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 765 ierr = MatSetUp(A);CHKERRQ(ierr); 766 ierr = MatSetLocalToGlobalMapping(A,ltog,ltog);CHKERRQ(ierr); 767 ierr = MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);CHKERRQ(ierr); 768 } 769 ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); 770 ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); 771 ierr = MatSetDM(A,da);CHKERRQ(ierr); 772 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 773 if (size > 1) { 774 /* change viewer to display matrix in natural ordering */ 775 ierr = MatShellSetOperation(A, MATOP_VIEW, (void (*)(void)) MatView_MPI_DA);CHKERRQ(ierr); 776 ierr = MatShellSetOperation(A, MATOP_LOAD, (void (*)(void)) MatLoad_MPI_DA);CHKERRQ(ierr); 777 } 778 *J = A; 779 PetscFunctionReturn(0); 780 } 781 782 /* ---------------------------------------------------------------------------------*/ 783 #undef __FUNCT__ 784 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ" 785 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J) 786 { 787 PetscErrorCode ierr; 788 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p; 789 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 790 MPI_Comm comm; 791 PetscScalar *values; 792 DMDABoundaryType bx,by; 793 ISLocalToGlobalMapping ltog,ltogb; 794 DMDAStencilType st; 795 796 PetscFunctionBegin; 797 /* 798 nc - number of components per grid point 799 col - number of colors needed in one direction for single component problem 800 801 */ 802 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 803 col = 2*s + 1; 804 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 805 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 806 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 807 808 ierr = PetscMalloc2(nc,PetscInt,&rows,col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 809 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 810 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 811 812 /* determine the matrix preallocation information */ 813 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 814 for (i=xs; i<xs+nx; i++) { 815 816 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 817 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 818 819 for (j=ys; j<ys+ny; j++) { 820 slot = i - gxs + gnx*(j - gys); 821 822 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 823 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 824 825 cnt = 0; 826 for (k=0; k<nc; k++) { 827 for (l=lstart; l<lend+1; l++) { 828 for (p=pstart; p<pend+1; p++) { 829 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 830 cols[cnt++] = k + nc*(slot + gnx*l + p); 831 } 832 } 833 } 834 rows[k] = k + nc*(slot); 835 } 836 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 837 } 838 } 839 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 840 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 841 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 842 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 843 844 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 845 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 846 847 /* 848 For each node in the grid: we get the neighbors in the local (on processor ordering 849 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 850 PETSc ordering. 851 */ 852 if (!da->prealloc_only) { 853 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 854 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 855 for (i=xs; i<xs+nx; i++) { 856 857 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 858 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 859 860 for (j=ys; j<ys+ny; j++) { 861 slot = i - gxs + gnx*(j - gys); 862 863 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 864 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 865 866 cnt = 0; 867 for (k=0; k<nc; k++) { 868 for (l=lstart; l<lend+1; l++) { 869 for (p=pstart; p<pend+1; p++) { 870 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 871 cols[cnt++] = k + nc*(slot + gnx*l + p); 872 } 873 } 874 } 875 rows[k] = k + nc*(slot); 876 } 877 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 878 } 879 } 880 ierr = PetscFree(values);CHKERRQ(ierr); 881 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 882 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 883 } 884 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 888 #undef __FUNCT__ 889 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ_Fill" 890 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J) 891 { 892 PetscErrorCode ierr; 893 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 894 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p; 895 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 896 DM_DA *dd = (DM_DA*)da->data; 897 PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill; 898 MPI_Comm comm; 899 PetscScalar *values; 900 DMDABoundaryType bx,by; 901 ISLocalToGlobalMapping ltog,ltogb; 902 DMDAStencilType st; 903 904 PetscFunctionBegin; 905 /* 906 nc - number of components per grid point 907 col - number of colors needed in one direction for single component problem 908 909 */ 910 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 911 col = 2*s + 1; 912 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 913 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 914 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 915 916 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 917 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 918 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 919 920 /* determine the matrix preallocation information */ 921 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 922 for (i=xs; i<xs+nx; i++) { 923 924 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 925 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 926 927 for (j=ys; j<ys+ny; j++) { 928 slot = i - gxs + gnx*(j - gys); 929 930 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 931 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 932 933 for (k=0; k<nc; k++) { 934 cnt = 0; 935 for (l=lstart; l<lend+1; l++) { 936 for (p=pstart; p<pend+1; p++) { 937 if (l || p) { 938 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 939 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) 940 cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 941 } 942 } else { 943 if (dfill) { 944 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) 945 cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 946 } else { 947 for (ifill_col=0; ifill_col<nc; ifill_col++) 948 cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 949 } 950 } 951 } 952 } 953 row = k + nc*(slot); 954 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 955 } 956 } 957 } 958 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 959 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 960 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 961 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 962 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 963 964 /* 965 For each node in the grid: we get the neighbors in the local (on processor ordering 966 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 967 PETSc ordering. 968 */ 969 if (!da->prealloc_only) { 970 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 971 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 972 for (i=xs; i<xs+nx; i++) { 973 974 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 975 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 976 977 for (j=ys; j<ys+ny; j++) { 978 slot = i - gxs + gnx*(j - gys); 979 980 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 981 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 982 983 for (k=0; k<nc; k++) { 984 cnt = 0; 985 for (l=lstart; l<lend+1; l++) { 986 for (p=pstart; p<pend+1; p++) { 987 if (l || p) { 988 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 989 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) 990 cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 991 } 992 } else { 993 if (dfill) { 994 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) 995 cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 996 } else { 997 for (ifill_col=0; ifill_col<nc; ifill_col++) 998 cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 999 } 1000 } 1001 } 1002 } 1003 row = k + nc*(slot); 1004 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1005 } 1006 } 1007 } 1008 ierr = PetscFree(values);CHKERRQ(ierr); 1009 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1010 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 } 1012 ierr = PetscFree(cols);CHKERRQ(ierr); 1013 PetscFunctionReturn(0); 1014 } 1015 1016 /* ---------------------------------------------------------------------------------*/ 1017 1018 #undef __FUNCT__ 1019 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ" 1020 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 1021 { 1022 PetscErrorCode ierr; 1023 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1024 PetscInt m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p,*dnz = PETSC_NULL,*onz = PETSC_NULL; 1025 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1026 MPI_Comm comm; 1027 PetscScalar *values; 1028 DMDABoundaryType bx,by,bz; 1029 ISLocalToGlobalMapping ltog,ltogb; 1030 DMDAStencilType st; 1031 1032 PetscFunctionBegin; 1033 /* 1034 nc - number of components per grid point 1035 col - number of colors needed in one direction for single component problem 1036 1037 */ 1038 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1039 col = 2*s + 1; 1040 1041 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1042 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1043 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1044 1045 ierr = PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 1046 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1047 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1048 1049 /* determine the matrix preallocation information */ 1050 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1051 for (i=xs; i<xs+nx; i++) { 1052 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1053 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1054 for (j=ys; j<ys+ny; j++) { 1055 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1056 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1057 for (k=zs; k<zs+nz; k++) { 1058 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1059 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1060 1061 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1062 1063 cnt = 0; 1064 for (l=0; l<nc; l++) { 1065 for (ii=istart; ii<iend+1; ii++) { 1066 for (jj=jstart; jj<jend+1; jj++) { 1067 for (kk=kstart; kk<kend+1; kk++) { 1068 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1069 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1070 } 1071 } 1072 } 1073 } 1074 rows[l] = l + nc*(slot); 1075 } 1076 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1077 } 1078 } 1079 } 1080 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1081 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1082 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1083 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1084 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1085 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1086 1087 /* 1088 For each node in the grid: we get the neighbors in the local (on processor ordering 1089 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1090 PETSc ordering. 1091 */ 1092 if (!da->prealloc_only) { 1093 ierr = PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1094 ierr = PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1095 for (i=xs; i<xs+nx; i++) { 1096 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1097 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1098 for (j=ys; j<ys+ny; j++) { 1099 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1100 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1101 for (k=zs; k<zs+nz; k++) { 1102 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1103 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1104 1105 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1106 1107 cnt = 0; 1108 for (l=0; l<nc; l++) { 1109 for (ii=istart; ii<iend+1; ii++) { 1110 for (jj=jstart; jj<jend+1; jj++) { 1111 for (kk=kstart; kk<kend+1; kk++) { 1112 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1113 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1114 } 1115 } 1116 } 1117 } 1118 rows[l] = l + nc*(slot); 1119 } 1120 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1121 } 1122 } 1123 } 1124 ierr = PetscFree(values);CHKERRQ(ierr); 1125 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1126 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1127 } 1128 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 /* ---------------------------------------------------------------------------------*/ 1133 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ_Fill" 1136 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J) 1137 { 1138 PetscErrorCode ierr; 1139 DM_DA *dd = (DM_DA*)da->data; 1140 PetscInt xs,nx,i,j,gxs,gnx,row,k,l; 1141 PetscInt m,dim,s,*cols = PETSC_NULL,nc,col,cnt,*ocols; 1142 PetscInt *ofill = dd->ofill; 1143 PetscScalar *values; 1144 DMDABoundaryType bx; 1145 ISLocalToGlobalMapping ltog,ltogb; 1146 PetscMPIInt rank,size; 1147 1148 PetscFunctionBegin; 1149 if (dd->bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions"); 1150 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 1151 ierr = MPI_Comm_size(((PetscObject)da)->comm,&size);CHKERRQ(ierr); 1152 1153 /* 1154 nc - number of components per grid point 1155 col - number of colors needed in one direction for single component problem 1156 1157 */ 1158 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1159 col = 2*s + 1; 1160 1161 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1162 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1163 1164 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1165 ierr = PetscMalloc2(nx*nc,PetscInt,&cols,nx*nc,PetscInt,&ocols);CHKERRQ(ierr); 1166 ierr = PetscMemzero(cols,nx*nc*sizeof(PetscInt));CHKERRQ(ierr); 1167 ierr = PetscMemzero(ocols,nx*nc*sizeof(PetscInt));CHKERRQ(ierr); 1168 1169 /* 1170 note should be smaller for first and last process with no periodic 1171 does not handle dfill 1172 */ 1173 cnt = 0; 1174 /* coupling with process to the left */ 1175 for (i=0; i<s; i++) { 1176 for (j=0; j<nc; j++){ 1177 ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j])); 1178 cols[cnt] = nc + (s + i)*(ofill[j+1] - ofill[j]); 1179 cnt++; 1180 } 1181 } 1182 for (i=s; i<nx-s; i++) { 1183 for (j=0; j<nc; j++){ 1184 cols[cnt] = nc + 2*s*(ofill[j+1] - ofill[j]); 1185 cnt++; 1186 } 1187 } 1188 /* coupling with process to the right */ 1189 for (i=nx-s; i<nx; i++){ 1190 for (j=0; j<nc; j++){ 1191 ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j])); 1192 cols[cnt] = nc + (s + nx - i - 1)*(ofill[j+1] - ofill[j]); 1193 cnt++; 1194 } 1195 } 1196 1197 ierr = MatSeqAIJSetPreallocation(J,0,cols);CHKERRQ(ierr); 1198 ierr = MatMPIAIJSetPreallocation(J,0,cols,0,ocols);CHKERRQ(ierr); 1199 ierr = PetscFree2(cols,ocols);CHKERRQ(ierr); 1200 1201 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1202 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1203 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1204 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1205 1206 /* 1207 For each node in the grid: we get the neighbors in the local (on processor ordering 1208 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1209 PETSc ordering. 1210 */ 1211 if (!da->prealloc_only) { 1212 ierr = PetscMalloc(col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1213 ierr = PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1214 ierr = PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1215 1216 row = xs*nc; 1217 /* coupling with process to the left */ 1218 for (i=xs; i<xs+s; i++) { 1219 for (j=0; j<nc; j++){ 1220 cnt = 0; 1221 if (rank) { 1222 for (l=0; l<s; l++) { 1223 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1224 } 1225 } 1226 for (k=0; k<nc; k++) { 1227 cols[cnt++] = i*nc + k; 1228 } 1229 for (l=0; l<s; l++) { 1230 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1231 } 1232 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1233 row++; 1234 } 1235 } 1236 for (i=xs+s; i<xs+nx-s; i++) { 1237 for (j=0; j<nc; j++){ 1238 cnt = 0; 1239 for (l=0; l<s; l++) { 1240 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1241 } 1242 for (k=0; k<nc; k++) { 1243 cols[cnt++] = i*nc + k; 1244 } 1245 for (l=0; l<s; l++) { 1246 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1247 } 1248 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1249 row++; 1250 } 1251 } 1252 /* coupling with process to the right */ 1253 for (i=xs+nx-s; i<xs+nx; i++){ 1254 for (j=0; j<nc; j++){ 1255 cnt = 0; 1256 for (l=0; l<s; l++) { 1257 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1258 } 1259 for (k=0; k<nc; k++) { 1260 cols[cnt++] = i*nc + k; 1261 } 1262 if (rank < size-1) { 1263 for (l=0; l<s; l++) { 1264 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1265 } 1266 } 1267 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1268 row++; 1269 } 1270 } 1271 ierr = PetscFree(values);CHKERRQ(ierr); 1272 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1273 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1274 ierr = PetscFree(cols);CHKERRQ(ierr); 1275 } 1276 PetscFunctionReturn(0); 1277 } 1278 1279 /* ---------------------------------------------------------------------------------*/ 1280 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ" 1283 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 1284 { 1285 PetscErrorCode ierr; 1286 PetscInt xs,nx,i,i1,slot,gxs,gnx; 1287 PetscInt m,dim,s,*cols = PETSC_NULL,nc,*rows = PETSC_NULL,col,cnt,l; 1288 PetscInt istart,iend; 1289 PetscScalar *values; 1290 DMDABoundaryType bx; 1291 ISLocalToGlobalMapping ltog,ltogb; 1292 1293 PetscFunctionBegin; 1294 /* 1295 nc - number of components per grid point 1296 col - number of colors needed in one direction for single component problem 1297 1298 */ 1299 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1300 col = 2*s + 1; 1301 1302 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1303 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1304 1305 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1306 ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 1307 ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 1308 1309 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1310 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1311 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1312 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1313 1314 /* 1315 For each node in the grid: we get the neighbors in the local (on processor ordering 1316 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1317 PETSc ordering. 1318 */ 1319 if (!da->prealloc_only) { 1320 ierr = PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 1321 ierr = PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1322 ierr = PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1323 for (i=xs; i<xs+nx; i++) { 1324 istart = PetscMax(-s,gxs - i); 1325 iend = PetscMin(s,gxs + gnx - i - 1); 1326 slot = i - gxs; 1327 1328 cnt = 0; 1329 for (l=0; l<nc; l++) { 1330 for (i1=istart; i1<iend+1; i1++) { 1331 cols[cnt++] = l + nc*(slot + i1); 1332 } 1333 rows[l] = l + nc*(slot); 1334 } 1335 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1336 } 1337 ierr = PetscFree(values);CHKERRQ(ierr); 1338 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1339 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1340 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1341 } 1342 PetscFunctionReturn(0); 1343 } 1344 1345 #undef __FUNCT__ 1346 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIBAIJ" 1347 PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 1348 { 1349 PetscErrorCode ierr; 1350 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1351 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1352 PetscInt istart,iend,jstart,jend,ii,jj; 1353 MPI_Comm comm; 1354 PetscScalar *values; 1355 DMDABoundaryType bx,by; 1356 DMDAStencilType st; 1357 ISLocalToGlobalMapping ltog,ltogb; 1358 1359 PetscFunctionBegin; 1360 /* 1361 nc - number of components per grid point 1362 col - number of colors needed in one direction for single component problem 1363 */ 1364 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1365 col = 2*s + 1; 1366 1367 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1368 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1369 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1370 1371 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1372 1373 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1374 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1375 1376 /* determine the matrix preallocation information */ 1377 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1378 for (i=xs; i<xs+nx; i++) { 1379 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1380 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1381 for (j=ys; j<ys+ny; j++) { 1382 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1383 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1384 slot = i - gxs + gnx*(j - gys); 1385 1386 /* Find block columns in block row */ 1387 cnt = 0; 1388 for (ii=istart; ii<iend+1; ii++) { 1389 for (jj=jstart; jj<jend+1; jj++) { 1390 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1391 cols[cnt++] = slot + ii + gnx*jj; 1392 } 1393 } 1394 } 1395 ierr = MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);CHKERRQ(ierr); 1396 } 1397 } 1398 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1399 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1400 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1401 1402 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1403 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1404 1405 /* 1406 For each node in the grid: we get the neighbors in the local (on processor ordering 1407 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1408 PETSc ordering. 1409 */ 1410 if (!da->prealloc_only) { 1411 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1412 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1413 for (i=xs; i<xs+nx; i++) { 1414 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1415 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1416 for (j=ys; j<ys+ny; j++) { 1417 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1418 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1419 slot = i - gxs + gnx*(j - gys); 1420 cnt = 0; 1421 for (ii=istart; ii<iend+1; ii++) { 1422 for (jj=jstart; jj<jend+1; jj++) { 1423 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1424 cols[cnt++] = slot + ii + gnx*jj; 1425 } 1426 } 1427 } 1428 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1429 } 1430 } 1431 ierr = PetscFree(values);CHKERRQ(ierr); 1432 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1433 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1434 } 1435 ierr = PetscFree(cols);CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIBAIJ" 1441 PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 1442 { 1443 PetscErrorCode ierr; 1444 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1445 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1446 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1447 MPI_Comm comm; 1448 PetscScalar *values; 1449 DMDABoundaryType bx,by,bz; 1450 DMDAStencilType st; 1451 ISLocalToGlobalMapping ltog,ltogb; 1452 1453 PetscFunctionBegin; 1454 /* 1455 nc - number of components per grid point 1456 col - number of colors needed in one direction for single component problem 1457 1458 */ 1459 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1460 col = 2*s + 1; 1461 1462 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1463 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1464 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1465 1466 ierr = PetscMalloc(col*col*col*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1467 1468 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1469 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1470 1471 /* determine the matrix preallocation information */ 1472 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1473 for (i=xs; i<xs+nx; i++) { 1474 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1475 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1476 for (j=ys; j<ys+ny; j++) { 1477 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1478 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1479 for (k=zs; k<zs+nz; k++) { 1480 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1481 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1482 1483 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1484 1485 /* Find block columns in block row */ 1486 cnt = 0; 1487 for (ii=istart; ii<iend+1; ii++) { 1488 for (jj=jstart; jj<jend+1; jj++) { 1489 for (kk=kstart; kk<kend+1; kk++) { 1490 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1491 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1492 } 1493 } 1494 } 1495 } 1496 ierr = MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);CHKERRQ(ierr); 1497 } 1498 } 1499 } 1500 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1501 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1502 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1503 1504 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1505 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1506 1507 /* 1508 For each node in the grid: we get the neighbors in the local (on processor ordering 1509 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1510 PETSc ordering. 1511 */ 1512 if (!da->prealloc_only) { 1513 ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1514 ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1515 for (i=xs; i<xs+nx; i++) { 1516 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1517 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1518 for (j=ys; j<ys+ny; j++) { 1519 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1520 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1521 for (k=zs; k<zs+nz; k++) { 1522 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1523 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1524 1525 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1526 1527 cnt = 0; 1528 for (ii=istart; ii<iend+1; ii++) { 1529 for (jj=jstart; jj<jend+1; jj++) { 1530 for (kk=kstart; kk<kend+1; kk++) { 1531 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1532 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1533 } 1534 } 1535 } 1536 } 1537 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1538 } 1539 } 1540 } 1541 ierr = PetscFree(values);CHKERRQ(ierr); 1542 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1543 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 } 1545 ierr = PetscFree(cols);CHKERRQ(ierr); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 #undef __FUNCT__ 1550 #define __FUNCT__ "L2GFilterUpperTriangular" 1551 /* 1552 This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 1553 identify in the local ordering with periodic domain. 1554 */ 1555 static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 1556 { 1557 PetscErrorCode ierr; 1558 PetscInt i,n; 1559 1560 PetscFunctionBegin; 1561 ierr = ISLocalToGlobalMappingApply(ltog,1,row,row);CHKERRQ(ierr); 1562 ierr = ISLocalToGlobalMappingApply(ltog,*cnt,col,col);CHKERRQ(ierr); 1563 for (i=0,n=0; i<*cnt; i++) { 1564 if (col[i] >= *row) col[n++] = col[i]; 1565 } 1566 *cnt = n; 1567 PetscFunctionReturn(0); 1568 } 1569 1570 #undef __FUNCT__ 1571 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPISBAIJ" 1572 PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 1573 { 1574 PetscErrorCode ierr; 1575 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1576 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1577 PetscInt istart,iend,jstart,jend,ii,jj; 1578 MPI_Comm comm; 1579 PetscScalar *values; 1580 DMDABoundaryType bx,by; 1581 DMDAStencilType st; 1582 ISLocalToGlobalMapping ltog,ltogb; 1583 1584 PetscFunctionBegin; 1585 /* 1586 nc - number of components per grid point 1587 col - number of colors needed in one direction for single component problem 1588 */ 1589 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1590 col = 2*s + 1; 1591 1592 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1593 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1594 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1595 1596 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1597 1598 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1599 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1600 1601 /* determine the matrix preallocation information */ 1602 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1603 for (i=xs; i<xs+nx; i++) { 1604 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1605 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1606 for (j=ys; j<ys+ny; j++) { 1607 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1608 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1609 slot = i - gxs + gnx*(j - gys); 1610 1611 /* Find block columns in block row */ 1612 cnt = 0; 1613 for (ii=istart; ii<iend+1; ii++) { 1614 for (jj=jstart; jj<jend+1; jj++) { 1615 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1616 cols[cnt++] = slot + ii + gnx*jj; 1617 } 1618 } 1619 } 1620 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1621 ierr = MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1622 } 1623 } 1624 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1625 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1626 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1627 1628 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1629 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1630 1631 /* 1632 For each node in the grid: we get the neighbors in the local (on processor ordering 1633 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1634 PETSc ordering. 1635 */ 1636 if (!da->prealloc_only) { 1637 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1638 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1639 for (i=xs; i<xs+nx; i++) { 1640 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1641 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1642 for (j=ys; j<ys+ny; j++) { 1643 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1644 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1645 slot = i - gxs + gnx*(j - gys); 1646 1647 /* Find block columns in block row */ 1648 cnt = 0; 1649 for (ii=istart; ii<iend+1; ii++) { 1650 for (jj=jstart; jj<jend+1; jj++) { 1651 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1652 cols[cnt++] = slot + ii + gnx*jj; 1653 } 1654 } 1655 } 1656 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1657 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1658 } 1659 } 1660 ierr = PetscFree(values);CHKERRQ(ierr); 1661 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1662 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1663 } 1664 ierr = PetscFree(cols);CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 #undef __FUNCT__ 1669 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPISBAIJ" 1670 PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 1671 { 1672 PetscErrorCode ierr; 1673 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1674 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1675 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1676 MPI_Comm comm; 1677 PetscScalar *values; 1678 DMDABoundaryType bx,by,bz; 1679 DMDAStencilType st; 1680 ISLocalToGlobalMapping ltog,ltogb; 1681 1682 PetscFunctionBegin; 1683 /* 1684 nc - number of components per grid point 1685 col - number of colors needed in one direction for single component problem 1686 */ 1687 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1688 col = 2*s + 1; 1689 1690 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1691 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1692 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1693 1694 /* create the matrix */ 1695 ierr = PetscMalloc(col*col*col*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1696 1697 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1698 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1699 1700 /* determine the matrix preallocation information */ 1701 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1702 for (i=xs; i<xs+nx; i++) { 1703 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1704 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1705 for (j=ys; j<ys+ny; j++) { 1706 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1707 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1708 for (k=zs; k<zs+nz; k++) { 1709 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1710 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1711 1712 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1713 1714 /* Find block columns in block row */ 1715 cnt = 0; 1716 for (ii=istart; ii<iend+1; ii++) { 1717 for (jj=jstart; jj<jend+1; jj++) { 1718 for (kk=kstart; kk<kend+1; kk++) { 1719 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1720 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1721 } 1722 } 1723 } 1724 } 1725 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1726 ierr = MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1727 } 1728 } 1729 } 1730 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1731 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1732 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1733 1734 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1735 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1736 1737 /* 1738 For each node in the grid: we get the neighbors in the local (on processor ordering 1739 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1740 PETSc ordering. 1741 */ 1742 if (!da->prealloc_only) { 1743 ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1744 ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1745 for (i=xs; i<xs+nx; i++) { 1746 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1747 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1748 for (j=ys; j<ys+ny; j++) { 1749 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1750 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1751 for (k=zs; k<zs+nz; k++) { 1752 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1753 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1754 1755 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1756 1757 cnt = 0; 1758 for (ii=istart; ii<iend+1; ii++) { 1759 for (jj=jstart; jj<jend+1; jj++) { 1760 for (kk=kstart; kk<kend+1; kk++) { 1761 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1762 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1763 } 1764 } 1765 } 1766 } 1767 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1768 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1769 } 1770 } 1771 } 1772 ierr = PetscFree(values);CHKERRQ(ierr); 1773 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1774 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1775 } 1776 ierr = PetscFree(cols);CHKERRQ(ierr); 1777 PetscFunctionReturn(0); 1778 } 1779 1780 /* ---------------------------------------------------------------------------------*/ 1781 1782 #undef __FUNCT__ 1783 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ_Fill" 1784 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 1785 { 1786 PetscErrorCode ierr; 1787 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1788 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz; 1789 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1790 DM_DA *dd = (DM_DA*)da->data; 1791 PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 1792 MPI_Comm comm; 1793 PetscScalar *values; 1794 DMDABoundaryType bx,by,bz; 1795 ISLocalToGlobalMapping ltog,ltogb; 1796 DMDAStencilType st; 1797 1798 PetscFunctionBegin; 1799 /* 1800 nc - number of components per grid point 1801 col - number of colors needed in one direction for single component problem 1802 1803 */ 1804 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1805 col = 2*s + 1; 1806 if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\ 1807 by 2*stencil_width + 1\n"); 1808 if (by == DMDA_BOUNDARY_PERIODIC && (n % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\ 1809 by 2*stencil_width + 1\n"); 1810 if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\ 1811 by 2*stencil_width + 1\n"); 1812 1813 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1814 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1815 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1816 1817 ierr = PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1818 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1819 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1820 1821 /* determine the matrix preallocation information */ 1822 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1823 1824 1825 for (i=xs; i<xs+nx; i++) { 1826 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1827 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1828 for (j=ys; j<ys+ny; j++) { 1829 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1830 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1831 for (k=zs; k<zs+nz; k++) { 1832 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1833 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1834 1835 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1836 1837 for (l=0; l<nc; l++) { 1838 cnt = 0; 1839 for (ii=istart; ii<iend+1; ii++) { 1840 for (jj=jstart; jj<jend+1; jj++) { 1841 for (kk=kstart; kk<kend+1; kk++) { 1842 if (ii || jj || kk) { 1843 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1844 for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) 1845 cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1846 } 1847 } else { 1848 if (dfill) { 1849 for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) 1850 cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1851 } else { 1852 for (ifill_col=0; ifill_col<nc; ifill_col++) 1853 cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1854 } 1855 } 1856 } 1857 } 1858 } 1859 row = l + nc*(slot); 1860 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1861 } 1862 } 1863 } 1864 } 1865 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1866 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1867 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1868 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1869 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1870 1871 /* 1872 For each node in the grid: we get the neighbors in the local (on processor ordering 1873 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1874 PETSc ordering. 1875 */ 1876 if (!da->prealloc_only) { 1877 ierr = PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1878 ierr = PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1879 for (i=xs; i<xs+nx; i++) { 1880 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1881 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1882 for (j=ys; j<ys+ny; j++) { 1883 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1884 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1885 for (k=zs; k<zs+nz; k++) { 1886 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1887 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1888 1889 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1890 1891 for (l=0; l<nc; l++) { 1892 cnt = 0; 1893 for (ii=istart; ii<iend+1; ii++) { 1894 for (jj=jstart; jj<jend+1; jj++) { 1895 for (kk=kstart; kk<kend+1; kk++) { 1896 if (ii || jj || kk) { 1897 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1898 for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) 1899 cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1900 } 1901 } else { 1902 if (dfill) { 1903 for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) 1904 cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1905 } else { 1906 for (ifill_col=0; ifill_col<nc; ifill_col++) 1907 cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1908 } 1909 } 1910 } 1911 } 1912 } 1913 row = l + nc*(slot); 1914 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1915 } 1916 } 1917 } 1918 } 1919 ierr = PetscFree(values);CHKERRQ(ierr); 1920 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1921 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1922 } 1923 ierr = PetscFree(cols);CHKERRQ(ierr); 1924 PetscFunctionReturn(0); 1925 } 1926