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 = 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 = 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 /* 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,maxcnt = 0,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 maxcnt = PetscMax(maxcnt,cnt); 960 ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 961 } 962 } 963 } 964 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 965 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 966 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 967 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 968 969 /* 970 For each node in the grid: we get the neighbors in the local (on processor ordering 971 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 972 PETSc ordering. 973 */ 974 if (!da->prealloc_only) { 975 ierr = PetscCalloc1(maxcnt,&values);CHKERRQ(ierr); 976 for (i=xs; i<xs+nx; i++) { 977 978 pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 979 pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 980 981 for (j=ys; j<ys+ny; j++) { 982 slot = i - gxs + gnx*(j - gys); 983 984 lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 985 lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 986 987 for (k=0; k<nc; k++) { 988 cnt = 0; 989 for (l=lstart; l<lend+1; l++) { 990 for (p=pstart; p<pend+1; p++) { 991 if (l || p) { 992 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 993 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 994 } 995 } else { 996 if (dfill) { 997 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 998 } else { 999 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 1000 } 1001 } 1002 } 1003 } 1004 row = k + nc*(slot); 1005 ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1006 } 1007 } 1008 } 1009 ierr = PetscFree(values);CHKERRQ(ierr); 1010 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1012 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1013 } 1014 ierr = PetscFree(cols);CHKERRQ(ierr); 1015 PetscFunctionReturn(0); 1016 } 1017 1018 /* ---------------------------------------------------------------------------------*/ 1019 1020 #undef __FUNCT__ 1021 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ" 1022 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 1023 { 1024 PetscErrorCode ierr; 1025 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1026 PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; 1027 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1028 MPI_Comm comm; 1029 PetscScalar *values; 1030 DMBoundaryType bx,by,bz; 1031 ISLocalToGlobalMapping ltog; 1032 DMDAStencilType st; 1033 1034 PetscFunctionBegin; 1035 /* 1036 nc - number of components per grid point 1037 col - number of colors needed in one direction for single component problem 1038 1039 */ 1040 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1041 col = 2*s + 1; 1042 1043 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1044 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1045 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1046 1047 ierr = PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);CHKERRQ(ierr); 1048 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1049 1050 /* determine the matrix preallocation information */ 1051 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1052 for (i=xs; i<xs+nx; i++) { 1053 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1054 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1055 for (j=ys; j<ys+ny; j++) { 1056 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1057 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1058 for (k=zs; k<zs+nz; k++) { 1059 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1060 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1061 1062 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1063 1064 cnt = 0; 1065 for (l=0; l<nc; l++) { 1066 for (ii=istart; ii<iend+1; ii++) { 1067 for (jj=jstart; jj<jend+1; jj++) { 1068 for (kk=kstart; kk<kend+1; kk++) { 1069 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1070 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1071 } 1072 } 1073 } 1074 } 1075 rows[l] = l + nc*(slot); 1076 } 1077 ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1078 } 1079 } 1080 } 1081 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1082 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 1083 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 1084 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1085 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1086 1087 /* 1088 For each node in the grid: we get the neighbors in the local (on processor ordering 1089 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1090 PETSc ordering. 1091 */ 1092 if (!da->prealloc_only) { 1093 ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 1094 for (i=xs; i<xs+nx; i++) { 1095 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1096 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1097 for (j=ys; j<ys+ny; j++) { 1098 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1099 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1100 for (k=zs; k<zs+nz; k++) { 1101 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1102 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1103 1104 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1105 1106 cnt = 0; 1107 for (l=0; l<nc; l++) { 1108 for (ii=istart; ii<iend+1; ii++) { 1109 for (jj=jstart; jj<jend+1; jj++) { 1110 for (kk=kstart; kk<kend+1; kk++) { 1111 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1112 cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1113 } 1114 } 1115 } 1116 } 1117 rows[l] = l + nc*(slot); 1118 } 1119 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1120 } 1121 } 1122 } 1123 ierr = PetscFree(values);CHKERRQ(ierr); 1124 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1125 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1126 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1127 } 1128 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 /* ---------------------------------------------------------------------------------*/ 1133 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ_Fill" 1136 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J) 1137 { 1138 PetscErrorCode ierr; 1139 DM_DA *dd = (DM_DA*)da->data; 1140 PetscInt xs,nx,i,j,gxs,gnx,row,k,l; 1141 PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols; 1142 PetscInt *ofill = dd->ofill,*dfill = dd->dfill; 1143 PetscScalar *values; 1144 DMBoundaryType bx; 1145 ISLocalToGlobalMapping ltog; 1146 PetscMPIInt rank,size; 1147 1148 PetscFunctionBegin; 1149 if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions"); 1150 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 1151 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 1152 1153 /* 1154 nc - number of components per grid point 1155 1156 */ 1157 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1158 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1159 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1160 1161 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1162 ierr = PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);CHKERRQ(ierr); 1163 1164 if (nx < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Need at least two grid points per process"); 1165 /* 1166 note should be smaller for first and last process with no periodic 1167 does not handle dfill 1168 */ 1169 cnt = 0; 1170 /* coupling with process to the left */ 1171 for (i=0; i<s; i++) { 1172 for (j=0; j<nc; j++) { 1173 ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j])); 1174 cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]); 1175 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1176 cnt++; 1177 } 1178 } 1179 for (i=s; i<nx-s; i++) { 1180 for (j=0; j<nc; j++) { 1181 cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]); 1182 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 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 maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1192 cnt++; 1193 } 1194 } 1195 1196 ierr = MatSeqAIJSetPreallocation(J,0,cols);CHKERRQ(ierr); 1197 ierr = MatMPIAIJSetPreallocation(J,0,cols,0,ocols);CHKERRQ(ierr); 1198 ierr = PetscFree2(cols,ocols);CHKERRQ(ierr); 1199 1200 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1201 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1202 1203 /* 1204 For each node in the grid: we get the neighbors in the local (on processor ordering 1205 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1206 PETSc ordering. 1207 */ 1208 if (!da->prealloc_only) { 1209 ierr = PetscCalloc2(maxcnt,&values,maxcnt,&cols);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 = PetscFree2(values,cols);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 } 1289 PetscFunctionReturn(0); 1290 } 1291 1292 /* ---------------------------------------------------------------------------------*/ 1293 1294 #undef __FUNCT__ 1295 #define __FUNCT__ "DMCreateMatrix_DA_1d_MPIAIJ" 1296 PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 1297 { 1298 PetscErrorCode ierr; 1299 PetscInt xs,nx,i,i1,slot,gxs,gnx; 1300 PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l; 1301 PetscInt istart,iend; 1302 PetscScalar *values; 1303 DMBoundaryType bx; 1304 ISLocalToGlobalMapping ltog; 1305 1306 PetscFunctionBegin; 1307 /* 1308 nc - number of components per grid point 1309 col - number of colors needed in one direction for single component problem 1310 1311 */ 1312 ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1313 col = 2*s + 1; 1314 1315 ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1316 ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1317 1318 ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 1319 ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 1320 ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 1321 1322 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1323 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1324 1325 /* 1326 For each node in the grid: we get the neighbors in the local (on processor ordering 1327 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1328 PETSc ordering. 1329 */ 1330 if (!da->prealloc_only) { 1331 ierr = PetscMalloc2(nc,&rows,col*nc*nc,&cols);CHKERRQ(ierr); 1332 ierr = PetscCalloc1(col*nc*nc,&values);CHKERRQ(ierr); 1333 for (i=xs; i<xs+nx; i++) { 1334 istart = PetscMax(-s,gxs - i); 1335 iend = PetscMin(s,gxs + gnx - i - 1); 1336 slot = i - gxs; 1337 1338 cnt = 0; 1339 for (l=0; l<nc; l++) { 1340 for (i1=istart; i1<iend+1; i1++) { 1341 cols[cnt++] = l + nc*(slot + i1); 1342 } 1343 rows[l] = l + nc*(slot); 1344 } 1345 ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1346 } 1347 ierr = PetscFree(values);CHKERRQ(ierr); 1348 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1349 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1350 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1351 ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1352 } 1353 PetscFunctionReturn(0); 1354 } 1355 1356 #undef __FUNCT__ 1357 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPIBAIJ" 1358 PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 1359 { 1360 PetscErrorCode ierr; 1361 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1362 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1363 PetscInt istart,iend,jstart,jend,ii,jj; 1364 MPI_Comm comm; 1365 PetscScalar *values; 1366 DMBoundaryType bx,by; 1367 DMDAStencilType st; 1368 ISLocalToGlobalMapping ltog; 1369 1370 PetscFunctionBegin; 1371 /* 1372 nc - number of components per grid point 1373 col - number of colors needed in one direction for single component problem 1374 */ 1375 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1376 col = 2*s + 1; 1377 1378 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1379 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1380 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1381 1382 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1383 1384 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1385 1386 /* determine the matrix preallocation information */ 1387 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1388 for (i=xs; i<xs+nx; i++) { 1389 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1390 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1391 for (j=ys; j<ys+ny; j++) { 1392 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1393 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1394 slot = i - gxs + gnx*(j - gys); 1395 1396 /* Find block columns in block row */ 1397 cnt = 0; 1398 for (ii=istart; ii<iend+1; ii++) { 1399 for (jj=jstart; jj<jend+1; jj++) { 1400 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1401 cols[cnt++] = slot + ii + gnx*jj; 1402 } 1403 } 1404 } 1405 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1406 } 1407 } 1408 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1409 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1410 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1411 1412 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1413 1414 /* 1415 For each node in the grid: we get the neighbors in the local (on processor ordering 1416 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1417 PETSc ordering. 1418 */ 1419 if (!da->prealloc_only) { 1420 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1421 for (i=xs; i<xs+nx; i++) { 1422 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1423 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1424 for (j=ys; j<ys+ny; j++) { 1425 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1426 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1427 slot = i - gxs + gnx*(j - gys); 1428 cnt = 0; 1429 for (ii=istart; ii<iend+1; ii++) { 1430 for (jj=jstart; jj<jend+1; jj++) { 1431 if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 1432 cols[cnt++] = slot + ii + gnx*jj; 1433 } 1434 } 1435 } 1436 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1437 } 1438 } 1439 ierr = PetscFree(values);CHKERRQ(ierr); 1440 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1441 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1442 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1443 } 1444 ierr = PetscFree(cols);CHKERRQ(ierr); 1445 PetscFunctionReturn(0); 1446 } 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIBAIJ" 1450 PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 1451 { 1452 PetscErrorCode ierr; 1453 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1454 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1455 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1456 MPI_Comm comm; 1457 PetscScalar *values; 1458 DMBoundaryType bx,by,bz; 1459 DMDAStencilType st; 1460 ISLocalToGlobalMapping ltog; 1461 1462 PetscFunctionBegin; 1463 /* 1464 nc - number of components per grid point 1465 col - number of colors needed in one direction for single component problem 1466 1467 */ 1468 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1469 col = 2*s + 1; 1470 1471 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1472 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1473 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1474 1475 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1476 1477 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1478 1479 /* determine the matrix preallocation information */ 1480 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1481 for (i=xs; i<xs+nx; i++) { 1482 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1483 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1484 for (j=ys; j<ys+ny; j++) { 1485 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1486 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1487 for (k=zs; k<zs+nz; k++) { 1488 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1489 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1490 1491 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1492 1493 /* Find block columns in block row */ 1494 cnt = 0; 1495 for (ii=istart; ii<iend+1; ii++) { 1496 for (jj=jstart; jj<jend+1; jj++) { 1497 for (kk=kstart; kk<kend+1; kk++) { 1498 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1499 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1500 } 1501 } 1502 } 1503 } 1504 ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1505 } 1506 } 1507 } 1508 ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1509 ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1510 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1511 1512 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1513 1514 /* 1515 For each node in the grid: we get the neighbors in the local (on processor ordering 1516 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1517 PETSc ordering. 1518 */ 1519 if (!da->prealloc_only) { 1520 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 1521 for (i=xs; i<xs+nx; i++) { 1522 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1523 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1524 for (j=ys; j<ys+ny; j++) { 1525 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1526 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1527 for (k=zs; k<zs+nz; k++) { 1528 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1529 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1530 1531 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1532 1533 cnt = 0; 1534 for (ii=istart; ii<iend+1; ii++) { 1535 for (jj=jstart; jj<jend+1; jj++) { 1536 for (kk=kstart; kk<kend+1; kk++) { 1537 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1538 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1539 } 1540 } 1541 } 1542 } 1543 ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1544 } 1545 } 1546 } 1547 ierr = PetscFree(values);CHKERRQ(ierr); 1548 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1549 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1550 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1551 } 1552 ierr = PetscFree(cols);CHKERRQ(ierr); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 #undef __FUNCT__ 1557 #define __FUNCT__ "L2GFilterUpperTriangular" 1558 /* 1559 This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 1560 identify in the local ordering with periodic domain. 1561 */ 1562 static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 1563 { 1564 PetscErrorCode ierr; 1565 PetscInt i,n; 1566 1567 PetscFunctionBegin; 1568 ierr = ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);CHKERRQ(ierr); 1569 ierr = ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);CHKERRQ(ierr); 1570 for (i=0,n=0; i<*cnt; i++) { 1571 if (col[i] >= *row) col[n++] = col[i]; 1572 } 1573 *cnt = n; 1574 PetscFunctionReturn(0); 1575 } 1576 1577 #undef __FUNCT__ 1578 #define __FUNCT__ "DMCreateMatrix_DA_2d_MPISBAIJ" 1579 PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 1580 { 1581 PetscErrorCode ierr; 1582 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1583 PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 1584 PetscInt istart,iend,jstart,jend,ii,jj; 1585 MPI_Comm comm; 1586 PetscScalar *values; 1587 DMBoundaryType bx,by; 1588 DMDAStencilType st; 1589 ISLocalToGlobalMapping ltog; 1590 1591 PetscFunctionBegin; 1592 /* 1593 nc - number of components per grid point 1594 col - number of colors needed in one direction for single component problem 1595 */ 1596 ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 1597 col = 2*s + 1; 1598 1599 ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1600 ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 1601 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1602 1603 ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 1604 1605 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1606 1607 /* determine the matrix preallocation information */ 1608 ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 1609 for (i=xs; i<xs+nx; i++) { 1610 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1611 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1612 for (j=ys; j<ys+ny; j++) { 1613 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1614 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1615 slot = i - gxs + gnx*(j - gys); 1616 1617 /* Find block columns in block row */ 1618 cnt = 0; 1619 for (ii=istart; ii<iend+1; ii++) { 1620 for (jj=jstart; jj<jend+1; jj++) { 1621 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1622 cols[cnt++] = slot + ii + gnx*jj; 1623 } 1624 } 1625 } 1626 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1627 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1628 } 1629 } 1630 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1631 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1632 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1633 1634 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1635 1636 /* 1637 For each node in the grid: we get the neighbors in the local (on processor ordering 1638 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1639 PETSc ordering. 1640 */ 1641 if (!da->prealloc_only) { 1642 ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 1643 for (i=xs; i<xs+nx; i++) { 1644 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1645 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1646 for (j=ys; j<ys+ny; j++) { 1647 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1648 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1649 slot = i - gxs + gnx*(j - gys); 1650 1651 /* Find block columns in block row */ 1652 cnt = 0; 1653 for (ii=istart; ii<iend+1; ii++) { 1654 for (jj=jstart; jj<jend+1; jj++) { 1655 if (st == DMDA_STENCIL_BOX || !ii || !jj) { 1656 cols[cnt++] = slot + ii + gnx*jj; 1657 } 1658 } 1659 } 1660 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1661 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1662 } 1663 } 1664 ierr = PetscFree(values);CHKERRQ(ierr); 1665 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1666 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1667 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1668 } 1669 ierr = PetscFree(cols);CHKERRQ(ierr); 1670 PetscFunctionReturn(0); 1671 } 1672 1673 #undef __FUNCT__ 1674 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPISBAIJ" 1675 PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 1676 { 1677 PetscErrorCode ierr; 1678 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1679 PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 1680 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1681 MPI_Comm comm; 1682 PetscScalar *values; 1683 DMBoundaryType bx,by,bz; 1684 DMDAStencilType st; 1685 ISLocalToGlobalMapping ltog; 1686 1687 PetscFunctionBegin; 1688 /* 1689 nc - number of components per grid point 1690 col - number of colors needed in one direction for single component problem 1691 */ 1692 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1693 col = 2*s + 1; 1694 1695 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1696 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1697 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1698 1699 /* create the matrix */ 1700 ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 1701 1702 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1703 1704 /* determine the matrix preallocation information */ 1705 ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1706 for (i=xs; i<xs+nx; i++) { 1707 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1708 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1709 for (j=ys; j<ys+ny; j++) { 1710 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1711 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1712 for (k=zs; k<zs+nz; k++) { 1713 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1714 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1715 1716 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1717 1718 /* Find block columns in block row */ 1719 cnt = 0; 1720 for (ii=istart; ii<iend+1; ii++) { 1721 for (jj=jstart; jj<jend+1; jj++) { 1722 for (kk=kstart; kk<kend+1; kk++) { 1723 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1724 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1725 } 1726 } 1727 } 1728 } 1729 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1730 ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 1731 } 1732 } 1733 } 1734 ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 1735 ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 1736 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1737 1738 ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1739 1740 /* 1741 For each node in the grid: we get the neighbors in the local (on processor ordering 1742 that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1743 PETSc ordering. 1744 */ 1745 if (!da->prealloc_only) { 1746 ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 1747 for (i=xs; i<xs+nx; i++) { 1748 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1749 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1750 for (j=ys; j<ys+ny; j++) { 1751 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1752 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1753 for (k=zs; k<zs+nz; k++) { 1754 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1755 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1756 1757 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1758 1759 cnt = 0; 1760 for (ii=istart; ii<iend+1; ii++) { 1761 for (jj=jstart; jj<jend+1; jj++) { 1762 for (kk=kstart; kk<kend+1; kk++) { 1763 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 1764 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 1765 } 1766 } 1767 } 1768 } 1769 ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1770 ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1771 } 1772 } 1773 } 1774 ierr = PetscFree(values);CHKERRQ(ierr); 1775 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1776 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1777 ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1778 } 1779 ierr = PetscFree(cols);CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 /* ---------------------------------------------------------------------------------*/ 1784 1785 #undef __FUNCT__ 1786 #define __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ_Fill" 1787 PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 1788 { 1789 PetscErrorCode ierr; 1790 PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1791 PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz; 1792 PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 1793 DM_DA *dd = (DM_DA*)da->data; 1794 PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 1795 MPI_Comm comm; 1796 PetscScalar *values; 1797 DMBoundaryType bx,by,bz; 1798 ISLocalToGlobalMapping ltog; 1799 DMDAStencilType st; 1800 1801 PetscFunctionBegin; 1802 /* 1803 nc - number of components per grid point 1804 col - number of colors needed in one direction for single component problem 1805 1806 */ 1807 ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 1808 col = 2*s + 1; 1809 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\ 1810 by 2*stencil_width + 1\n"); 1811 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\ 1812 by 2*stencil_width + 1\n"); 1813 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\ 1814 by 2*stencil_width + 1\n"); 1815 1816 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1817 ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 1818 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1819 1820 ierr = PetscMalloc1(col*col*col*nc,&cols);CHKERRQ(ierr); 1821 ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1822 1823 /* determine the matrix preallocation information */ 1824 ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 1825 1826 1827 for (i=xs; i<xs+nx; i++) { 1828 istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1829 iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 1830 for (j=ys; j<ys+ny; j++) { 1831 jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1832 jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 1833 for (k=zs; k<zs+nz; k++) { 1834 kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1835 kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 1836 1837 slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 1838 1839 for (l=0; l<nc; l++) { 1840 cnt = 0; 1841 for (ii=istart; ii<iend+1; ii++) { 1842 for (jj=jstart; jj<jend+1; jj++) { 1843 for (kk=kstart; kk<kend+1; kk++) { 1844 if (ii || jj || kk) { 1845 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 1846 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); 1847 } 1848 } else { 1849 if (dfill) { 1850 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); 1851 } else { 1852 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 1853 } 1854 } 1855 } 1856 } 1857 } 1858 row = l + nc*(slot); 1859 maxcnt = PetscMax(maxcnt,cnt); 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(maxcnt,&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