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_MPISELL(DM,Mat); 490 extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM,Mat); 491 extern PetscErrorCode DMCreateMatrix_DA_IS(DM,Mat); 492 493 /*@C 494 MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix 495 496 Logically Collective on Mat 497 498 Input Parameters: 499 + mat - the matrix 500 - da - the da 501 502 Level: intermediate 503 504 @*/ 505 PetscErrorCode MatSetupDM(Mat mat,DM da) 506 { 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 511 PetscValidHeaderSpecific(da,DM_CLASSID,1); 512 ierr = PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));CHKERRQ(ierr); 513 PetscFunctionReturn(0); 514 } 515 516 PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer) 517 { 518 DM da; 519 PetscErrorCode ierr; 520 const char *prefix; 521 Mat Anatural; 522 AO ao; 523 PetscInt rstart,rend,*petsc,i; 524 IS is; 525 MPI_Comm comm; 526 PetscViewerFormat format; 527 528 PetscFunctionBegin; 529 /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */ 530 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 531 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 532 533 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 534 ierr = MatGetDM(A, &da);CHKERRQ(ierr); 535 if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 536 537 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 538 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 539 ierr = PetscMalloc1(rend-rstart,&petsc);CHKERRQ(ierr); 540 for (i=rstart; i<rend; i++) petsc[i-rstart] = i; 541 ierr = AOApplicationToPetsc(ao,rend-rstart,petsc);CHKERRQ(ierr); 542 ierr = ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 543 544 /* call viewer on natural ordering */ 545 ierr = MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);CHKERRQ(ierr); 546 ierr = ISDestroy(&is);CHKERRQ(ierr); 547 ierr = PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);CHKERRQ(ierr); 548 ierr = PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);CHKERRQ(ierr); 549 ierr = PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);CHKERRQ(ierr); 550 ierr = MatView(Anatural,viewer);CHKERRQ(ierr); 551 ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer) 556 { 557 DM da; 558 PetscErrorCode ierr; 559 Mat Anatural,Aapp; 560 AO ao; 561 PetscInt rstart,rend,*app,i,m,n,M,N; 562 IS is; 563 MPI_Comm comm; 564 565 PetscFunctionBegin; 566 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 567 ierr = MatGetDM(A, &da);CHKERRQ(ierr); 568 if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 569 570 /* Load the matrix in natural ordering */ 571 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Anatural);CHKERRQ(ierr); 572 ierr = MatSetType(Anatural,((PetscObject)A)->type_name);CHKERRQ(ierr); 573 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 574 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 575 ierr = MatSetSizes(Anatural,m,n,M,N);CHKERRQ(ierr); 576 ierr = MatLoad(Anatural,viewer);CHKERRQ(ierr); 577 578 /* Map natural ordering to application ordering and create IS */ 579 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 580 ierr = MatGetOwnershipRange(Anatural,&rstart,&rend);CHKERRQ(ierr); 581 ierr = PetscMalloc1(rend-rstart,&app);CHKERRQ(ierr); 582 for (i=rstart; i<rend; i++) app[i-rstart] = i; 583 ierr = AOPetscToApplication(ao,rend-rstart,app);CHKERRQ(ierr); 584 ierr = ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 585 586 /* Do permutation and replace header */ 587 ierr = MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);CHKERRQ(ierr); 588 ierr = MatHeaderReplace(A,&Aapp);CHKERRQ(ierr); 589 ierr = ISDestroy(&is);CHKERRQ(ierr); 590 ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 594 PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J) 595 { 596 PetscErrorCode ierr; 597 PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P; 598 Mat A; 599 MPI_Comm comm; 600 MatType Atype; 601 PetscSection section, sectionGlobal; 602 void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL,(*sell)(void)=NULL,(*is)(void)=NULL; 603 MatType mtype; 604 PetscMPIInt size; 605 DM_DA *dd = (DM_DA*)da->data; 606 607 PetscFunctionBegin; 608 ierr = MatInitializePackage();CHKERRQ(ierr); 609 mtype = da->mattype; 610 611 ierr = DMGetDefaultSection(da, §ion);CHKERRQ(ierr); 612 if (section) { 613 PetscInt bs = -1; 614 PetscInt localSize; 615 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric; 616 617 ierr = DMGetDefaultGlobalSection(da, §ionGlobal);CHKERRQ(ierr); 618 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);CHKERRQ(ierr); 619 ierr = MatCreate(PetscObjectComm((PetscObject)da),&A);CHKERRQ(ierr); 620 ierr = MatSetSizes(A,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 621 ierr = MatSetType(A,mtype);CHKERRQ(ierr); 622 ierr = PetscStrcmp(mtype,MATSHELL,&isShell);CHKERRQ(ierr); 623 ierr = PetscStrcmp(mtype,MATBAIJ,&isBlock);CHKERRQ(ierr); 624 ierr = PetscStrcmp(mtype,MATSEQBAIJ,&isSeqBlock);CHKERRQ(ierr); 625 ierr = PetscStrcmp(mtype,MATMPIBAIJ,&isMPIBlock);CHKERRQ(ierr); 626 ierr = PetscStrcmp(mtype,MATSBAIJ,&isSymBlock);CHKERRQ(ierr); 627 ierr = PetscStrcmp(mtype,MATSEQSBAIJ,&isSymSeqBlock);CHKERRQ(ierr); 628 ierr = PetscStrcmp(mtype,MATMPISBAIJ,&isSymMPIBlock);CHKERRQ(ierr); 629 /* Check for symmetric storage */ 630 isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock); 631 if (isSymmetric) { 632 ierr = MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr); 633 } 634 if (!isShell) { 635 PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal; 636 637 if (bs < 0) { 638 if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) { 639 PetscInt pStart, pEnd, p, dof; 640 641 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 642 for (p = pStart; p < pEnd; ++p) { 643 ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr); 644 if (dof) { 645 bs = dof; 646 break; 647 } 648 } 649 } else { 650 bs = 1; 651 } 652 /* Must have same blocksize on all procs (some might have no points) */ 653 bsLocal = bs; 654 ierr = MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));CHKERRQ(ierr); 655 } 656 ierr = PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);CHKERRQ(ierr); 657 /* ierr = DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);CHKERRQ(ierr); */ 658 ierr = PetscFree4(dnz, onz, dnzu, onzu);CHKERRQ(ierr); 659 } 660 } 661 /* 662 m 663 ------------------------------------------------------ 664 | | 665 | | 666 | ---------------------- | 667 | | | | 668 n | ny | | | 669 | | | | 670 | .--------------------- | 671 | (xs,ys) nx | 672 | . | 673 | (gxs,gys) | 674 | | 675 ----------------------------------------------------- 676 */ 677 678 /* 679 nc - number of components per grid point 680 col - number of colors needed in one direction for single component problem 681 682 */ 683 M = dd->M; 684 N = dd->N; 685 P = dd->P; 686 dim = da->dim; 687 dof = dd->w; 688 /* ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); */ 689 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 690 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 691 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 692 ierr = MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);CHKERRQ(ierr); 693 ierr = MatSetType(A,mtype);CHKERRQ(ierr); 694 ierr = MatSetDM(A,da);CHKERRQ(ierr); 695 if (da->structure_only) { 696 ierr = MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);CHKERRQ(ierr); 697 } 698 ierr = MatGetType(A,&Atype);CHKERRQ(ierr); 699 /* 700 We do not provide a getmatrix function in the DMDA operations because 701 the basic DMDA does not know about matrices. We think of DMDA as being more 702 more low-level than matrices. This is kind of cheating but, cause sometimes 703 we think of DMDA has higher level than matrices. 704 705 We could switch based on Atype (or mtype), but we do not since the 706 specialized setting routines depend only the particular preallocation 707 details of the matrix, not the type itself. 708 */ 709 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 710 if (!aij) { 711 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 712 } 713 if (!aij) { 714 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 715 if (!baij) { 716 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 717 } 718 if (!baij) { 719 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 720 if (!sbaij) { 721 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 722 } 723 if (!sbaij) { 724 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);CHKERRQ(ierr); 725 if (!sell) { 726 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);CHKERRQ(ierr); 727 } 728 } 729 if (!sell) { 730 ierr = PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);CHKERRQ(ierr); 731 } 732 } 733 } 734 if (aij) { 735 if (dim == 1) { 736 if (dd->ofill) { 737 ierr = DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 738 } else { 739 ierr = DMCreateMatrix_DA_1d_MPIAIJ(da,A);CHKERRQ(ierr); 740 } 741 } else if (dim == 2) { 742 if (dd->ofill) { 743 ierr = DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 744 } else { 745 ierr = DMCreateMatrix_DA_2d_MPIAIJ(da,A);CHKERRQ(ierr); 746 } 747 } else if (dim == 3) { 748 if (dd->ofill) { 749 ierr = DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 750 } else { 751 ierr = DMCreateMatrix_DA_3d_MPIAIJ(da,A);CHKERRQ(ierr); 752 } 753 } 754 } else if (baij) { 755 if (dim == 2) { 756 ierr = DMCreateMatrix_DA_2d_MPIBAIJ(da,A);CHKERRQ(ierr); 757 } else if (dim == 3) { 758 ierr = DMCreateMatrix_DA_3d_MPIBAIJ(da,A);CHKERRQ(ierr); 759 } 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); 760 } else if (sbaij) { 761 if (dim == 2) { 762 ierr = DMCreateMatrix_DA_2d_MPISBAIJ(da,A);CHKERRQ(ierr); 763 } else if (dim == 3) { 764 ierr = DMCreateMatrix_DA_3d_MPISBAIJ(da,A);CHKERRQ(ierr); 765 } 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); 766 } else if (sell) { 767 if (dim == 2) { 768 ierr = DMCreateMatrix_DA_2d_MPISELL(da,A);CHKERRQ(ierr); 769 } else if (dim == 3) { 770 ierr = DMCreateMatrix_DA_3d_MPISELL(da,A);CHKERRQ(ierr); 771 } 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); 772 } else if (is) { 773 ierr = DMCreateMatrix_DA_IS(da,A);CHKERRQ(ierr); 774 } else { 775 ISLocalToGlobalMapping ltog; 776 777 ierr = MatSetBlockSize(A,dof);CHKERRQ(ierr); 778 ierr = MatSetUp(A);CHKERRQ(ierr); 779 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 780 ierr = MatSetLocalToGlobalMapping(A,ltog,ltog);CHKERRQ(ierr); 781 } 782 ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); 783 ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); 784 ierr = MatSetDM(A,da);CHKERRQ(ierr); 785 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 786 if (size > 1) { 787 /* change viewer to display matrix in natural ordering */ 788 ierr = MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);CHKERRQ(ierr); 789 ierr = MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);CHKERRQ(ierr); 790 } 791 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 792 *J = A; 793 PetscFunctionReturn(0); 794 } 795 796 /* ---------------------------------------------------------------------------------*/ 797 PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J) 798 { 799 DM_DA *da = (DM_DA*)dm->data; 800 Mat lJ; 801 ISLocalToGlobalMapping ltog; 802 IS is_loc_filt, is_glob; 803 const PetscInt *e_loc; 804 PetscInt i,nel,nen,dnz,nv,dof,dim; 805 PetscErrorCode ierr; 806 807 /* The l2g map of DMDA has all ghosted nodes, and e_loc is a subset of all the local nodes (including the ghosted) 808 We need to filter the local indices that are represented through the DMDAGetElements decomposition 809 This is because the size of the local matrices in MATIS is the local size of the l2g map */ 810 PetscFunctionBegin; 811 dof = da->w; 812 dim = dm->dim; 813 ierr = DMGetLocalToGlobalMapping(dm,<og);CHKERRQ(ierr); 814 ierr = ISLocalToGlobalMappingGetSize(ltog,&nv);CHKERRQ(ierr); 815 ierr = DMDAGetElements(dm,&nel,&nen,&e_loc);CHKERRQ(ierr); /* this will throw an error if the stencil type is not DMDA_STENCIL_BOX */ 816 ierr = ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nel*nen,e_loc,PETSC_COPY_VALUES,&is_loc_filt);CHKERRQ(ierr); 817 ierr = DMDARestoreElements(dm,&nel,&nen,&e_loc);CHKERRQ(ierr); 818 ierr = ISSortRemoveDups(is_loc_filt);CHKERRQ(ierr); 819 ierr = ISLocalToGlobalMappingApplyIS(ltog,is_loc_filt,&is_glob);CHKERRQ(ierr); 820 ierr = ISLocalToGlobalMappingCreateIS(is_glob,<og);CHKERRQ(ierr); 821 ierr = ISDestroy(&is_glob);CHKERRQ(ierr); 822 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 823 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 824 /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */ 825 ierr = MatISGetLocalMat(J,&lJ);CHKERRQ(ierr); 826 ierr = ISCreateStride(PetscObjectComm((PetscObject)lJ),nv,0,1,&is_glob);CHKERRQ(ierr); 827 ierr = ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);CHKERRQ(ierr); 828 ierr = ISDestroy(&is_loc_filt);CHKERRQ(ierr); 829 ierr = ISGlobalToLocalMappingApplyIS(ltog,IS_GTOLM_MASK,is_glob,&is_loc_filt);CHKERRQ(ierr); 830 ierr = ISDestroy(&is_glob);CHKERRQ(ierr); 831 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 832 ierr = ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);CHKERRQ(ierr); 833 ierr = ISDestroy(&is_loc_filt);CHKERRQ(ierr); 834 ierr = MatSetLocalToGlobalMapping(lJ,ltog,ltog);CHKERRQ(ierr); 835 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 836 /* Preallocation (not exact) */ 837 switch (da->elementtype) { 838 case DMDA_ELEMENT_P1: 839 case DMDA_ELEMENT_Q1: 840 dnz = 1; 841 for (i=0; i<dim; i++) dnz *= 3; 842 dnz *= dof; 843 break; 844 default: 845 SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled element type %d",da->elementtype); 846 break; 847 } 848 ierr = MatSeqAIJSetPreallocation(lJ,dnz,NULL);CHKERRQ(ierr); 849 ierr = MatSeqBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);CHKERRQ(ierr); 850 ierr = MatSeqSBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);CHKERRQ(ierr); 851 ierr = MatISRestoreLocalMat(J,&lJ);CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854 855 PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J) 856 { 857 PetscErrorCode ierr; 858 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; 859 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 860 MPI_Comm comm; 861 PetscScalar *values; 862 DMBoundaryType bx,by; 863 ISLocalToGlobalMapping ltog; 864 DMDAStencilType st; 865 866 PetscFunctionBegin; 867 /* 868 nc - number of components per grid point 869 col - number of colors needed in one direction for single component problem 870 871 */ 872 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 873 col = 2*s + 1; 874 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 875 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 876 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 877 878 ierr = PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);CHKERRQ(ierr); 879 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 880 881 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 882 /* determine the matrix preallocation information */ 883 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 884 for (i=xs; i<xs+nx; i++) { 885 886 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 887 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 888 889 for (j=ys; j<ys+ny; j++) { 890 slot = i - gxs + gnx*(j - gys); 891 892 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 893 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 894 895 cnt = 0; 896 for (k=0; k<nc; k++) { 897 for (l=lstart; l<lend+1; l++) { 898 for (p=pstart; p<pend+1; p++) { 899 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 900 cols[cnt++] = k + nc*(slot + gnx*l + p); 901 } 902 } 903 } 904 rows[k] = k + nc*(slot); 905 } 906 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 907 } 908 } 909 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 910 ierr = MatSeqSELLSetPreallocation(J,0,dnz);CHKERRQ(ierr); 911 ierr = MatMPISELLSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 912 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 913 914 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 915 916 /* 917 For each node in the grid: we get the neighbors in the local (on processor ordering 918 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 919 PETSc ordering. 920 */ 921 if (!da->prealloc_only) { 922 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 923 for (i=xs; i<xs+nx; i++) { 924 925 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 926 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 927 928 for (j=ys; j<ys+ny; j++) { 929 slot = i - gxs + gnx*(j - gys); 930 931 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 932 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 933 934 cnt = 0; 935 for (k=0; k<nc; k++) { 936 for (l=lstart; l<lend+1; l++) { 937 for (p=pstart; p<pend+1; p++) { 938 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 939 cols[cnt++] = k + nc*(slot + gnx*l + p); 940 } 941 } 942 } 943 rows[k] = k + nc*(slot); 944 } 945 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 946 } 947 } 948 ierr = PetscFree(values);CHKERRQ(ierr); 949 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 950 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 951 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 952 } 953 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 954 PetscFunctionReturn(0); 955 } 956 957 PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da,Mat J) 958 { 959 PetscErrorCode ierr; 960 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 961 PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; 962 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P; 963 MPI_Comm comm; 964 PetscScalar *values; 965 DMBoundaryType bx,by,bz; 966 ISLocalToGlobalMapping ltog; 967 DMDAStencilType st; 968 969 PetscFunctionBegin; 970 /* 971 nc - number of components per grid point 972 col - number of colors needed in one direction for single component problem 973 974 */ 975 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 976 col = 2*s + 1; 977 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 978 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 979 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 980 981 ierr = PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);CHKERRQ(ierr); 982 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 983 984 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 985 /* determine the matrix preallocation information */ 986 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 987 for (i=xs; i<xs+nx; i++) { 988 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 989 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 990 for (j=ys; j<ys+ny; j++) { 991 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 992 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 993 for (k=zs; k<zs+nz; k++) { 994 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 995 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 996 997 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 998 999 cnt = 0; 1000 for (l=0; l<nc; l++) { 1001 for (ii=istart; ii<iend+1; ii++) { 1002 for (jj=jstart; jj<jend+1; jj++) { 1003 for (kk=kstart; kk<kend+1; kk++) { 1004 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1005 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1006 } 1007 } 1008 } 1009 } 1010 rows[l] = l + nc*(slot); 1011 } 1012 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1013 } 1014 } 1015 } 1016 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1017 ierr = MatSeqSELLSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1018 ierr = MatMPISELLSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1019 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1020 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1021 1022 /* 1023 For each node in the grid: we get the neighbors in the local (on processor ordering 1024 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1025 PETSc ordering. 1026 */ 1027 if (!da->prealloc_only) { 1028 ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 1029 for (i=xs; i<xs+nx; i++) { 1030 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1031 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1032 for (j=ys; j<ys+ny; j++) { 1033 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1034 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1035 for (k=zs; k<zs+nz; k++) { 1036 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1037 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1038 1039 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1040 1041 cnt = 0; 1042 for (l=0; l<nc; l++) { 1043 for (ii=istart; ii<iend+1; ii++) { 1044 for (jj=jstart; jj<jend+1; jj++) { 1045 for (kk=kstart; kk<kend+1; kk++) { 1046 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1047 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1048 } 1049 } 1050 } 1051 } 1052 rows[l] = l + nc*(slot); 1053 } 1054 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1055 } 1056 } 1057 } 1058 ierr = PetscFree(values);CHKERRQ(ierr); 1059 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1060 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1061 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1062 } 1063 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1064 PetscFunctionReturn(0); 1065 } 1066 1067 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J) 1068 { 1069 PetscErrorCode ierr; 1070 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; 1071 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 1072 MPI_Comm comm; 1073 PetscScalar *values; 1074 DMBoundaryType bx,by; 1075 ISLocalToGlobalMapping ltog; 1076 DMDAStencilType st; 1077 PetscBool removedups = PETSC_FALSE; 1078 1079 PetscFunctionBegin; 1080 /* 1081 nc - number of components per grid point 1082 col - number of colors needed in one direction for single component problem 1083 1084 */ 1085 ierr = DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1086 col = 2*s + 1; 1087 /* 1088 With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 1089 because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 1090 */ 1091 if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 1092 if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 1093 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1094 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1095 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1096 1097 ierr = PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);CHKERRQ(ierr); 1098 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1099 1100 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1101 /* determine the matrix preallocation information */ 1102 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 1103 for (i=xs; i<xs+nx; i++) { 1104 1105 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1106 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1107 1108 for (j=ys; j<ys+ny; j++) { 1109 slot = i - gxs + gnx*(j - gys); 1110 1111 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1112 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1113 1114 cnt = 0; 1115 for (k=0; k<nc; k++) { 1116 for (l=lstart; l<lend+1; l++) { 1117 for (p=pstart; p<pend+1; p++) { 1118 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 1119 cols[cnt++] = k + nc*(slot + gnx*l + p); 1120 } 1121 } 1122 } 1123 rows[k] = k + nc*(slot); 1124 } 1125 if (removedups) { 1126 ierr = MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1127 } else { 1128 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1129 } 1130 } 1131 } 1132 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1133 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1134 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1135 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1136 1137 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1138 1139 /* 1140 For each node in the grid: we get the neighbors in the local (on processor ordering 1141 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1142 PETSc ordering. 1143 */ 1144 if (!da->prealloc_only) { 1145 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1146 for (i=xs; i<xs+nx; i++) { 1147 1148 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1149 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1150 1151 for (j=ys; j<ys+ny; j++) { 1152 slot = i - gxs + gnx*(j - gys); 1153 1154 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1155 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1156 1157 cnt = 0; 1158 for (k=0; k<nc; k++) { 1159 for (l=lstart; l<lend+1; l++) { 1160 for (p=pstart; p<pend+1; p++) { 1161 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 1162 cols[cnt++] = k + nc*(slot + gnx*l + p); 1163 } 1164 } 1165 } 1166 rows[k] = k + nc*(slot); 1167 } 1168 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1169 } 1170 } 1171 ierr = PetscFree(values);CHKERRQ(ierr); 1172 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1173 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1174 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1175 } 1176 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J) 1181 { 1182 PetscErrorCode ierr; 1183 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1184 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N; 1185 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 1186 DM_DA *dd = (DM_DA*)da->data; 1187 PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill; 1188 MPI_Comm comm; 1189 PetscScalar *values; 1190 DMBoundaryType bx,by; 1191 ISLocalToGlobalMapping ltog; 1192 DMDAStencilType st; 1193 PetscBool removedups = PETSC_FALSE; 1194 1195 PetscFunctionBegin; 1196 /* 1197 nc - number of components per grid point 1198 col - number of colors needed in one direction for single component problem 1199 1200 */ 1201 ierr = DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1202 col = 2*s + 1; 1203 /* 1204 With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 1205 because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 1206 */ 1207 if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 1208 if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 1209 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1210 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1211 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1212 1213 ierr = PetscMalloc1(col*col*nc,&cols);CHKERRQ(ierr); 1214 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1215 1216 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1217 /* determine the matrix preallocation information */ 1218 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 1219 for (i=xs; i<xs+nx; i++) { 1220 1221 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1222 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1223 1224 for (j=ys; j<ys+ny; j++) { 1225 slot = i - gxs + gnx*(j - gys); 1226 1227 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1228 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1229 1230 for (k=0; k<nc; k++) { 1231 cnt = 0; 1232 for (l=lstart; l<lend+1; l++) { 1233 for (p=pstart; p<pend+1; p++) { 1234 if (l || p) { 1235 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 1236 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 1237 } 1238 } else { 1239 if (dfill) { 1240 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 1241 } else { 1242 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 1243 } 1244 } 1245 } 1246 } 1247 row = k + nc*(slot); 1248 maxcnt = PetscMax(maxcnt,cnt); 1249 if (removedups) { 1250 ierr = MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1251 } else { 1252 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1253 } 1254 } 1255 } 1256 } 1257 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1258 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1259 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1260 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1261 1262 /* 1263 For each node in the grid: we get the neighbors in the local (on processor ordering 1264 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1265 PETSc ordering. 1266 */ 1267 if (!da->prealloc_only) { 1268 ierr = PetscCalloc1(maxcnt,&values);CHKERRQ(ierr); 1269 for (i=xs; i<xs+nx; i++) { 1270 1271 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1272 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1273 1274 for (j=ys; j<ys+ny; j++) { 1275 slot = i - gxs + gnx*(j - gys); 1276 1277 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1278 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1279 1280 for (k=0; k<nc; k++) { 1281 cnt = 0; 1282 for (l=lstart; l<lend+1; l++) { 1283 for (p=pstart; p<pend+1; p++) { 1284 if (l || p) { 1285 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 1286 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 1287 } 1288 } else { 1289 if (dfill) { 1290 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 1291 } else { 1292 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 1293 } 1294 } 1295 } 1296 } 1297 row = k + nc*(slot); 1298 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1299 } 1300 } 1301 } 1302 ierr = PetscFree(values);CHKERRQ(ierr); 1303 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1304 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1305 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1306 } 1307 ierr = PetscFree(cols);CHKERRQ(ierr); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 /* ---------------------------------------------------------------------------------*/ 1312 1313 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 1314 { 1315 PetscErrorCode ierr; 1316 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1317 PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; 1318 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P; 1319 MPI_Comm comm; 1320 PetscScalar *values; 1321 DMBoundaryType bx,by,bz; 1322 ISLocalToGlobalMapping ltog; 1323 DMDAStencilType st; 1324 PetscBool removedups = PETSC_FALSE; 1325 1326 PetscFunctionBegin; 1327 /* 1328 nc - number of components per grid point 1329 col - number of colors needed in one direction for single component problem 1330 1331 */ 1332 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1333 col = 2*s + 1; 1334 1335 /* 1336 With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 1337 because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 1338 */ 1339 if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 1340 if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 1341 if (P == 1 && 2*s >= p) removedups = PETSC_TRUE; 1342 1343 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1344 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1345 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1346 1347 ierr = PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);CHKERRQ(ierr); 1348 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1349 1350 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1351 /* determine the matrix preallocation information */ 1352 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1353 for (i=xs; i<xs+nx; i++) { 1354 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1355 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1356 for (j=ys; j<ys+ny; j++) { 1357 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1358 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1359 for (k=zs; k<zs+nz; k++) { 1360 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1361 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1362 1363 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1364 1365 cnt = 0; 1366 for (l=0; l<nc; l++) { 1367 for (ii=istart; ii<iend+1; ii++) { 1368 for (jj=jstart; jj<jend+1; jj++) { 1369 for (kk=kstart; kk<kend+1; kk++) { 1370 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1371 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1372 } 1373 } 1374 } 1375 } 1376 rows[l] = l + nc*(slot); 1377 } 1378 if (removedups) { 1379 ierr = MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1380 } else { 1381 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1382 } 1383 } 1384 } 1385 } 1386 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1387 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1388 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1389 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1390 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1391 1392 /* 1393 For each node in the grid: we get the neighbors in the local (on processor ordering 1394 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1395 PETSc ordering. 1396 */ 1397 if (!da->prealloc_only) { 1398 ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 1399 for (i=xs; i<xs+nx; i++) { 1400 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1401 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1402 for (j=ys; j<ys+ny; j++) { 1403 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1404 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1405 for (k=zs; k<zs+nz; k++) { 1406 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1407 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1408 1409 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1410 1411 cnt = 0; 1412 for (l=0; l<nc; l++) { 1413 for (ii=istart; ii<iend+1; ii++) { 1414 for (jj=jstart; jj<jend+1; jj++) { 1415 for (kk=kstart; kk<kend+1; kk++) { 1416 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1417 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1418 } 1419 } 1420 } 1421 } 1422 rows[l] = l + nc*(slot); 1423 } 1424 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1425 } 1426 } 1427 } 1428 ierr = PetscFree(values);CHKERRQ(ierr); 1429 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1430 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1431 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1432 } 1433 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1434 PetscFunctionReturn(0); 1435 } 1436 1437 /* ---------------------------------------------------------------------------------*/ 1438 1439 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J) 1440 { 1441 PetscErrorCode ierr; 1442 DM_DA *dd = (DM_DA*)da->data; 1443 PetscInt xs,nx,i,j,gxs,gnx,row,k,l; 1444 PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols; 1445 PetscInt *ofill = dd->ofill,*dfill = dd->dfill; 1446 PetscScalar *values; 1447 DMBoundaryType bx; 1448 ISLocalToGlobalMapping ltog; 1449 PetscMPIInt rank,size; 1450 1451 PetscFunctionBegin; 1452 if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions"); 1453 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 1454 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 1455 1456 /* 1457 nc - number of components per grid point 1458 1459 */ 1460 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1461 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1462 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1463 1464 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1465 ierr = PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);CHKERRQ(ierr); 1466 1467 /* 1468 note should be smaller for first and last process with no periodic 1469 does not handle dfill 1470 */ 1471 cnt = 0; 1472 /* coupling with process to the left */ 1473 for (i=0; i<s; i++) { 1474 for (j=0; j<nc; j++) { 1475 ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j])); 1476 cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]); 1477 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1478 cnt++; 1479 } 1480 } 1481 for (i=s; i<nx-s; i++) { 1482 for (j=0; j<nc; j++) { 1483 cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]); 1484 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1485 cnt++; 1486 } 1487 } 1488 /* coupling with process to the right */ 1489 for (i=nx-s; i<nx; i++) { 1490 for (j=0; j<nc; j++) { 1491 ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j])); 1492 cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]); 1493 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1494 cnt++; 1495 } 1496 } 1497 1498 ierr = MatSeqAIJSetPreallocation(J,0,cols);CHKERRQ(ierr); 1499 ierr = MatMPIAIJSetPreallocation(J,0,cols,0,ocols);CHKERRQ(ierr); 1500 ierr = PetscFree2(cols,ocols);CHKERRQ(ierr); 1501 1502 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1503 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1504 1505 /* 1506 For each node in the grid: we get the neighbors in the local (on processor ordering 1507 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1508 PETSc ordering. 1509 */ 1510 if (!da->prealloc_only) { 1511 ierr = PetscCalloc2(maxcnt,&values,maxcnt,&cols);CHKERRQ(ierr); 1512 1513 row = xs*nc; 1514 /* coupling with process to the left */ 1515 for (i=xs; i<xs+s; i++) { 1516 for (j=0; j<nc; j++) { 1517 cnt = 0; 1518 if (rank) { 1519 for (l=0; l<s; l++) { 1520 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1521 } 1522 } 1523 if (dfill) { 1524 for (k=dfill[j]; k<dfill[j+1]; k++) { 1525 cols[cnt++] = i*nc + dfill[k]; 1526 } 1527 } else { 1528 for (k=0; k<nc; k++) { 1529 cols[cnt++] = i*nc + k; 1530 } 1531 } 1532 for (l=0; l<s; l++) { 1533 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1534 } 1535 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1536 row++; 1537 } 1538 } 1539 for (i=xs+s; i<xs+nx-s; i++) { 1540 for (j=0; j<nc; j++) { 1541 cnt = 0; 1542 for (l=0; l<s; l++) { 1543 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1544 } 1545 if (dfill) { 1546 for (k=dfill[j]; k<dfill[j+1]; k++) { 1547 cols[cnt++] = i*nc + dfill[k]; 1548 } 1549 } else { 1550 for (k=0; k<nc; k++) { 1551 cols[cnt++] = i*nc + k; 1552 } 1553 } 1554 for (l=0; l<s; l++) { 1555 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1556 } 1557 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1558 row++; 1559 } 1560 } 1561 /* coupling with process to the right */ 1562 for (i=xs+nx-s; i<xs+nx; i++) { 1563 for (j=0; j<nc; j++) { 1564 cnt = 0; 1565 for (l=0; l<s; l++) { 1566 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1567 } 1568 if (dfill) { 1569 for (k=dfill[j]; k<dfill[j+1]; k++) { 1570 cols[cnt++] = i*nc + dfill[k]; 1571 } 1572 } else { 1573 for (k=0; k<nc; k++) { 1574 cols[cnt++] = i*nc + k; 1575 } 1576 } 1577 if (rank < size-1) { 1578 for (l=0; l<s; l++) { 1579 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1580 } 1581 } 1582 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1583 row++; 1584 } 1585 } 1586 ierr = PetscFree2(values,cols);CHKERRQ(ierr); 1587 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1588 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1589 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1590 } 1591 PetscFunctionReturn(0); 1592 } 1593 1594 /* ---------------------------------------------------------------------------------*/ 1595 1596 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 1597 { 1598 PetscErrorCode ierr; 1599 PetscInt xs,nx,i,i1,slot,gxs,gnx; 1600 PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l; 1601 PetscInt istart,iend; 1602 PetscScalar *values; 1603 DMBoundaryType bx; 1604 ISLocalToGlobalMapping ltog; 1605 1606 PetscFunctionBegin; 1607 /* 1608 nc - number of components per grid point 1609 col - number of colors needed in one direction for single component problem 1610 1611 */ 1612 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1613 col = 2*s + 1; 1614 1615 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1616 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1617 1618 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1619 ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 1620 ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 1621 1622 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1623 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1624 1625 /* 1626 For each node in the grid: we get the neighbors in the local (on processor ordering 1627 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1628 PETSc ordering. 1629 */ 1630 if (!da->prealloc_only) { 1631 ierr = PetscMalloc2(nc,&rows,col*nc*nc,&cols);CHKERRQ(ierr); 1632 ierr = PetscCalloc1(col*nc*nc,&values);CHKERRQ(ierr); 1633 for (i=xs; i<xs+nx; i++) { 1634 istart = PetscMax(-s,gxs - i); 1635 iend = PetscMin(s,gxs + gnx - i - 1); 1636 slot = i - gxs; 1637 1638 cnt = 0; 1639 for (l=0; l<nc; l++) { 1640 for (i1=istart; i1<iend+1; i1++) { 1641 cols[cnt++] = l + nc*(slot + i1); 1642 } 1643 rows[l] = l + nc*(slot); 1644 } 1645 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1646 } 1647 ierr = PetscFree(values);CHKERRQ(ierr); 1648 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1649 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1650 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1651 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1652 } 1653 PetscFunctionReturn(0); 1654 } 1655 1656 PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 1657 { 1658 PetscErrorCode ierr; 1659 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1660 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1661 PetscInt istart,iend,jstart,jend,ii,jj; 1662 MPI_Comm comm; 1663 PetscScalar *values; 1664 DMBoundaryType bx,by; 1665 DMDAStencilType st; 1666 ISLocalToGlobalMapping ltog; 1667 1668 PetscFunctionBegin; 1669 /* 1670 nc - number of components per grid point 1671 col - number of colors needed in one direction for single component problem 1672 */ 1673 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1674 col = 2*s + 1; 1675 1676 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1677 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1678 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1679 1680 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1681 1682 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1683 1684 /* determine the matrix preallocation information */ 1685 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1686 for (i=xs; i<xs+nx; i++) { 1687 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1688 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1689 for (j=ys; j<ys+ny; j++) { 1690 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1691 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1692 slot = i - gxs + gnx*(j - gys); 1693 1694 /* Find block columns in block row */ 1695 cnt = 0; 1696 for (ii=istart; ii<iend+1; ii++) { 1697 for (jj=jstart; jj<jend+1; jj++) { 1698 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1699 cols[cnt++] = slot + ii + gnx*jj; 1700 } 1701 } 1702 } 1703 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1704 } 1705 } 1706 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1707 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1708 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1709 1710 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1711 1712 /* 1713 For each node in the grid: we get the neighbors in the local (on processor ordering 1714 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1715 PETSc ordering. 1716 */ 1717 if (!da->prealloc_only) { 1718 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1719 for (i=xs; i<xs+nx; i++) { 1720 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1721 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1722 for (j=ys; j<ys+ny; j++) { 1723 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1724 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1725 slot = i - gxs + gnx*(j - gys); 1726 cnt = 0; 1727 for (ii=istart; ii<iend+1; ii++) { 1728 for (jj=jstart; jj<jend+1; jj++) { 1729 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1730 cols[cnt++] = slot + ii + gnx*jj; 1731 } 1732 } 1733 } 1734 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1735 } 1736 } 1737 ierr = PetscFree(values);CHKERRQ(ierr); 1738 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1739 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1740 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1741 } 1742 ierr = PetscFree(cols);CHKERRQ(ierr); 1743 PetscFunctionReturn(0); 1744 } 1745 1746 PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 1747 { 1748 PetscErrorCode ierr; 1749 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1750 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1751 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1752 MPI_Comm comm; 1753 PetscScalar *values; 1754 DMBoundaryType bx,by,bz; 1755 DMDAStencilType st; 1756 ISLocalToGlobalMapping ltog; 1757 1758 PetscFunctionBegin; 1759 /* 1760 nc - number of components per grid point 1761 col - number of colors needed in one direction for single component problem 1762 1763 */ 1764 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1765 col = 2*s + 1; 1766 1767 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1768 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1769 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1770 1771 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1772 1773 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1774 1775 /* determine the matrix preallocation information */ 1776 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1777 for (i=xs; i<xs+nx; i++) { 1778 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1779 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1780 for (j=ys; j<ys+ny; j++) { 1781 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1782 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1783 for (k=zs; k<zs+nz; k++) { 1784 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1785 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1786 1787 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1788 1789 /* Find block columns in block row */ 1790 cnt = 0; 1791 for (ii=istart; ii<iend+1; ii++) { 1792 for (jj=jstart; jj<jend+1; jj++) { 1793 for (kk=kstart; kk<kend+1; kk++) { 1794 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1795 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1796 } 1797 } 1798 } 1799 } 1800 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1801 } 1802 } 1803 } 1804 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1805 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1806 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1807 1808 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1809 1810 /* 1811 For each node in the grid: we get the neighbors in the local (on processor ordering 1812 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1813 PETSc ordering. 1814 */ 1815 if (!da->prealloc_only) { 1816 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 1817 for (i=xs; i<xs+nx; i++) { 1818 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1819 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1820 for (j=ys; j<ys+ny; j++) { 1821 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1822 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1823 for (k=zs; k<zs+nz; k++) { 1824 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1825 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1826 1827 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1828 1829 cnt = 0; 1830 for (ii=istart; ii<iend+1; ii++) { 1831 for (jj=jstart; jj<jend+1; jj++) { 1832 for (kk=kstart; kk<kend+1; kk++) { 1833 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1834 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1835 } 1836 } 1837 } 1838 } 1839 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1840 } 1841 } 1842 } 1843 ierr = PetscFree(values);CHKERRQ(ierr); 1844 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1845 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1846 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1847 } 1848 ierr = PetscFree(cols);CHKERRQ(ierr); 1849 PetscFunctionReturn(0); 1850 } 1851 1852 /* 1853 This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 1854 identify in the local ordering with periodic domain. 1855 */ 1856 static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 1857 { 1858 PetscErrorCode ierr; 1859 PetscInt i,n; 1860 1861 PetscFunctionBegin; 1862 ierr = ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);CHKERRQ(ierr); 1863 ierr = ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);CHKERRQ(ierr); 1864 for (i=0,n=0; i<*cnt; i++) { 1865 if (col[i] >= *row) col[n++] = col[i]; 1866 } 1867 *cnt = n; 1868 PetscFunctionReturn(0); 1869 } 1870 1871 PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 1872 { 1873 PetscErrorCode ierr; 1874 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1875 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1876 PetscInt istart,iend,jstart,jend,ii,jj; 1877 MPI_Comm comm; 1878 PetscScalar *values; 1879 DMBoundaryType bx,by; 1880 DMDAStencilType st; 1881 ISLocalToGlobalMapping ltog; 1882 1883 PetscFunctionBegin; 1884 /* 1885 nc - number of components per grid point 1886 col - number of colors needed in one direction for single component problem 1887 */ 1888 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1889 col = 2*s + 1; 1890 1891 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1892 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1893 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1894 1895 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1896 1897 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1898 1899 /* determine the matrix preallocation information */ 1900 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1901 for (i=xs; i<xs+nx; i++) { 1902 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1903 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1904 for (j=ys; j<ys+ny; j++) { 1905 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1906 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1907 slot = i - gxs + gnx*(j - gys); 1908 1909 /* Find block columns in block row */ 1910 cnt = 0; 1911 for (ii=istart; ii<iend+1; ii++) { 1912 for (jj=jstart; jj<jend+1; jj++) { 1913 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1914 cols[cnt++] = slot + ii + gnx*jj; 1915 } 1916 } 1917 } 1918 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1919 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1920 } 1921 } 1922 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1923 ierr = MatMPISBAIJSetPreallocation(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*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 slot = i - gxs + gnx*(j - gys); 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 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1948 cols[cnt++] = slot + ii + gnx*jj; 1949 } 1950 } 1951 } 1952 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1953 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1954 } 1955 } 1956 ierr = PetscFree(values);CHKERRQ(ierr); 1957 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1958 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1959 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1960 } 1961 ierr = PetscFree(cols);CHKERRQ(ierr); 1962 PetscFunctionReturn(0); 1963 } 1964 1965 PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 1966 { 1967 PetscErrorCode ierr; 1968 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1969 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1970 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1971 MPI_Comm comm; 1972 PetscScalar *values; 1973 DMBoundaryType bx,by,bz; 1974 DMDAStencilType st; 1975 ISLocalToGlobalMapping ltog; 1976 1977 PetscFunctionBegin; 1978 /* 1979 nc - number of components per grid point 1980 col - number of colors needed in one direction for single component problem 1981 */ 1982 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1983 col = 2*s + 1; 1984 1985 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1986 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1987 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1988 1989 /* create the matrix */ 1990 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1991 1992 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1993 1994 /* determine the matrix preallocation information */ 1995 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1996 for (i=xs; i<xs+nx; i++) { 1997 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1998 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1999 for (j=ys; j<ys+ny; j++) { 2000 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 2001 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 2002 for (k=zs; k<zs+nz; k++) { 2003 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 2004 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 2005 2006 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 2007 2008 /* Find block columns in block row */ 2009 cnt = 0; 2010 for (ii=istart; ii<iend+1; ii++) { 2011 for (jj=jstart; jj<jend+1; jj++) { 2012 for (kk=kstart; kk<kend+1; kk++) { 2013 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 2014 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 2015 } 2016 } 2017 } 2018 } 2019 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 2020 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 2021 } 2022 } 2023 } 2024 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 2025 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 2026 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 2027 2028 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 2029 2030 /* 2031 For each node in the grid: we get the neighbors in the local (on processor ordering 2032 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 2033 PETSc ordering. 2034 */ 2035 if (!da->prealloc_only) { 2036 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 2037 for (i=xs; i<xs+nx; i++) { 2038 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 2039 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 2040 for (j=ys; j<ys+ny; j++) { 2041 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 2042 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 2043 for (k=zs; k<zs+nz; k++) { 2044 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 2045 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 2046 2047 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 2048 2049 cnt = 0; 2050 for (ii=istart; ii<iend+1; ii++) { 2051 for (jj=jstart; jj<jend+1; jj++) { 2052 for (kk=kstart; kk<kend+1; kk++) { 2053 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 2054 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 2055 } 2056 } 2057 } 2058 } 2059 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 2060 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 2061 } 2062 } 2063 } 2064 ierr = PetscFree(values);CHKERRQ(ierr); 2065 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2066 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2067 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2068 } 2069 ierr = PetscFree(cols);CHKERRQ(ierr); 2070 PetscFunctionReturn(0); 2071 } 2072 2073 /* ---------------------------------------------------------------------------------*/ 2074 2075 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 2076 { 2077 PetscErrorCode ierr; 2078 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 2079 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz; 2080 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P; 2081 DM_DA *dd = (DM_DA*)da->data; 2082 PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 2083 MPI_Comm comm; 2084 PetscScalar *values; 2085 DMBoundaryType bx,by,bz; 2086 ISLocalToGlobalMapping ltog; 2087 DMDAStencilType st; 2088 PetscBool removedups = PETSC_FALSE; 2089 2090 PetscFunctionBegin; 2091 /* 2092 nc - number of components per grid point 2093 col - number of colors needed in one direction for single component problem 2094 2095 */ 2096 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 2097 col = 2*s + 1; 2098 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\ 2099 by 2*stencil_width + 1\n"); 2100 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\ 2101 by 2*stencil_width + 1\n"); 2102 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\ 2103 by 2*stencil_width + 1\n"); 2104 2105 /* 2106 With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 2107 because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 2108 */ 2109 if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 2110 if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 2111 if (P == 1 && 2*s >= p) removedups = PETSC_TRUE; 2112 2113 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 2114 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 2115 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 2116 2117 ierr = PetscMalloc1(col*col*col*nc,&cols);CHKERRQ(ierr); 2118 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 2119 2120 /* determine the matrix preallocation information */ 2121 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 2122 2123 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 2124 for (i=xs; i<xs+nx; i++) { 2125 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 2126 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 2127 for (j=ys; j<ys+ny; j++) { 2128 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 2129 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 2130 for (k=zs; k<zs+nz; k++) { 2131 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 2132 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 2133 2134 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 2135 2136 for (l=0; l<nc; l++) { 2137 cnt = 0; 2138 for (ii=istart; ii<iend+1; ii++) { 2139 for (jj=jstart; jj<jend+1; jj++) { 2140 for (kk=kstart; kk<kend+1; kk++) { 2141 if (ii || jj || kk) { 2142 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 2143 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); 2144 } 2145 } else { 2146 if (dfill) { 2147 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); 2148 } else { 2149 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 2150 } 2151 } 2152 } 2153 } 2154 } 2155 row = l + nc*(slot); 2156 maxcnt = PetscMax(maxcnt,cnt); 2157 if (removedups) { 2158 ierr = MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 2159 } else { 2160 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 2161 } 2162 } 2163 } 2164 } 2165 } 2166 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 2167 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 2168 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 2169 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 2170 2171 /* 2172 For each node in the grid: we get the neighbors in the local (on processor ordering 2173 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 2174 PETSc ordering. 2175 */ 2176 if (!da->prealloc_only) { 2177 ierr = PetscCalloc1(maxcnt,&values);CHKERRQ(ierr); 2178 for (i=xs; i<xs+nx; i++) { 2179 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 2180 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 2181 for (j=ys; j<ys+ny; j++) { 2182 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 2183 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 2184 for (k=zs; k<zs+nz; k++) { 2185 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 2186 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 2187 2188 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 2189 2190 for (l=0; l<nc; l++) { 2191 cnt = 0; 2192 for (ii=istart; ii<iend+1; ii++) { 2193 for (jj=jstart; jj<jend+1; jj++) { 2194 for (kk=kstart; kk<kend+1; kk++) { 2195 if (ii || jj || kk) { 2196 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 2197 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); 2198 } 2199 } else { 2200 if (dfill) { 2201 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); 2202 } else { 2203 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 2204 } 2205 } 2206 } 2207 } 2208 } 2209 row = l + nc*(slot); 2210 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 2211 } 2212 } 2213 } 2214 } 2215 ierr = PetscFree(values);CHKERRQ(ierr); 2216 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2217 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2218 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2219 } 2220 ierr = PetscFree(cols);CHKERRQ(ierr); 2221 PetscFunctionReturn(0); 2222 } 2223