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