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