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