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