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