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