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