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