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