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; 776 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 777 ierr = MatSetUp(A);CHKERRQ(ierr); 778 ierr = MatSetLocalToGlobalMapping(A,ltog,ltog);CHKERRQ(ierr); 779 } 780 ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); 781 ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); 782 ierr = MatSetDM(A,da);CHKERRQ(ierr); 783 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 784 if (size > 1) { 785 /* change viewer to display matrix in natural ordering */ 786 ierr = MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);CHKERRQ(ierr); 787 ierr = MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);CHKERRQ(ierr); 788 } 789 *J = A; 790 PetscFunctionReturn(0); 791 } 792 793 /* ---------------------------------------------------------------------------------*/ 794 #undef __FUNCT__ 795 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ" 796 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J) 797 { 798 PetscErrorCode ierr; 799 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; 800 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 801 MPI_Comm comm; 802 PetscScalar *values; 803 DMBoundaryType bx,by; 804 ISLocalToGlobalMapping ltog; 805 DMDAStencilType st; 806 807 PetscFunctionBegin; 808 /* 809 nc - number of components per grid point 810 col - number of colors needed in one direction for single component problem 811 812 */ 813 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 814 col = 2*s + 1; 815 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 816 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 817 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 818 819 ierr = PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);CHKERRQ(ierr); 820 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 821 822 /* determine the matrix preallocation information */ 823 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 824 for (i=xs; i<xs+nx; i++) { 825 826 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 827 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 828 829 for (j=ys; j<ys+ny; j++) { 830 slot = i - gxs + gnx*(j - gys); 831 832 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 833 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 834 835 cnt = 0; 836 for (k=0; k<nc; k++) { 837 for (l=lstart; l<lend+1; l++) { 838 for (p=pstart; p<pend+1; p++) { 839 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 840 cols[cnt++] = k + nc*(slot + gnx*l + p); 841 } 842 } 843 } 844 rows[k] = k + nc*(slot); 845 } 846 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 847 } 848 } 849 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 850 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 851 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 852 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 853 854 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 855 856 /* 857 For each node in the grid: we get the neighbors in the local (on processor ordering 858 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 859 PETSc ordering. 860 */ 861 if (!da->prealloc_only) { 862 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 863 for (i=xs; i<xs+nx; i++) { 864 865 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 866 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 867 868 for (j=ys; j<ys+ny; j++) { 869 slot = i - gxs + gnx*(j - gys); 870 871 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 872 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 873 874 cnt = 0; 875 for (k=0; k<nc; k++) { 876 for (l=lstart; l<lend+1; l++) { 877 for (p=pstart; p<pend+1; p++) { 878 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 879 cols[cnt++] = k + nc*(slot + gnx*l + p); 880 } 881 } 882 } 883 rows[k] = k + nc*(slot); 884 } 885 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 886 } 887 } 888 ierr = PetscFree(values);CHKERRQ(ierr); 889 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 890 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 891 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 892 } 893 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 #undef __FUNCT__ 898 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIAIJ_Fill" 899 PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J) 900 { 901 PetscErrorCode ierr; 902 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 903 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p; 904 PetscInt lstart,lend,pstart,pend,*dnz,*onz; 905 DM_DA *dd = (DM_DA*)da->data; 906 PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill; 907 MPI_Comm comm; 908 PetscScalar *values; 909 DMBoundaryType bx,by; 910 ISLocalToGlobalMapping ltog; 911 DMDAStencilType st; 912 913 PetscFunctionBegin; 914 /* 915 nc - number of components per grid point 916 col - number of colors needed in one direction for single component problem 917 918 */ 919 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 920 col = 2*s + 1; 921 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 922 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 923 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 924 925 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 926 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 927 928 /* determine the matrix preallocation information */ 929 ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 930 for (i=xs; i<xs+nx; i++) { 931 932 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 933 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 934 935 for (j=ys; j<ys+ny; j++) { 936 slot = i - gxs + gnx*(j - gys); 937 938 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 939 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 940 941 for (k=0; k<nc; k++) { 942 cnt = 0; 943 for (l=lstart; l<lend+1; l++) { 944 for (p=pstart; p<pend+1; p++) { 945 if (l || p) { 946 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 947 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 948 } 949 } else { 950 if (dfill) { 951 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 952 } else { 953 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 954 } 955 } 956 } 957 } 958 row = k + nc*(slot); 959 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 960 } 961 } 962 } 963 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 964 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 965 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 966 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 967 968 /* 969 For each node in the grid: we get the neighbors in the local (on processor ordering 970 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 971 PETSc ordering. 972 */ 973 if (!da->prealloc_only) { 974 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 975 for (i=xs; i<xs+nx; i++) { 976 977 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 978 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 979 980 for (j=ys; j<ys+ny; j++) { 981 slot = i - gxs + gnx*(j - gys); 982 983 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 984 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 985 986 for (k=0; k<nc; k++) { 987 cnt = 0; 988 for (l=lstart; l<lend+1; l++) { 989 for (p=pstart; p<pend+1; p++) { 990 if (l || p) { 991 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 992 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 993 } 994 } else { 995 if (dfill) { 996 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 997 } else { 998 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 999 } 1000 } 1001 } 1002 } 1003 row = k + nc*(slot); 1004 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1005 } 1006 } 1007 } 1008 ierr = PetscFree(values);CHKERRQ(ierr); 1009 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1010 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1012 } 1013 ierr = PetscFree(cols);CHKERRQ(ierr); 1014 PetscFunctionReturn(0); 1015 } 1016 1017 /* ---------------------------------------------------------------------------------*/ 1018 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ" 1021 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 1022 { 1023 PetscErrorCode ierr; 1024 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1025 PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; 1026 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1027 MPI_Comm comm; 1028 PetscScalar *values; 1029 DMBoundaryType bx,by,bz; 1030 ISLocalToGlobalMapping ltog; 1031 DMDAStencilType st; 1032 1033 PetscFunctionBegin; 1034 /* 1035 nc - number of components per grid point 1036 col - number of colors needed in one direction for single component problem 1037 1038 */ 1039 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1040 col = 2*s + 1; 1041 1042 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1043 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1044 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1045 1046 ierr = PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);CHKERRQ(ierr); 1047 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1048 1049 /* determine the matrix preallocation information */ 1050 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1051 for (i=xs; i<xs+nx; i++) { 1052 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1053 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1054 for (j=ys; j<ys+ny; j++) { 1055 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1056 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1057 for (k=zs; k<zs+nz; k++) { 1058 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1059 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1060 1061 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1062 1063 cnt = 0; 1064 for (l=0; l<nc; l++) { 1065 for (ii=istart; ii<iend+1; ii++) { 1066 for (jj=jstart; jj<jend+1; jj++) { 1067 for (kk=kstart; kk<kend+1; kk++) { 1068 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1069 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1070 } 1071 } 1072 } 1073 } 1074 rows[l] = l + nc*(slot); 1075 } 1076 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1077 } 1078 } 1079 } 1080 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1081 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1082 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1083 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1084 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1085 1086 /* 1087 For each node in the grid: we get the neighbors in the local (on processor ordering 1088 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1089 PETSc ordering. 1090 */ 1091 if (!da->prealloc_only) { 1092 ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 1093 for (i=xs; i<xs+nx; i++) { 1094 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1095 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1096 for (j=ys; j<ys+ny; j++) { 1097 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1098 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1099 for (k=zs; k<zs+nz; k++) { 1100 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1101 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1102 1103 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1104 1105 cnt = 0; 1106 for (l=0; l<nc; l++) { 1107 for (ii=istart; ii<iend+1; ii++) { 1108 for (jj=jstart; jj<jend+1; jj++) { 1109 for (kk=kstart; kk<kend+1; kk++) { 1110 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1111 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1112 } 1113 } 1114 } 1115 } 1116 rows[l] = l + nc*(slot); 1117 } 1118 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1119 } 1120 } 1121 } 1122 ierr = PetscFree(values);CHKERRQ(ierr); 1123 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1124 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1125 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1126 } 1127 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1128 PetscFunctionReturn(0); 1129 } 1130 1131 /* ---------------------------------------------------------------------------------*/ 1132 1133 #undef __FUNCT__ 1134 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ_Fill" 1135 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J) 1136 { 1137 PetscErrorCode ierr; 1138 DM_DA *dd = (DM_DA*)da->data; 1139 PetscInt xs,nx,i,j,gxs,gnx,row,k,l; 1140 PetscInt m,dim,s,*cols = NULL,nc,col,cnt,*ocols; 1141 PetscInt *ofill = dd->ofill,*dfill = dd->dfill; 1142 PetscScalar *values; 1143 DMBoundaryType bx; 1144 ISLocalToGlobalMapping ltog; 1145 PetscMPIInt rank,size; 1146 1147 PetscFunctionBegin; 1148 if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions"); 1149 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 1150 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 1151 1152 /* 1153 nc - number of components per grid point 1154 col - number of colors needed in one direction for single component problem 1155 1156 */ 1157 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1158 col = 2*s + 1; 1159 1160 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1161 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1162 1163 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1164 ierr = PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);CHKERRQ(ierr); 1165 1166 if (nx < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Need at least two grid points per process"); 1167 /* 1168 note should be smaller for first and last process with no periodic 1169 does not handle dfill 1170 */ 1171 cnt = 0; 1172 /* coupling with process to the left */ 1173 for (i=0; i<s; i++) { 1174 for (j=0; j<nc; j++) { 1175 ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j])); 1176 cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]); 1177 cnt++; 1178 } 1179 } 1180 for (i=s; i<nx-s; i++) { 1181 for (j=0; j<nc; j++) { 1182 cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]); 1183 cnt++; 1184 } 1185 } 1186 /* coupling with process to the right */ 1187 for (i=nx-s; i<nx; i++) { 1188 for (j=0; j<nc; j++) { 1189 ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j])); 1190 cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]); 1191 cnt++; 1192 } 1193 } 1194 1195 ierr = MatSeqAIJSetPreallocation(J,0,cols);CHKERRQ(ierr); 1196 ierr = MatMPIAIJSetPreallocation(J,0,cols,0,ocols);CHKERRQ(ierr); 1197 ierr = PetscFree2(cols,ocols);CHKERRQ(ierr); 1198 1199 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1200 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1201 1202 /* 1203 For each node in the grid: we get the neighbors in the local (on processor ordering 1204 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1205 PETSc ordering. 1206 */ 1207 if (!da->prealloc_only) { 1208 ierr = PetscMalloc1(col*nc*nc,&cols);CHKERRQ(ierr); 1209 ierr = PetscCalloc1(col*nc*nc,&values);CHKERRQ(ierr); 1210 1211 row = xs*nc; 1212 /* coupling with process to the left */ 1213 for (i=xs; i<xs+s; i++) { 1214 for (j=0; j<nc; j++) { 1215 cnt = 0; 1216 if (rank) { 1217 for (l=0; l<s; l++) { 1218 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1219 } 1220 } 1221 if (dfill) { 1222 for (k=dfill[j]; k<dfill[j+1]; k++) { 1223 cols[cnt++] = i*nc + dfill[k]; 1224 } 1225 } else { 1226 for (k=0; k<nc; k++) { 1227 cols[cnt++] = i*nc + k; 1228 } 1229 } 1230 for (l=0; l<s; l++) { 1231 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1232 } 1233 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1234 row++; 1235 } 1236 } 1237 for (i=xs+s; i<xs+nx-s; i++) { 1238 for (j=0; j<nc; j++) { 1239 cnt = 0; 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 if (dfill) { 1244 for (k=dfill[j]; k<dfill[j+1]; k++) { 1245 cols[cnt++] = i*nc + dfill[k]; 1246 } 1247 } else { 1248 for (k=0; k<nc; k++) { 1249 cols[cnt++] = i*nc + k; 1250 } 1251 } 1252 for (l=0; l<s; l++) { 1253 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1254 } 1255 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1256 row++; 1257 } 1258 } 1259 /* coupling with process to the right */ 1260 for (i=xs+nx-s; i<xs+nx; i++) { 1261 for (j=0; j<nc; j++) { 1262 cnt = 0; 1263 for (l=0; l<s; l++) { 1264 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1265 } 1266 if (dfill) { 1267 for (k=dfill[j]; k<dfill[j+1]; k++) { 1268 cols[cnt++] = i*nc + dfill[k]; 1269 } 1270 } else { 1271 for (k=0; k<nc; k++) { 1272 cols[cnt++] = i*nc + k; 1273 } 1274 } 1275 if (rank < size-1) { 1276 for (l=0; l<s; l++) { 1277 for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1278 } 1279 } 1280 ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1281 row++; 1282 } 1283 } 1284 ierr = PetscFree(values);CHKERRQ(ierr); 1285 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1286 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1287 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1288 ierr = PetscFree(cols);CHKERRQ(ierr); 1289 } 1290 PetscFunctionReturn(0); 1291 } 1292 1293 /* ---------------------------------------------------------------------------------*/ 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ" 1297 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 1298 { 1299 PetscErrorCode ierr; 1300 PetscInt xs,nx,i,i1,slot,gxs,gnx; 1301 PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l; 1302 PetscInt istart,iend; 1303 PetscScalar *values; 1304 DMBoundaryType bx; 1305 ISLocalToGlobalMapping ltog; 1306 1307 PetscFunctionBegin; 1308 /* 1309 nc - number of components per grid point 1310 col - number of colors needed in one direction for single component problem 1311 1312 */ 1313 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1314 col = 2*s + 1; 1315 1316 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1317 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1318 1319 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1320 ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 1321 ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 1322 1323 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1324 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1325 1326 /* 1327 For each node in the grid: we get the neighbors in the local (on processor ordering 1328 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1329 PETSc ordering. 1330 */ 1331 if (!da->prealloc_only) { 1332 ierr = PetscMalloc2(nc,&rows,col*nc*nc,&cols);CHKERRQ(ierr); 1333 ierr = PetscCalloc1(col*nc*nc,&values);CHKERRQ(ierr); 1334 for (i=xs; i<xs+nx; i++) { 1335 istart = PetscMax(-s,gxs - i); 1336 iend = PetscMin(s,gxs + gnx - i - 1); 1337 slot = i - gxs; 1338 1339 cnt = 0; 1340 for (l=0; l<nc; l++) { 1341 for (i1=istart; i1<iend+1; i1++) { 1342 cols[cnt++] = l + nc*(slot + i1); 1343 } 1344 rows[l] = l + nc*(slot); 1345 } 1346 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1347 } 1348 ierr = PetscFree(values);CHKERRQ(ierr); 1349 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1350 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1351 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1352 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1353 } 1354 PetscFunctionReturn(0); 1355 } 1356 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIBAIJ" 1359 PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 1360 { 1361 PetscErrorCode ierr; 1362 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1363 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1364 PetscInt istart,iend,jstart,jend,ii,jj; 1365 MPI_Comm comm; 1366 PetscScalar *values; 1367 DMBoundaryType bx,by; 1368 DMDAStencilType st; 1369 ISLocalToGlobalMapping ltog; 1370 1371 PetscFunctionBegin; 1372 /* 1373 nc - number of components per grid point 1374 col - number of colors needed in one direction for single component problem 1375 */ 1376 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1377 col = 2*s + 1; 1378 1379 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1380 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1381 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1382 1383 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1384 1385 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1386 1387 /* determine the matrix preallocation information */ 1388 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1389 for (i=xs; i<xs+nx; i++) { 1390 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1391 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1392 for (j=ys; j<ys+ny; j++) { 1393 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1394 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1395 slot = i - gxs + gnx*(j - gys); 1396 1397 /* Find block columns in block row */ 1398 cnt = 0; 1399 for (ii=istart; ii<iend+1; ii++) { 1400 for (jj=jstart; jj<jend+1; jj++) { 1401 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1402 cols[cnt++] = slot + ii + gnx*jj; 1403 } 1404 } 1405 } 1406 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1407 } 1408 } 1409 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1410 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1411 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1412 1413 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1414 1415 /* 1416 For each node in the grid: we get the neighbors in the local (on processor ordering 1417 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1418 PETSc ordering. 1419 */ 1420 if (!da->prealloc_only) { 1421 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1422 for (i=xs; i<xs+nx; i++) { 1423 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1424 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1425 for (j=ys; j<ys+ny; j++) { 1426 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1427 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1428 slot = i - gxs + gnx*(j - gys); 1429 cnt = 0; 1430 for (ii=istart; ii<iend+1; ii++) { 1431 for (jj=jstart; jj<jend+1; jj++) { 1432 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1433 cols[cnt++] = slot + ii + gnx*jj; 1434 } 1435 } 1436 } 1437 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1438 } 1439 } 1440 ierr = PetscFree(values);CHKERRQ(ierr); 1441 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1442 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1443 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1444 } 1445 ierr = PetscFree(cols);CHKERRQ(ierr); 1446 PetscFunctionReturn(0); 1447 } 1448 1449 #undef __FUNCT__ 1450 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIBAIJ" 1451 PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 1452 { 1453 PetscErrorCode ierr; 1454 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1455 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1456 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1457 MPI_Comm comm; 1458 PetscScalar *values; 1459 DMBoundaryType bx,by,bz; 1460 DMDAStencilType st; 1461 ISLocalToGlobalMapping ltog; 1462 1463 PetscFunctionBegin; 1464 /* 1465 nc - number of components per grid point 1466 col - number of colors needed in one direction for single component problem 1467 1468 */ 1469 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1470 col = 2*s + 1; 1471 1472 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1473 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1474 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1475 1476 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1477 1478 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1479 1480 /* determine the matrix preallocation information */ 1481 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1482 for (i=xs; i<xs+nx; i++) { 1483 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1484 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1485 for (j=ys; j<ys+ny; j++) { 1486 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1487 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1488 for (k=zs; k<zs+nz; k++) { 1489 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1490 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1491 1492 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1493 1494 /* Find block columns in block row */ 1495 cnt = 0; 1496 for (ii=istart; ii<iend+1; ii++) { 1497 for (jj=jstart; jj<jend+1; jj++) { 1498 for (kk=kstart; kk<kend+1; kk++) { 1499 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1500 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1501 } 1502 } 1503 } 1504 } 1505 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1506 } 1507 } 1508 } 1509 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1510 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1511 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1512 1513 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1514 1515 /* 1516 For each node in the grid: we get the neighbors in the local (on processor ordering 1517 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1518 PETSc ordering. 1519 */ 1520 if (!da->prealloc_only) { 1521 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 1522 for (i=xs; i<xs+nx; i++) { 1523 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1524 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1525 for (j=ys; j<ys+ny; j++) { 1526 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1527 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1528 for (k=zs; k<zs+nz; k++) { 1529 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1530 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1531 1532 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1533 1534 cnt = 0; 1535 for (ii=istart; ii<iend+1; ii++) { 1536 for (jj=jstart; jj<jend+1; jj++) { 1537 for (kk=kstart; kk<kend+1; kk++) { 1538 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1539 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1540 } 1541 } 1542 } 1543 } 1544 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1545 } 1546 } 1547 } 1548 ierr = PetscFree(values);CHKERRQ(ierr); 1549 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1550 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1551 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1552 } 1553 ierr = PetscFree(cols);CHKERRQ(ierr); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "L2GFilterUpperTriangular" 1559 /* 1560 This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 1561 identify in the local ordering with periodic domain. 1562 */ 1563 static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 1564 { 1565 PetscErrorCode ierr; 1566 PetscInt i,n; 1567 1568 PetscFunctionBegin; 1569 ierr = ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);CHKERRQ(ierr); 1570 ierr = ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);CHKERRQ(ierr); 1571 for (i=0,n=0; i<*cnt; i++) { 1572 if (col[i] >= *row) col[n++] = col[i]; 1573 } 1574 *cnt = n; 1575 PetscFunctionReturn(0); 1576 } 1577 1578 #undef __FUNCT__ 1579 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPISBAIJ" 1580 PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 1581 { 1582 PetscErrorCode ierr; 1583 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1584 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1585 PetscInt istart,iend,jstart,jend,ii,jj; 1586 MPI_Comm comm; 1587 PetscScalar *values; 1588 DMBoundaryType bx,by; 1589 DMDAStencilType st; 1590 ISLocalToGlobalMapping ltog; 1591 1592 PetscFunctionBegin; 1593 /* 1594 nc - number of components per grid point 1595 col - number of colors needed in one direction for single component problem 1596 */ 1597 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1598 col = 2*s + 1; 1599 1600 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1601 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1602 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1603 1604 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1605 1606 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1607 1608 /* determine the matrix preallocation information */ 1609 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1610 for (i=xs; i<xs+nx; i++) { 1611 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1612 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1613 for (j=ys; j<ys+ny; j++) { 1614 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1615 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1616 slot = i - gxs + gnx*(j - gys); 1617 1618 /* Find block columns in block row */ 1619 cnt = 0; 1620 for (ii=istart; ii<iend+1; ii++) { 1621 for (jj=jstart; jj<jend+1; jj++) { 1622 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1623 cols[cnt++] = slot + ii + gnx*jj; 1624 } 1625 } 1626 } 1627 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1628 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1629 } 1630 } 1631 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1632 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1633 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1634 1635 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1636 1637 /* 1638 For each node in the grid: we get the neighbors in the local (on processor ordering 1639 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1640 PETSc ordering. 1641 */ 1642 if (!da->prealloc_only) { 1643 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1644 for (i=xs; i<xs+nx; i++) { 1645 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1646 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1647 for (j=ys; j<ys+ny; j++) { 1648 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1649 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1650 slot = i - gxs + gnx*(j - gys); 1651 1652 /* Find block columns in block row */ 1653 cnt = 0; 1654 for (ii=istart; ii<iend+1; ii++) { 1655 for (jj=jstart; jj<jend+1; jj++) { 1656 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1657 cols[cnt++] = slot + ii + gnx*jj; 1658 } 1659 } 1660 } 1661 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1662 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1663 } 1664 } 1665 ierr = PetscFree(values);CHKERRQ(ierr); 1666 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1667 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1668 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1669 } 1670 ierr = PetscFree(cols);CHKERRQ(ierr); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPISBAIJ" 1676 PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 1677 { 1678 PetscErrorCode ierr; 1679 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1680 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1681 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1682 MPI_Comm comm; 1683 PetscScalar *values; 1684 DMBoundaryType bx,by,bz; 1685 DMDAStencilType st; 1686 ISLocalToGlobalMapping ltog; 1687 1688 PetscFunctionBegin; 1689 /* 1690 nc - number of components per grid point 1691 col - number of colors needed in one direction for single component problem 1692 */ 1693 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1694 col = 2*s + 1; 1695 1696 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1697 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1698 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1699 1700 /* create the matrix */ 1701 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1702 1703 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1704 1705 /* determine the matrix preallocation information */ 1706 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1707 for (i=xs; i<xs+nx; i++) { 1708 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1709 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1710 for (j=ys; j<ys+ny; j++) { 1711 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1712 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1713 for (k=zs; k<zs+nz; k++) { 1714 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1715 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1716 1717 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1718 1719 /* Find block columns in block row */ 1720 cnt = 0; 1721 for (ii=istart; ii<iend+1; ii++) { 1722 for (jj=jstart; jj<jend+1; jj++) { 1723 for (kk=kstart; kk<kend+1; kk++) { 1724 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1725 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1726 } 1727 } 1728 } 1729 } 1730 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1731 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1732 } 1733 } 1734 } 1735 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1736 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1737 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1738 1739 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1740 1741 /* 1742 For each node in the grid: we get the neighbors in the local (on processor ordering 1743 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1744 PETSc ordering. 1745 */ 1746 if (!da->prealloc_only) { 1747 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 1748 for (i=xs; i<xs+nx; i++) { 1749 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1750 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1751 for (j=ys; j<ys+ny; j++) { 1752 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1753 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1754 for (k=zs; k<zs+nz; k++) { 1755 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1756 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1757 1758 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1759 1760 cnt = 0; 1761 for (ii=istart; ii<iend+1; ii++) { 1762 for (jj=jstart; jj<jend+1; jj++) { 1763 for (kk=kstart; kk<kend+1; kk++) { 1764 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1765 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1766 } 1767 } 1768 } 1769 } 1770 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1771 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1772 } 1773 } 1774 } 1775 ierr = PetscFree(values);CHKERRQ(ierr); 1776 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1777 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1778 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1779 } 1780 ierr = PetscFree(cols);CHKERRQ(ierr); 1781 PetscFunctionReturn(0); 1782 } 1783 1784 /* ---------------------------------------------------------------------------------*/ 1785 1786 #undef __FUNCT__ 1787 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ_Fill" 1788 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 1789 { 1790 PetscErrorCode ierr; 1791 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1792 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz; 1793 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1794 DM_DA *dd = (DM_DA*)da->data; 1795 PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 1796 MPI_Comm comm; 1797 PetscScalar *values; 1798 DMBoundaryType bx,by,bz; 1799 ISLocalToGlobalMapping ltog; 1800 DMDAStencilType st; 1801 1802 PetscFunctionBegin; 1803 /* 1804 nc - number of components per grid point 1805 col - number of colors needed in one direction for single component problem 1806 1807 */ 1808 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1809 col = 2*s + 1; 1810 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\ 1811 by 2*stencil_width + 1\n"); 1812 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\ 1813 by 2*stencil_width + 1\n"); 1814 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\ 1815 by 2*stencil_width + 1\n"); 1816 1817 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1818 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1819 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1820 1821 ierr = PetscMalloc1(col*col*col*nc,&cols);CHKERRQ(ierr); 1822 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1823 1824 /* determine the matrix preallocation information */ 1825 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1826 1827 1828 for (i=xs; i<xs+nx; i++) { 1829 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1830 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1831 for (j=ys; j<ys+ny; j++) { 1832 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1833 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1834 for (k=zs; k<zs+nz; k++) { 1835 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1836 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1837 1838 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1839 1840 for (l=0; l<nc; l++) { 1841 cnt = 0; 1842 for (ii=istart; ii<iend+1; ii++) { 1843 for (jj=jstart; jj<jend+1; jj++) { 1844 for (kk=kstart; kk<kend+1; kk++) { 1845 if (ii || jj || kk) { 1846 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1847 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); 1848 } 1849 } else { 1850 if (dfill) { 1851 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); 1852 } else { 1853 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1854 } 1855 } 1856 } 1857 } 1858 } 1859 row = l + nc*(slot); 1860 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1861 } 1862 } 1863 } 1864 } 1865 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1866 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1867 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1868 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1869 1870 /* 1871 For each node in the grid: we get the neighbors in the local (on processor ordering 1872 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1873 PETSc ordering. 1874 */ 1875 if (!da->prealloc_only) { 1876 ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 1877 for (i=xs; i<xs+nx; i++) { 1878 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1879 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1880 for (j=ys; j<ys+ny; j++) { 1881 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1882 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1883 for (k=zs; k<zs+nz; k++) { 1884 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1885 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1886 1887 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1888 1889 for (l=0; l<nc; l++) { 1890 cnt = 0; 1891 for (ii=istart; ii<iend+1; ii++) { 1892 for (jj=jstart; jj<jend+1; jj++) { 1893 for (kk=kstart; kk<kend+1; kk++) { 1894 if (ii || jj || kk) { 1895 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1896 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); 1897 } 1898 } else { 1899 if (dfill) { 1900 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); 1901 } else { 1902 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1903 } 1904 } 1905 } 1906 } 1907 } 1908 row = l + nc*(slot); 1909 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1910 } 1911 } 1912 } 1913 } 1914 ierr = PetscFree(values);CHKERRQ(ierr); 1915 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1916 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1917 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1918 } 1919 ierr = PetscFree(cols);CHKERRQ(ierr); 1920 PetscFunctionReturn(0); 1921 } 1922