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