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__ "MatSetupDM" 475 /*@C 476 MatSetupDM - 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 MatSetupDM(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,"MatSetupDM_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 = MatGetDM(A, &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 = MatGetDM(A, &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 M = dd->M; 687 N = dd->N; 688 P = dd->P; 689 dim = dd->dim; 690 dof = dd->w; 691 /* ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); */ 692 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 693 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 694 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 695 ierr = MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);CHKERRQ(ierr); 696 ierr = MatSetType(A,(MatType)ttype);CHKERRQ(ierr); 697 ierr = MatSetDM(A,da);CHKERRQ(ierr); 698 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 699 ierr = MatGetType(A,&Atype);CHKERRQ(ierr); 700 /* 701 We do not provide a getmatrix function in the DMDA operations because 702 the basic DMDA does not know about matrices. We think of DMDA as being more 703 more low-level than matrices. This is kind of cheating but, cause sometimes 704 we think of DMDA has higher level than matrices. 705 706 We could switch based on Atype (or mtype), but we do not since the 707 specialized setting routines depend only the particular preallocation 708 details of the matrix, not the type itself. 709 */ 710 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 711 if (!aij) { 712 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 713 } 714 if (!aij) { 715 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 716 if (!baij) { 717 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 718 } 719 if (!baij){ 720 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 721 if (!sbaij) { 722 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 723 } 724 } 725 } 726 if (aij) { 727 if (dim == 1) { 728 ierr = DMCreateMatrix_DA_1d_MPIAIJ(da,A);CHKERRQ(ierr); 729 } else if (dim == 2) { 730 if (dd->ofill) { 731 ierr = DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 732 } else { 733 ierr = DMCreateMatrix_DA_2d_MPIAIJ(da,A);CHKERRQ(ierr); 734 } 735 } else if (dim == 3) { 736 if (dd->ofill) { 737 ierr = DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 738 } else { 739 ierr = DMCreateMatrix_DA_3d_MPIAIJ(da,A);CHKERRQ(ierr); 740 } 741 } 742 } else if (baij) { 743 if (dim == 2) { 744 ierr = DMCreateMatrix_DA_2d_MPIBAIJ(da,A);CHKERRQ(ierr); 745 } else if (dim == 3) { 746 ierr = DMCreateMatrix_DA_3d_MPIBAIJ(da,A);CHKERRQ(ierr); 747 } 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); 748 } else if (sbaij) { 749 if (dim == 2) { 750 ierr = DMCreateMatrix_DA_2d_MPISBAIJ(da,A);CHKERRQ(ierr); 751 } else if (dim == 3) { 752 ierr = DMCreateMatrix_DA_3d_MPISBAIJ(da,A);CHKERRQ(ierr); 753 } 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); 754 } else { 755 ISLocalToGlobalMapping ltog,ltogb; 756 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 757 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 758 ierr = MatSetUp(A);CHKERRQ(ierr); 759 ierr = MatSetLocalToGlobalMapping(A,ltog,ltog);CHKERRQ(ierr); 760 ierr = MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);CHKERRQ(ierr); 761 } 762 ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); 763 ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); 764 ierr = MatSetDM(A,da);CHKERRQ(ierr); 765 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 766 if (size > 1) { 767 /* change viewer to display matrix in natural ordering */ 768 ierr = MatShellSetOperation(A, MATOP_VIEW, (void (*)(void)) MatView_MPI_DA);CHKERRQ(ierr); 769 ierr = MatShellSetOperation(A, MATOP_LOAD, (void (*)(void)) MatLoad_MPI_DA);CHKERRQ(ierr); 770 } 771 *J = A; 772 PetscFunctionReturn(0); 773 } 774 775 /* ---------------------------------------------------------------------------------*/ 776 #undef __FUNCT__ 777 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ" 778 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J) 779 { 780 PetscErrorCode ierr; 781 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; 782 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 783 MPI_Comm comm; 784 PetscScalar *values; 785 DMDABoundaryType bx,by; 786 ISLocalToGlobalMapping ltog,ltogb; 787 DMDAStencilType st; 788 789 PetscFunctionBegin; 790 /* 791 nc - number of components per grid point 792 col - number of colors needed in one direction for single component problem 793 794 */ 795 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 796 col = 2*s + 1; 797 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 798 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 799 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 800 801 ierr = PetscMalloc2(nc,PetscInt,&rows,col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 802 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 803 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 804 805 /* determine the matrix preallocation information */ 806 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 807 for (i=xs; i<xs+nx; i++) { 808 809 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 810 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 811 812 for (j=ys; j<ys+ny; j++) { 813 slot = i - gxs + gnx*(j - gys); 814 815 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 816 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 817 818 cnt = 0; 819 for (k=0; k<nc; k++) { 820 for (l=lstart; l<lend+1; l++) { 821 for (p=pstart; p<pend+1; p++) { 822 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 823 cols[cnt++] = k + nc*(slot + gnx*l + p); 824 } 825 } 826 } 827 rows[k] = k + nc*(slot); 828 } 829 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 830 } 831 } 832 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 833 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 834 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 835 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 836 837 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 838 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 839 840 /* 841 For each node in the grid: we get the neighbors in the local (on processor ordering 842 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 843 PETSc ordering. 844 */ 845 if (!da->prealloc_only) { 846 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 847 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 848 for (i=xs; i<xs+nx; i++) { 849 850 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 851 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 852 853 for (j=ys; j<ys+ny; j++) { 854 slot = i - gxs + gnx*(j - gys); 855 856 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 857 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 858 859 cnt = 0; 860 for (k=0; k<nc; k++) { 861 for (l=lstart; l<lend+1; l++) { 862 for (p=pstart; p<pend+1; p++) { 863 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 864 cols[cnt++] = k + nc*(slot + gnx*l + p); 865 } 866 } 867 } 868 rows[k] = k + nc*(slot); 869 } 870 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 871 } 872 } 873 ierr = PetscFree(values);CHKERRQ(ierr); 874 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 875 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 876 } 877 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881 #undef __FUNCT__ 882 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ_Fill" 883 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J) 884 { 885 PetscErrorCode ierr; 886 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 887 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p; 888 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 889 DM_DA *dd = (DM_DA*)da->data; 890 PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill; 891 MPI_Comm comm; 892 PetscScalar *values; 893 DMDABoundaryType bx,by; 894 ISLocalToGlobalMapping ltog,ltogb; 895 DMDAStencilType st; 896 897 PetscFunctionBegin; 898 /* 899 nc - number of components per grid point 900 col - number of colors needed in one direction for single component problem 901 902 */ 903 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 904 col = 2*s + 1; 905 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 906 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 907 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 908 909 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 910 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 911 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 912 913 /* determine the matrix preallocation information */ 914 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 915 for (i=xs; i<xs+nx; i++) { 916 917 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 918 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 919 920 for (j=ys; j<ys+ny; j++) { 921 slot = i - gxs + gnx*(j - gys); 922 923 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 924 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 925 926 for (k=0; k<nc; k++) { 927 cnt = 0; 928 for (l=lstart; l<lend+1; l++) { 929 for (p=pstart; p<pend+1; p++) { 930 if (l || p) { 931 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 932 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) 933 cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 934 } 935 } else { 936 if (dfill) { 937 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) 938 cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 939 } else { 940 for (ifill_col=0; ifill_col<nc; ifill_col++) 941 cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 942 } 943 } 944 } 945 } 946 row = k + nc*(slot); 947 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 948 } 949 } 950 } 951 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 952 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 953 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 954 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 955 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 956 957 /* 958 For each node in the grid: we get the neighbors in the local (on processor ordering 959 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 960 PETSc ordering. 961 */ 962 if (!da->prealloc_only) { 963 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 964 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 965 for (i=xs; i<xs+nx; i++) { 966 967 pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 968 pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 969 970 for (j=ys; j<ys+ny; j++) { 971 slot = i - gxs + gnx*(j - gys); 972 973 lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 974 lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 975 976 for (k=0; k<nc; k++) { 977 cnt = 0; 978 for (l=lstart; l<lend+1; l++) { 979 for (p=pstart; p<pend+1; p++) { 980 if (l || p) { 981 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 982 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) 983 cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 984 } 985 } else { 986 if (dfill) { 987 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) 988 cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 989 } else { 990 for (ifill_col=0; ifill_col<nc; ifill_col++) 991 cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 992 } 993 } 994 } 995 } 996 row = k + nc*(slot); 997 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 998 } 999 } 1000 } 1001 ierr = PetscFree(values);CHKERRQ(ierr); 1002 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1003 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1004 } 1005 ierr = PetscFree(cols);CHKERRQ(ierr); 1006 PetscFunctionReturn(0); 1007 } 1008 1009 /* ---------------------------------------------------------------------------------*/ 1010 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ" 1013 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 1014 { 1015 PetscErrorCode ierr; 1016 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1017 PetscInt m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p,*dnz = PETSC_NULL,*onz = PETSC_NULL; 1018 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1019 MPI_Comm comm; 1020 PetscScalar *values; 1021 DMDABoundaryType bx,by,bz; 1022 ISLocalToGlobalMapping ltog,ltogb; 1023 DMDAStencilType st; 1024 1025 PetscFunctionBegin; 1026 /* 1027 nc - number of components per grid point 1028 col - number of colors needed in one direction for single component problem 1029 1030 */ 1031 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1032 col = 2*s + 1; 1033 1034 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1035 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1036 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1037 1038 ierr = PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 1039 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1040 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1041 1042 /* determine the matrix preallocation information */ 1043 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1044 for (i=xs; i<xs+nx; i++) { 1045 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1046 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1047 for (j=ys; j<ys+ny; j++) { 1048 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1049 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1050 for (k=zs; k<zs+nz; k++) { 1051 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1052 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1053 1054 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1055 1056 cnt = 0; 1057 for (l=0; l<nc; l++) { 1058 for (ii=istart; ii<iend+1; ii++) { 1059 for (jj=jstart; jj<jend+1; jj++) { 1060 for (kk=kstart; kk<kend+1; kk++) { 1061 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1062 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1063 } 1064 } 1065 } 1066 } 1067 rows[l] = l + nc*(slot); 1068 } 1069 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1070 } 1071 } 1072 } 1073 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1074 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1075 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1076 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1077 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1078 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1079 1080 /* 1081 For each node in the grid: we get the neighbors in the local (on processor ordering 1082 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1083 PETSc ordering. 1084 */ 1085 if (!da->prealloc_only) { 1086 ierr = PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1087 ierr = PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1088 for (i=xs; i<xs+nx; i++) { 1089 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1090 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1091 for (j=ys; j<ys+ny; j++) { 1092 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1093 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1094 for (k=zs; k<zs+nz; k++) { 1095 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1096 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1097 1098 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1099 1100 cnt = 0; 1101 for (l=0; l<nc; l++) { 1102 for (ii=istart; ii<iend+1; ii++) { 1103 for (jj=jstart; jj<jend+1; jj++) { 1104 for (kk=kstart; kk<kend+1; kk++) { 1105 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1106 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1107 } 1108 } 1109 } 1110 } 1111 rows[l] = l + nc*(slot); 1112 } 1113 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1114 } 1115 } 1116 } 1117 ierr = PetscFree(values);CHKERRQ(ierr); 1118 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1119 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1120 } 1121 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 /* ---------------------------------------------------------------------------------*/ 1126 1127 #undef __FUNCT__ 1128 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ" 1129 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 1130 { 1131 PetscErrorCode ierr; 1132 PetscInt xs,nx,i,i1,slot,gxs,gnx; 1133 PetscInt m,dim,s,*cols = PETSC_NULL,nc,*rows = PETSC_NULL,col,cnt,l; 1134 PetscInt istart,iend; 1135 PetscScalar *values; 1136 DMDABoundaryType bx; 1137 ISLocalToGlobalMapping ltog,ltogb; 1138 1139 PetscFunctionBegin; 1140 /* 1141 nc - number of components per grid point 1142 col - number of colors needed in one direction for single component problem 1143 1144 */ 1145 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1146 col = 2*s + 1; 1147 1148 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1149 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1150 1151 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1152 ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 1153 ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 1154 ierr = PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); 1155 1156 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1157 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1158 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1159 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1160 1161 /* 1162 For each node in the grid: we get the neighbors in the local (on processor ordering 1163 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1164 PETSc ordering. 1165 */ 1166 if (!da->prealloc_only) { 1167 ierr = PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1168 ierr = PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1169 for (i=xs; i<xs+nx; i++) { 1170 istart = PetscMax(-s,gxs - i); 1171 iend = PetscMin(s,gxs + gnx - i - 1); 1172 slot = i - gxs; 1173 1174 cnt = 0; 1175 for (l=0; l<nc; l++) { 1176 for (i1=istart; i1<iend+1; i1++) { 1177 cols[cnt++] = l + nc*(slot + i1); 1178 } 1179 rows[l] = l + nc*(slot); 1180 } 1181 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1182 } 1183 ierr = PetscFree(values);CHKERRQ(ierr); 1184 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1185 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1186 } 1187 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1188 PetscFunctionReturn(0); 1189 } 1190 1191 #undef __FUNCT__ 1192 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIBAIJ" 1193 PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 1194 { 1195 PetscErrorCode ierr; 1196 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1197 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1198 PetscInt istart,iend,jstart,jend,ii,jj; 1199 MPI_Comm comm; 1200 PetscScalar *values; 1201 DMDABoundaryType bx,by; 1202 DMDAStencilType st; 1203 ISLocalToGlobalMapping ltog,ltogb; 1204 1205 PetscFunctionBegin; 1206 /* 1207 nc - number of components per grid point 1208 col - number of colors needed in one direction for single component problem 1209 */ 1210 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1211 col = 2*s + 1; 1212 1213 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1214 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1215 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1216 1217 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1218 1219 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1220 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1221 1222 /* determine the matrix preallocation information */ 1223 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1224 for (i=xs; i<xs+nx; i++) { 1225 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1226 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1227 for (j=ys; j<ys+ny; j++) { 1228 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1229 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1230 slot = i - gxs + gnx*(j - gys); 1231 1232 /* Find block columns in block row */ 1233 cnt = 0; 1234 for (ii=istart; ii<iend+1; ii++) { 1235 for (jj=jstart; jj<jend+1; jj++) { 1236 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1237 cols[cnt++] = slot + ii + gnx*jj; 1238 } 1239 } 1240 } 1241 ierr = MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);CHKERRQ(ierr); 1242 } 1243 } 1244 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1245 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1246 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1247 1248 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1249 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1250 1251 /* 1252 For each node in the grid: we get the neighbors in the local (on processor ordering 1253 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1254 PETSc ordering. 1255 */ 1256 if (!da->prealloc_only) { 1257 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1258 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1259 for (i=xs; i<xs+nx; i++) { 1260 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1261 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1262 for (j=ys; j<ys+ny; j++) { 1263 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1264 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1265 slot = i - gxs + gnx*(j - gys); 1266 cnt = 0; 1267 for (ii=istart; ii<iend+1; ii++) { 1268 for (jj=jstart; jj<jend+1; jj++) { 1269 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1270 cols[cnt++] = slot + ii + gnx*jj; 1271 } 1272 } 1273 } 1274 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1275 } 1276 } 1277 ierr = PetscFree(values);CHKERRQ(ierr); 1278 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1279 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1280 } 1281 ierr = PetscFree(cols);CHKERRQ(ierr); 1282 PetscFunctionReturn(0); 1283 } 1284 1285 #undef __FUNCT__ 1286 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIBAIJ" 1287 PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 1288 { 1289 PetscErrorCode ierr; 1290 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1291 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1292 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1293 MPI_Comm comm; 1294 PetscScalar *values; 1295 DMDABoundaryType bx,by,bz; 1296 DMDAStencilType st; 1297 ISLocalToGlobalMapping ltog,ltogb; 1298 1299 PetscFunctionBegin; 1300 /* 1301 nc - number of components per grid point 1302 col - number of colors needed in one direction for single component problem 1303 1304 */ 1305 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1306 col = 2*s + 1; 1307 1308 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1309 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1310 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1311 1312 ierr = PetscMalloc(col*col*col*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1313 1314 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1315 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1316 1317 /* determine the matrix preallocation information */ 1318 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1319 for (i=xs; i<xs+nx; i++) { 1320 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1321 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1322 for (j=ys; j<ys+ny; j++) { 1323 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1324 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1325 for (k=zs; k<zs+nz; k++) { 1326 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1327 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1328 1329 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1330 1331 /* Find block columns in block row */ 1332 cnt = 0; 1333 for (ii=istart; ii<iend+1; ii++) { 1334 for (jj=jstart; jj<jend+1; jj++) { 1335 for (kk=kstart; kk<kend+1; kk++) { 1336 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1337 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1338 } 1339 } 1340 } 1341 } 1342 ierr = MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);CHKERRQ(ierr); 1343 } 1344 } 1345 } 1346 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1347 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1348 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1349 1350 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1351 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1352 1353 /* 1354 For each node in the grid: we get the neighbors in the local (on processor ordering 1355 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1356 PETSc ordering. 1357 */ 1358 if (!da->prealloc_only) { 1359 ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1360 ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1361 for (i=xs; i<xs+nx; i++) { 1362 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1363 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1364 for (j=ys; j<ys+ny; j++) { 1365 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1366 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1367 for (k=zs; k<zs+nz; k++) { 1368 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1369 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1370 1371 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1372 1373 cnt = 0; 1374 for (ii=istart; ii<iend+1; ii++) { 1375 for (jj=jstart; jj<jend+1; jj++) { 1376 for (kk=kstart; kk<kend+1; kk++) { 1377 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1378 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1379 } 1380 } 1381 } 1382 } 1383 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1384 } 1385 } 1386 } 1387 ierr = PetscFree(values);CHKERRQ(ierr); 1388 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1389 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1390 } 1391 ierr = PetscFree(cols);CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "L2GFilterUpperTriangular" 1397 /* 1398 This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 1399 identify in the local ordering with periodic domain. 1400 */ 1401 static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 1402 { 1403 PetscErrorCode ierr; 1404 PetscInt i,n; 1405 1406 PetscFunctionBegin; 1407 ierr = ISLocalToGlobalMappingApply(ltog,1,row,row);CHKERRQ(ierr); 1408 ierr = ISLocalToGlobalMappingApply(ltog,*cnt,col,col);CHKERRQ(ierr); 1409 for (i=0,n=0; i<*cnt; i++) { 1410 if (col[i] >= *row) col[n++] = col[i]; 1411 } 1412 *cnt = n; 1413 PetscFunctionReturn(0); 1414 } 1415 1416 #undef __FUNCT__ 1417 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPISBAIJ" 1418 PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 1419 { 1420 PetscErrorCode ierr; 1421 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1422 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1423 PetscInt istart,iend,jstart,jend,ii,jj; 1424 MPI_Comm comm; 1425 PetscScalar *values; 1426 DMDABoundaryType bx,by; 1427 DMDAStencilType st; 1428 ISLocalToGlobalMapping ltog,ltogb; 1429 1430 PetscFunctionBegin; 1431 /* 1432 nc - number of components per grid point 1433 col - number of colors needed in one direction for single component problem 1434 */ 1435 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1436 col = 2*s + 1; 1437 1438 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1439 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1440 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1441 1442 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1443 1444 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1445 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1446 1447 /* determine the matrix preallocation information */ 1448 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1449 for (i=xs; i<xs+nx; i++) { 1450 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1451 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1452 for (j=ys; j<ys+ny; j++) { 1453 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1454 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1455 slot = i - gxs + gnx*(j - gys); 1456 1457 /* Find block columns in block row */ 1458 cnt = 0; 1459 for (ii=istart; ii<iend+1; ii++) { 1460 for (jj=jstart; jj<jend+1; jj++) { 1461 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1462 cols[cnt++] = slot + ii + gnx*jj; 1463 } 1464 } 1465 } 1466 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1467 ierr = MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1468 } 1469 } 1470 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1471 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1472 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1473 1474 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1475 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1476 1477 /* 1478 For each node in the grid: we get the neighbors in the local (on processor ordering 1479 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1480 PETSc ordering. 1481 */ 1482 if (!da->prealloc_only) { 1483 ierr = PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1484 ierr = PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1485 for (i=xs; i<xs+nx; i++) { 1486 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1487 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1488 for (j=ys; j<ys+ny; j++) { 1489 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1490 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1491 slot = i - gxs + gnx*(j - gys); 1492 1493 /* Find block columns in block row */ 1494 cnt = 0; 1495 for (ii=istart; ii<iend+1; ii++) { 1496 for (jj=jstart; jj<jend+1; jj++) { 1497 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1498 cols[cnt++] = slot + ii + gnx*jj; 1499 } 1500 } 1501 } 1502 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1503 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1504 } 1505 } 1506 ierr = PetscFree(values);CHKERRQ(ierr); 1507 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1508 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1509 } 1510 ierr = PetscFree(cols);CHKERRQ(ierr); 1511 PetscFunctionReturn(0); 1512 } 1513 1514 #undef __FUNCT__ 1515 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPISBAIJ" 1516 PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 1517 { 1518 PetscErrorCode ierr; 1519 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1520 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1521 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1522 MPI_Comm comm; 1523 PetscScalar *values; 1524 DMDABoundaryType bx,by,bz; 1525 DMDAStencilType st; 1526 ISLocalToGlobalMapping ltog,ltogb; 1527 1528 PetscFunctionBegin; 1529 /* 1530 nc - number of components per grid point 1531 col - number of colors needed in one direction for single component problem 1532 */ 1533 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1534 col = 2*s + 1; 1535 1536 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1537 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1538 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1539 1540 /* create the matrix */ 1541 ierr = PetscMalloc(col*col*col*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1542 1543 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1544 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1545 1546 /* determine the matrix preallocation information */ 1547 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1548 for (i=xs; i<xs+nx; i++) { 1549 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1550 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1551 for (j=ys; j<ys+ny; j++) { 1552 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1553 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1554 for (k=zs; k<zs+nz; k++) { 1555 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1556 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1557 1558 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1559 1560 /* Find block columns in block row */ 1561 cnt = 0; 1562 for (ii=istart; ii<iend+1; ii++) { 1563 for (jj=jstart; jj<jend+1; jj++) { 1564 for (kk=kstart; kk<kend+1; kk++) { 1565 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1566 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1567 } 1568 } 1569 } 1570 } 1571 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1572 ierr = MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1573 } 1574 } 1575 } 1576 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1577 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1578 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1579 1580 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1581 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1582 1583 /* 1584 For each node in the grid: we get the neighbors in the local (on processor ordering 1585 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1586 PETSc ordering. 1587 */ 1588 if (!da->prealloc_only) { 1589 ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1590 ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1591 for (i=xs; i<xs+nx; i++) { 1592 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1593 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1594 for (j=ys; j<ys+ny; j++) { 1595 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1596 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1597 for (k=zs; k<zs+nz; k++) { 1598 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1599 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1600 1601 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1602 1603 cnt = 0; 1604 for (ii=istart; ii<iend+1; ii++) { 1605 for (jj=jstart; jj<jend+1; jj++) { 1606 for (kk=kstart; kk<kend+1; kk++) { 1607 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1608 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1609 } 1610 } 1611 } 1612 } 1613 ierr = L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);CHKERRQ(ierr); 1614 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1615 } 1616 } 1617 } 1618 ierr = PetscFree(values);CHKERRQ(ierr); 1619 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1620 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1621 } 1622 ierr = PetscFree(cols);CHKERRQ(ierr); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 /* ---------------------------------------------------------------------------------*/ 1627 1628 #undef __FUNCT__ 1629 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ_Fill" 1630 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 1631 { 1632 PetscErrorCode ierr; 1633 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1634 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz; 1635 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1636 DM_DA *dd = (DM_DA*)da->data; 1637 PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 1638 MPI_Comm comm; 1639 PetscScalar *values; 1640 DMDABoundaryType bx,by,bz; 1641 ISLocalToGlobalMapping ltog,ltogb; 1642 DMDAStencilType st; 1643 1644 PetscFunctionBegin; 1645 /* 1646 nc - number of components per grid point 1647 col - number of colors needed in one direction for single component problem 1648 1649 */ 1650 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1651 col = 2*s + 1; 1652 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\ 1653 by 2*stencil_width + 1\n"); 1654 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\ 1655 by 2*stencil_width + 1\n"); 1656 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\ 1657 by 2*stencil_width + 1\n"); 1658 1659 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1660 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1661 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1662 1663 ierr = PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1664 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1665 ierr = DMGetLocalToGlobalMappingBlock(da,<ogb);CHKERRQ(ierr); 1666 1667 /* determine the matrix preallocation information */ 1668 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1669 1670 1671 for (i=xs; i<xs+nx; i++) { 1672 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1673 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1674 for (j=ys; j<ys+ny; j++) { 1675 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1676 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1677 for (k=zs; k<zs+nz; k++) { 1678 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1679 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1680 1681 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1682 1683 for (l=0; l<nc; l++) { 1684 cnt = 0; 1685 for (ii=istart; ii<iend+1; ii++) { 1686 for (jj=jstart; jj<jend+1; jj++) { 1687 for (kk=kstart; kk<kend+1; kk++) { 1688 if (ii || jj || kk) { 1689 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1690 for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) 1691 cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1692 } 1693 } else { 1694 if (dfill) { 1695 for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) 1696 cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1697 } else { 1698 for (ifill_col=0; ifill_col<nc; ifill_col++) 1699 cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1700 } 1701 } 1702 } 1703 } 1704 } 1705 row = l + nc*(slot); 1706 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1707 } 1708 } 1709 } 1710 } 1711 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1712 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1713 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1714 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1715 ierr = MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);CHKERRQ(ierr); 1716 1717 /* 1718 For each node in the grid: we get the neighbors in the local (on processor ordering 1719 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1720 PETSc ordering. 1721 */ 1722 if (!da->prealloc_only) { 1723 ierr = PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1724 ierr = PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));CHKERRQ(ierr); 1725 for (i=xs; i<xs+nx; i++) { 1726 istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1727 iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1728 for (j=ys; j<ys+ny; j++) { 1729 jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1730 jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1731 for (k=zs; k<zs+nz; k++) { 1732 kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1733 kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1734 1735 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1736 1737 for (l=0; l<nc; l++) { 1738 cnt = 0; 1739 for (ii=istart; ii<iend+1; ii++) { 1740 for (jj=jstart; jj<jend+1; jj++) { 1741 for (kk=kstart; kk<kend+1; kk++) { 1742 if (ii || jj || kk) { 1743 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1744 for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) 1745 cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1746 } 1747 } else { 1748 if (dfill) { 1749 for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) 1750 cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1751 } else { 1752 for (ifill_col=0; ifill_col<nc; ifill_col++) 1753 cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1754 } 1755 } 1756 } 1757 } 1758 } 1759 row = l + nc*(slot); 1760 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1761 } 1762 } 1763 } 1764 } 1765 ierr = PetscFree(values);CHKERRQ(ierr); 1766 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1767 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1768 } 1769 ierr = PetscFree(cols);CHKERRQ(ierr); 1770 PetscFunctionReturn(0); 1771 } 1772