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