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