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