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