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