Actual source code: fdda.c

petsc-3.10.3 2018-12-18
Report Typos and Errors

  2:  #include <petsc/private/dmdaimpl.h>
  3:  #include <petscmat.h>

  5: extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring*);
  6: extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring*);
  7: extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring*);
  8: extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring*);

 10: /*
 11:    For ghost i that may be negative or greater than the upper bound this
 12:   maps it into the 0:m-1 range using periodicity
 13: */
 14: #define SetInRange(i,m) ((i < 0) ? m+i : ((i >= m) ? i-m : i))

 16: static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill)
 17: {
 19:   PetscInt       i,j,nz,*fill;

 22:   if (!dfill) return(0);

 24:   /* count number nonzeros */
 25:   nz = 0;
 26:   for (i=0; i<w; i++) {
 27:     for (j=0; j<w; j++) {
 28:       if (dfill[w*i+j]) nz++;
 29:     }
 30:   }
 31:   PetscMalloc1(nz + w + 1,&fill);
 32:   /* construct modified CSR storage of nonzero structure */
 33:   /*  fill[0 -- w] marks starts of each row of column indices (and end of last row)
 34:    so fill[1] - fill[0] gives number of nonzeros in first row etc */
 35:   nz = w + 1;
 36:   for (i=0; i<w; i++) {
 37:     fill[i] = nz;
 38:     for (j=0; j<w; j++) {
 39:       if (dfill[w*i+j]) {
 40:         fill[nz] = j;
 41:         nz++;
 42:       }
 43:     }
 44:   }
 45:   fill[w] = nz;

 47:   *rfill = fill;
 48:   return(0);
 49: }


 52: static PetscErrorCode DMDASetBlockFillsSparse_Private(const PetscInt *dfillsparse,PetscInt w,PetscInt **rfill)
 53: {
 55:   PetscInt       nz;

 58:   if (!dfillsparse) return(0);

 60:   /* Determine number of non-zeros */
 61:   nz = (dfillsparse[w] - w - 1);

 63:   /* Allocate space for our copy of the given sparse matrix representation. */
 64:   PetscMalloc1(nz + w + 1,rfill);
 65:   PetscMemcpy(*rfill,dfillsparse,(nz+w+1)*sizeof(PetscInt));
 66:   return(0);
 67: }


 70: static PetscErrorCode DMDASetBlockFills_Private2(DM_DA *dd)
 71: {
 73:   PetscInt       i,k,cnt = 1;


 77:   /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of
 78:    columns to the left with any nonzeros in them plus 1 */
 79:   PetscCalloc1(dd->w,&dd->ofillcols);
 80:   for (i=0; i<dd->w; i++) {
 81:     for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
 82:   }
 83:   for (i=0; i<dd->w; i++) {
 84:     if (dd->ofillcols[i]) {
 85:       dd->ofillcols[i] = cnt++;
 86:     }
 87:   }
 88:   return(0);
 89: }



 93: /*@
 94:     DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
 95:     of the matrix returned by DMCreateMatrix().

 97:     Logically Collective on DMDA

 99:     Input Parameter:
100: +   da - the distributed array
101: .   dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
102: -   ofill - the fill pattern in the off-diagonal blocks


105:     Level: developer

107:     Notes:
108:     This only makes sense when you are doing multicomponent problems but using the
109:        MPIAIJ matrix format

111:            The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
112:        representing coupling and 0 entries for missing coupling. For example
113: $             dfill[9] = {1, 0, 0,
114: $                         1, 1, 0,
115: $                         0, 1, 1}
116:        means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
117:        itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
118:        diagonal block).

120:      DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
121:      can be represented in the dfill, ofill format

123:    Contributed by Glenn Hammond

125: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()

127: @*/
128: PetscErrorCode  DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
129: {
130:   DM_DA          *dd = (DM_DA*)da->data;

134:   /* save the given dfill and ofill information */
135:   DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
136:   DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);

138:   /* count nonzeros in ofill columns */
139:   DMDASetBlockFills_Private2(dd);

141:   return(0);
142: }


145: /*@
146:     DMDASetBlockFillsSparse - Sets the fill pattern in each block for a multi-component problem
147:     of the matrix returned by DMCreateMatrix(), using sparse representations
148:     of fill patterns.

150:     Logically Collective on DMDA

152:     Input Parameter:
153: +   da - the distributed array
154: .   dfill - the sparse fill pattern in the diagonal block (may be NULL, means use dense block)
155: -   ofill - the sparse fill pattern in the off-diagonal blocks


158:     Level: developer

160:     Notes: This only makes sense when you are doing multicomponent problems but using the
161:        MPIAIJ matrix format

163:            The format for dfill and ofill is a sparse representation of a
164:            dof-by-dof matrix with 1 entries representing coupling and 0 entries
165:            for missing coupling.  The sparse representation is a 1 dimensional
166:            array of length nz + dof + 1, where nz is the number of non-zeros in
167:            the matrix.  The first dof entries in the array give the
168:            starting array indices of each row's items in the rest of the array,
169:            the dof+1st item contains the value nz + dof + 1 (i.e. the entire length of the array)
170:            and the remaining nz items give the column indices of each of
171:            the 1s within the logical 2D matrix.  Each row's items within
172:            the array are the column indices of the 1s within that row
173:            of the 2D matrix.  PETSc developers may recognize that this is the
174:            same format as that computed by the DMDASetBlockFills_Private()
175:            function from a dense 2D matrix representation.

177:      DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
178:      can be represented in the dfill, ofill format

180:    Contributed by Philip C. Roth

182: .seealso DMDASetBlockFills(), DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()

184: @*/
185: PetscErrorCode  DMDASetBlockFillsSparse(DM da,const PetscInt *dfillsparse,const PetscInt *ofillsparse)
186: {
187:   DM_DA          *dd = (DM_DA*)da->data;

191:   /* save the given dfill and ofill information */
192:   DMDASetBlockFillsSparse_Private(dfillsparse,dd->w,&dd->dfill);
193:   DMDASetBlockFillsSparse_Private(ofillsparse,dd->w,&dd->ofill);

195:   /* count nonzeros in ofill columns */
196:   DMDASetBlockFills_Private2(dd);

198:   return(0);
199: }


202: PetscErrorCode  DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
203: {
204:   PetscErrorCode   ierr;
205:   PetscInt         dim,m,n,p,nc;
206:   DMBoundaryType bx,by,bz;
207:   MPI_Comm         comm;
208:   PetscMPIInt      size;
209:   PetscBool        isBAIJ;
210:   DM_DA            *dd = (DM_DA*)da->data;

213:   /*
214:                                   m
215:           ------------------------------------------------------
216:          |                                                     |
217:          |                                                     |
218:          |               ----------------------                |
219:          |               |                    |                |
220:       n  |           yn  |                    |                |
221:          |               |                    |                |
222:          |               .---------------------                |
223:          |             (xs,ys)     xn                          |
224:          |            .                                        |
225:          |         (gxs,gys)                                   |
226:          |                                                     |
227:           -----------------------------------------------------
228:   */

230:   /*
231:          nc - number of components per grid point
232:          col - number of colors needed in one direction for single component problem

234:   */
235:   DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);

237:   PetscObjectGetComm((PetscObject)da,&comm);
238:   MPI_Comm_size(comm,&size);
239:   if (ctype == IS_COLORING_LOCAL) {
240:     if (size == 1) {
241:       ctype = IS_COLORING_GLOBAL;
242:     } else if (dim > 1) {
243:       if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
244:         SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_LOCAL cannot be used for periodic boundary condition having both ends of the domain  on the same process");
245:       }
246:     }
247:   }

249:   /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
250:      matrices is for the blocks, not the individual matrix elements  */
251:   PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);
252:   if (!isBAIJ) {PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);}
253:   if (!isBAIJ) {PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);}
254:   if (isBAIJ) {
255:     dd->w  = 1;
256:     dd->xs = dd->xs/nc;
257:     dd->xe = dd->xe/nc;
258:     dd->Xs = dd->Xs/nc;
259:     dd->Xe = dd->Xe/nc;
260:   }

262:   /*
263:      We do not provide a getcoloring function in the DMDA operations because
264:    the basic DMDA does not know about matrices. We think of DMDA as being more
265:    more low-level then matrices.
266:   */
267:   if (dim == 1) {
268:     DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
269:   } else if (dim == 2) {
270:      DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
271:   } else if (dim == 3) {
272:      DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
273:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
274:   if (isBAIJ) {
275:     dd->w  = nc;
276:     dd->xs = dd->xs*nc;
277:     dd->xe = dd->xe*nc;
278:     dd->Xs = dd->Xs*nc;
279:     dd->Xe = dd->Xe*nc;
280:   }
281:   return(0);
282: }

284: /* ---------------------------------------------------------------------------------*/

286: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
287: {
288:   PetscErrorCode   ierr;
289:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
290:   PetscInt         ncolors;
291:   MPI_Comm         comm;
292:   DMBoundaryType bx,by;
293:   DMDAStencilType  st;
294:   ISColoringValue  *colors;
295:   DM_DA            *dd = (DM_DA*)da->data;

298:   /*
299:          nc - number of components per grid point
300:          col - number of colors needed in one direction for single component problem

302:   */
303:   DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
304:   col  = 2*s + 1;
305:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
306:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
307:   PetscObjectGetComm((PetscObject)da,&comm);

309:   /* special case as taught to us by Paul Hovland */
310:   if (st == DMDA_STENCIL_STAR && s == 1) {
311:     DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
312:   } else {

314:     if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\
315:                                                             by 2*stencil_width + 1 (%d)\n", m, col);
316:     if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\
317:                                                             by 2*stencil_width + 1 (%d)\n", n, col);
318:     if (ctype == IS_COLORING_GLOBAL) {
319:       if (!dd->localcoloring) {
320:         PetscMalloc1(nc*nx*ny,&colors);
321:         ii   = 0;
322:         for (j=ys; j<ys+ny; j++) {
323:           for (i=xs; i<xs+nx; i++) {
324:             for (k=0; k<nc; k++) {
325:               colors[ii++] = k + nc*((i % col) + col*(j % col));
326:             }
327:           }
328:         }
329:         ncolors = nc + nc*(col-1 + col*(col-1));
330:         ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
331:       }
332:       *coloring = dd->localcoloring;
333:     } else if (ctype == IS_COLORING_LOCAL) {
334:       if (!dd->ghostedcoloring) {
335:         PetscMalloc1(nc*gnx*gny,&colors);
336:         ii   = 0;
337:         for (j=gys; j<gys+gny; j++) {
338:           for (i=gxs; i<gxs+gnx; i++) {
339:             for (k=0; k<nc; k++) {
340:               /* the complicated stuff is to handle periodic boundaries */
341:               colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
342:             }
343:           }
344:         }
345:         ncolors = nc + nc*(col - 1 + col*(col-1));
346:         ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
347:         /* PetscIntView(ncolors,(PetscInt*)colors,0); */

349:         ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
350:       }
351:       *coloring = dd->ghostedcoloring;
352:     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
353:   }
354:   ISColoringReference(*coloring);
355:   return(0);
356: }

358: /* ---------------------------------------------------------------------------------*/

360: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
361: {
362:   PetscErrorCode   ierr;
363:   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;
364:   PetscInt         ncolors;
365:   MPI_Comm         comm;
366:   DMBoundaryType bx,by,bz;
367:   DMDAStencilType  st;
368:   ISColoringValue  *colors;
369:   DM_DA            *dd = (DM_DA*)da->data;

372:   /*
373:          nc - number of components per grid point
374:          col - number of colors needed in one direction for single component problem

376:   */
377:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
378:   col  = 2*s + 1;
379:   if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
380:                                                          by 2*stencil_width + 1\n");
381:   if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
382:                                                          by 2*stencil_width + 1\n");
383:   if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
384:                                                          by 2*stencil_width + 1\n");

386:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
387:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
388:   PetscObjectGetComm((PetscObject)da,&comm);

390:   /* create the coloring */
391:   if (ctype == IS_COLORING_GLOBAL) {
392:     if (!dd->localcoloring) {
393:       PetscMalloc1(nc*nx*ny*nz,&colors);
394:       ii   = 0;
395:       for (k=zs; k<zs+nz; k++) {
396:         for (j=ys; j<ys+ny; j++) {
397:           for (i=xs; i<xs+nx; i++) {
398:             for (l=0; l<nc; l++) {
399:               colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
400:             }
401:           }
402:         }
403:       }
404:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
405:       ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
406:     }
407:     *coloring = dd->localcoloring;
408:   } else if (ctype == IS_COLORING_LOCAL) {
409:     if (!dd->ghostedcoloring) {
410:       PetscMalloc1(nc*gnx*gny*gnz,&colors);
411:       ii   = 0;
412:       for (k=gzs; k<gzs+gnz; k++) {
413:         for (j=gys; j<gys+gny; j++) {
414:           for (i=gxs; i<gxs+gnx; i++) {
415:             for (l=0; l<nc; l++) {
416:               /* the complicated stuff is to handle periodic boundaries */
417:               colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
418:             }
419:           }
420:         }
421:       }
422:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
423:       ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
424:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
425:     }
426:     *coloring = dd->ghostedcoloring;
427:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
428:   ISColoringReference(*coloring);
429:   return(0);
430: }

432: /* ---------------------------------------------------------------------------------*/

434: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
435: {
436:   PetscErrorCode   ierr;
437:   PetscInt         xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
438:   PetscInt         ncolors;
439:   MPI_Comm         comm;
440:   DMBoundaryType bx;
441:   ISColoringValue  *colors;
442:   DM_DA            *dd = (DM_DA*)da->data;

445:   /*
446:          nc - number of components per grid point
447:          col - number of colors needed in one direction for single component problem

449:   */
450:   DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
451:   col  = 2*s + 1;

453:   if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\
454:                                                           by 2*stencil_width + 1 %d\n",(int)m,(int)col);

456:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
457:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
458:   PetscObjectGetComm((PetscObject)da,&comm);

460:   /* create the coloring */
461:   if (ctype == IS_COLORING_GLOBAL) {
462:     if (!dd->localcoloring) {
463:       PetscMalloc1(nc*nx,&colors);
464:       if (dd->ofillcols) {
465:         PetscInt tc = 0;
466:         for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
467:         i1 = 0;
468:         for (i=xs; i<xs+nx; i++) {
469:           for (l=0; l<nc; l++) {
470:             if (dd->ofillcols[l] && (i % col)) {
471:               colors[i1++] =  nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
472:             } else {
473:               colors[i1++] = l;
474:             }
475:           }
476:         }
477:         ncolors = nc + 2*s*tc;
478:       } else {
479:         i1 = 0;
480:         for (i=xs; i<xs+nx; i++) {
481:           for (l=0; l<nc; l++) {
482:             colors[i1++] = l + nc*(i % col);
483:           }
484:         }
485:         ncolors = nc + nc*(col-1);
486:       }
487:       ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
488:     }
489:     *coloring = dd->localcoloring;
490:   } else if (ctype == IS_COLORING_LOCAL) {
491:     if (!dd->ghostedcoloring) {
492:       PetscMalloc1(nc*gnx,&colors);
493:       i1   = 0;
494:       for (i=gxs; i<gxs+gnx; i++) {
495:         for (l=0; l<nc; l++) {
496:           /* the complicated stuff is to handle periodic boundaries */
497:           colors[i1++] = l + nc*(SetInRange(i,m) % col);
498:         }
499:       }
500:       ncolors = nc + nc*(col-1);
501:       ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
502:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
503:     }
504:     *coloring = dd->ghostedcoloring;
505:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
506:   ISColoringReference(*coloring);
507:   return(0);
508: }

510: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
511: {
512:   PetscErrorCode   ierr;
513:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
514:   PetscInt         ncolors;
515:   MPI_Comm         comm;
516:   DMBoundaryType bx,by;
517:   ISColoringValue  *colors;
518:   DM_DA            *dd = (DM_DA*)da->data;

521:   /*
522:          nc - number of components per grid point
523:          col - number of colors needed in one direction for single component problem

525:   */
526:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
527:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
528:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
529:   PetscObjectGetComm((PetscObject)da,&comm);

531:   if (bx == DM_BOUNDARY_PERIODIC && (m % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible by 5\n");
532:   if (by == DM_BOUNDARY_PERIODIC && (n % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible by 5\n");

534:   /* create the coloring */
535:   if (ctype == IS_COLORING_GLOBAL) {
536:     if (!dd->localcoloring) {
537:       PetscMalloc1(nc*nx*ny,&colors);
538:       ii   = 0;
539:       for (j=ys; j<ys+ny; j++) {
540:         for (i=xs; i<xs+nx; i++) {
541:           for (k=0; k<nc; k++) {
542:             colors[ii++] = k + nc*((3*j+i) % 5);
543:           }
544:         }
545:       }
546:       ncolors = 5*nc;
547:       ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
548:     }
549:     *coloring = dd->localcoloring;
550:   } else if (ctype == IS_COLORING_LOCAL) {
551:     if (!dd->ghostedcoloring) {
552:       PetscMalloc1(nc*gnx*gny,&colors);
553:       ii = 0;
554:       for (j=gys; j<gys+gny; j++) {
555:         for (i=gxs; i<gxs+gnx; i++) {
556:           for (k=0; k<nc; k++) {
557:             colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
558:           }
559:         }
560:       }
561:       ncolors = 5*nc;
562:       ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
563:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
564:     }
565:     *coloring = dd->ghostedcoloring;
566:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
567:   return(0);
568: }

570: /* =========================================================================== */
571: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
572: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
573: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
574: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
575: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
576: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
577: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
578: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
579: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
580: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
581: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM,Mat);
582: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM,Mat);
583: extern PetscErrorCode DMCreateMatrix_DA_IS(DM,Mat);

585: /*@C
586:    MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix

588:    Logically Collective on Mat

590:    Input Parameters:
591: +  mat - the matrix
592: -  da - the da

594:    Level: intermediate

596: @*/
597: PetscErrorCode MatSetupDM(Mat mat,DM da)
598: {

604:   PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
605:   return(0);
606: }

608: PetscErrorCode  MatView_MPI_DA(Mat A,PetscViewer viewer)
609: {
610:   DM                da;
611:   PetscErrorCode    ierr;
612:   const char        *prefix;
613:   Mat               Anatural;
614:   AO                ao;
615:   PetscInt          rstart,rend,*petsc,i;
616:   IS                is;
617:   MPI_Comm          comm;
618:   PetscViewerFormat format;

621:   /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
622:   PetscViewerGetFormat(viewer,&format);
623:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);

625:   PetscObjectGetComm((PetscObject)A,&comm);
626:   MatGetDM(A, &da);
627:   if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");

629:   DMDAGetAO(da,&ao);
630:   MatGetOwnershipRange(A,&rstart,&rend);
631:   PetscMalloc1(rend-rstart,&petsc);
632:   for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
633:   AOApplicationToPetsc(ao,rend-rstart,petsc);
634:   ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);

636:   /* call viewer on natural ordering */
637:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
638:   ISDestroy(&is);
639:   PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
640:   PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
641:   PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
642:   ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
643:   MatView(Anatural,viewer);
644:   ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
645:   MatDestroy(&Anatural);
646:   return(0);
647: }

649: PetscErrorCode  MatLoad_MPI_DA(Mat A,PetscViewer viewer)
650: {
651:   DM             da;
653:   Mat            Anatural,Aapp;
654:   AO             ao;
655:   PetscInt       rstart,rend,*app,i,m,n,M,N;
656:   IS             is;
657:   MPI_Comm       comm;

660:   PetscObjectGetComm((PetscObject)A,&comm);
661:   MatGetDM(A, &da);
662:   if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");

664:   /* Load the matrix in natural ordering */
665:   MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
666:   MatSetType(Anatural,((PetscObject)A)->type_name);
667:   MatGetSize(A,&M,&N);
668:   MatGetLocalSize(A,&m,&n);
669:   MatSetSizes(Anatural,m,n,M,N);
670:   MatLoad(Anatural,viewer);

672:   /* Map natural ordering to application ordering and create IS */
673:   DMDAGetAO(da,&ao);
674:   MatGetOwnershipRange(Anatural,&rstart,&rend);
675:   PetscMalloc1(rend-rstart,&app);
676:   for (i=rstart; i<rend; i++) app[i-rstart] = i;
677:   AOPetscToApplication(ao,rend-rstart,app);
678:   ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);

680:   /* Do permutation and replace header */
681:   MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
682:   MatHeaderReplace(A,&Aapp);
683:   ISDestroy(&is);
684:   MatDestroy(&Anatural);
685:   return(0);
686: }

688: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
689: {
691:   PetscInt       dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
692:   Mat            A;
693:   MPI_Comm       comm;
694:   MatType        Atype;
695:   PetscSection   section, sectionGlobal;
696:   void           (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL,(*sell)(void)=NULL,(*is)(void)=NULL;
697:   MatType        mtype;
698:   PetscMPIInt    size;
699:   DM_DA          *dd = (DM_DA*)da->data;

702:   MatInitializePackage();
703:   mtype = da->mattype;

705:   DMGetSection(da, &section);
706:   if (section) {
707:     PetscInt  bs = -1;
708:     PetscInt  localSize;
709:     PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;

711:     DMGetGlobalSection(da, &sectionGlobal);
712:     PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
713:     MatCreate(PetscObjectComm((PetscObject)da),&A);
714:     MatSetSizes(A,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
715:     MatSetType(A,mtype);
716:     PetscStrcmp(mtype,MATSHELL,&isShell);
717:     PetscStrcmp(mtype,MATBAIJ,&isBlock);
718:     PetscStrcmp(mtype,MATSEQBAIJ,&isSeqBlock);
719:     PetscStrcmp(mtype,MATMPIBAIJ,&isMPIBlock);
720:     PetscStrcmp(mtype,MATSBAIJ,&isSymBlock);
721:     PetscStrcmp(mtype,MATSEQSBAIJ,&isSymSeqBlock);
722:     PetscStrcmp(mtype,MATMPISBAIJ,&isSymMPIBlock);
723:     /* Check for symmetric storage */
724:     isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
725:     if (isSymmetric) {
726:       MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
727:     }
728:     if (!isShell) {
729:       PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;

731:       if (bs < 0) {
732:         if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
733:           PetscInt pStart, pEnd, p, dof;

735:           PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
736:           for (p = pStart; p < pEnd; ++p) {
737:             PetscSectionGetDof(sectionGlobal, p, &dof);
738:             if (dof) {
739:               bs = dof;
740:               break;
741:             }
742:           }
743:         } else {
744:           bs = 1;
745:         }
746:         /* Must have same blocksize on all procs (some might have no points) */
747:         bsLocal = bs;
748:         MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
749:       }
750:       PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
751:       /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
752:       PetscFree4(dnz, onz, dnzu, onzu);
753:     }
754:   }
755:   /*
756:                                   m
757:           ------------------------------------------------------
758:          |                                                     |
759:          |                                                     |
760:          |               ----------------------                |
761:          |               |                    |                |
762:       n  |           ny  |                    |                |
763:          |               |                    |                |
764:          |               .---------------------                |
765:          |             (xs,ys)     nx                          |
766:          |            .                                        |
767:          |         (gxs,gys)                                   |
768:          |                                                     |
769:           -----------------------------------------------------
770:   */

772:   /*
773:          nc - number of components per grid point
774:          col - number of colors needed in one direction for single component problem

776:   */
777:   M   = dd->M;
778:   N   = dd->N;
779:   P   = dd->P;
780:   dim = da->dim;
781:   dof = dd->w;
782:   /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
783:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
784:   PetscObjectGetComm((PetscObject)da,&comm);
785:   MatCreate(comm,&A);
786:   MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
787:   MatSetType(A,mtype);
788:   MatSetFromOptions(A);
789:   MatSetDM(A,da);
790:   if (da->structure_only) {
791:     MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);
792:   }
793:   MatGetType(A,&Atype);
794:   /*
795:      We do not provide a getmatrix function in the DMDA operations because
796:    the basic DMDA does not know about matrices. We think of DMDA as being more
797:    more low-level than matrices. This is kind of cheating but, cause sometimes
798:    we think of DMDA has higher level than matrices.

800:      We could switch based on Atype (or mtype), but we do not since the
801:    specialized setting routines depend only on the particular preallocation
802:    details of the matrix, not the type itself.
803:   */
804:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
805:   if (!aij) {
806:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
807:   }
808:   if (!aij) {
809:     PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
810:     if (!baij) {
811:       PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
812:     }
813:     if (!baij) {
814:       PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
815:       if (!sbaij) {
816:         PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
817:       }
818:       if (!sbaij) {
819:         PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);
820:         if (!sell) {
821:           PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);
822:         }
823:       }
824:       if (!sell) {
825:         PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
826:       }
827:     }
828:   }
829:   if (aij) {
830:     if (dim == 1) {
831:       if (dd->ofill) {
832:         DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
833:       } else {
834:         DMCreateMatrix_DA_1d_MPIAIJ(da,A);
835:       }
836:     } else if (dim == 2) {
837:       if (dd->ofill) {
838:         DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
839:       } else {
840:         DMCreateMatrix_DA_2d_MPIAIJ(da,A);
841:       }
842:     } else if (dim == 3) {
843:       if (dd->ofill) {
844:         DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
845:       } else {
846:         DMCreateMatrix_DA_3d_MPIAIJ(da,A);
847:       }
848:     }
849:   } else if (baij) {
850:     if (dim == 2) {
851:       DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
852:     } else if (dim == 3) {
853:       DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
854:     } 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);
855:   } else if (sbaij) {
856:     if (dim == 2) {
857:       DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
858:     } else if (dim == 3) {
859:       DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
860:     } 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);
861:   } else if (sell) {
862:      if (dim == 2) {
863:        DMCreateMatrix_DA_2d_MPISELL(da,A);
864:      } else if (dim == 3) {
865:        DMCreateMatrix_DA_3d_MPISELL(da,A);
866:      } 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);
867:   } else if (is) {
868:     DMCreateMatrix_DA_IS(da,A);
869:   } else {
870:     ISLocalToGlobalMapping ltog;

872:     MatSetBlockSize(A,dof);
873:     MatSetUp(A);
874:     DMGetLocalToGlobalMapping(da,&ltog);
875:     MatSetLocalToGlobalMapping(A,ltog,ltog);
876:   }
877:   DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
878:   MatSetStencil(A,dim,dims,starts,dof);
879:   MatSetDM(A,da);
880:   MPI_Comm_size(comm,&size);
881:   if (size > 1) {
882:     /* change viewer to display matrix in natural ordering */
883:     MatSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
884:     MatSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
885:   }
886:   *J = A;
887:   return(0);
888: }

890: /* ---------------------------------------------------------------------------------*/
891: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]);

893: PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J)
894: {
895:   DM_DA                  *da = (DM_DA*)dm->data;
896:   Mat                    lJ;
897:   ISLocalToGlobalMapping ltog;
898:   IS                     is_loc_filt, is_glob;
899:   const PetscInt         *e_loc,*idx;
900:   PetscInt               nel,nen,nv,dof,dim,*gidx,nb;
901:   PetscBool              flg;
902:   PetscErrorCode         ierr;

904:   /* The l2g map of DMDA has all ghosted nodes, and e_loc is a subset of all the local nodes (including the ghosted)
905:      We need to filter the local indices that are represented through the DMDAGetElements decomposition
906:      This is because the size of the local matrices in MATIS is the local size of the l2g map */
908:   dof  = da->w;
909:   dim  = dm->dim;

911:   MatSetBlockSize(J,dof);

913:   /* get local elements indices in local DMDA numbering */
914:   DMDAGetElements(dm,&nel,&nen,&e_loc); /* this will throw an error if the stencil type is not DMDA_STENCIL_BOX */
915:   ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nel*nen,e_loc,PETSC_COPY_VALUES,&is_loc_filt);
916:   DMDARestoreElements(dm,&nel,&nen,&e_loc);

918:   /* obtain a consistent local ordering for MATIS */
919:   ISSortRemoveDups(is_loc_filt);
920:   ISBlockGetLocalSize(is_loc_filt,&nb);
921:   DMGetLocalToGlobalMapping(dm,&ltog);
922:   ISLocalToGlobalMappingGetSize(ltog,&nv);
923:   PetscMalloc1(PetscMax(nb,nv/dof),&gidx);
924:   ISBlockGetIndices(is_loc_filt,&idx);
925:   ISLocalToGlobalMappingApplyBlock(ltog,nb,idx,gidx);
926:   ISBlockRestoreIndices(is_loc_filt,&idx);
927:   ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nb,gidx,PETSC_USE_POINTER,&is_glob);
928:   ISLocalToGlobalMappingCreateIS(is_glob,&ltog);
929:   ISDestroy(&is_glob);
930:   MatSetLocalToGlobalMapping(J,ltog,ltog);
931:   ISLocalToGlobalMappingDestroy(&ltog);

933:   /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */
934:   MatISGetLocalMat(J,&lJ);
935:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
936:   ISDestroy(&is_loc_filt);
937:   ISCreateStride(PetscObjectComm((PetscObject)lJ),nv/dof,0,1,&is_glob);
938:   ISGetIndices(is_glob,&idx);
939:   ISGlobalToLocalMappingApplyBlock(ltog,IS_GTOLM_MASK,nv/dof,idx,&nb,gidx);
940:   ISRestoreIndices(is_glob,&idx);
941:   ISDestroy(&is_glob);
942:   ISLocalToGlobalMappingDestroy(&ltog);
943:   ISCreateBlock(PETSC_COMM_SELF,dof,nb,gidx,PETSC_USE_POINTER,&is_loc_filt);
944:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
945:   ISDestroy(&is_loc_filt);
946:   MatSetLocalToGlobalMapping(lJ,ltog,ltog);
947:   ISLocalToGlobalMappingDestroy(&ltog);
948:   PetscFree(gidx);

950:   /* Preallocation (not exact): we reuse the preallocation routines of the assembled version  */
951:   flg = dm->prealloc_only;
952:   dm->prealloc_only = PETSC_TRUE;
953:   switch (dim) {
954:   case 1:
955:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
956:     DMCreateMatrix_DA_1d_MPIAIJ(dm,J);
957:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
958:     break;
959:   case 2:
960:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
961:     DMCreateMatrix_DA_2d_MPIAIJ(dm,J);
962:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
963:     break;
964:   case 3:
965:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
966:     DMCreateMatrix_DA_3d_MPIAIJ(dm,J);
967:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
968:     break;
969:   default:
970:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dimension %d",dim);
971:     break;
972:   }
973:   dm->prealloc_only = flg;
974:   return(0);
975: }

977: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J)
978: {
979:   PetscErrorCode         ierr;
980:   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;
981:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
982:   MPI_Comm               comm;
983:   PetscScalar            *values;
984:   DMBoundaryType         bx,by;
985:   ISLocalToGlobalMapping ltog;
986:   DMDAStencilType        st;

989:   /*
990:          nc - number of components per grid point
991:          col - number of colors needed in one direction for single component problem

993:   */
994:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
995:   col  = 2*s + 1;
996:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
997:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
998:   PetscObjectGetComm((PetscObject)da,&comm);

1000:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1001:   DMGetLocalToGlobalMapping(da,&ltog);

1003:   MatSetBlockSize(J,nc);
1004:   /* determine the matrix preallocation information */
1005:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1006:   for (i=xs; i<xs+nx; i++) {

1008:     pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1009:     pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1011:     for (j=ys; j<ys+ny; j++) {
1012:       slot = i - gxs + gnx*(j - gys);

1014:       lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1015:       lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1017:       cnt = 0;
1018:       for (k=0; k<nc; k++) {
1019:         for (l=lstart; l<lend+1; l++) {
1020:           for (p=pstart; p<pend+1; p++) {
1021:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1022:               cols[cnt++] = k + nc*(slot + gnx*l + p);
1023:             }
1024:           }
1025:         }
1026:         rows[k] = k + nc*(slot);
1027:       }
1028:       MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1029:     }
1030:   }
1031:   MatSetBlockSize(J,nc);
1032:   MatSeqSELLSetPreallocation(J,0,dnz);
1033:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1034:   MatPreallocateFinalize(dnz,onz);

1036:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1038:   /*
1039:     For each node in the grid: we get the neighbors in the local (on processor ordering
1040:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1041:     PETSc ordering.
1042:   */
1043:   if (!da->prealloc_only) {
1044:     PetscCalloc1(col*col*nc*nc,&values);
1045:     for (i=xs; i<xs+nx; i++) {

1047:       pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1048:       pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1050:       for (j=ys; j<ys+ny; j++) {
1051:         slot = i - gxs + gnx*(j - gys);

1053:         lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1054:         lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1056:         cnt = 0;
1057:         for (k=0; k<nc; k++) {
1058:           for (l=lstart; l<lend+1; l++) {
1059:             for (p=pstart; p<pend+1; p++) {
1060:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1061:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
1062:               }
1063:             }
1064:           }
1065:           rows[k] = k + nc*(slot);
1066:         }
1067:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1068:       }
1069:     }
1070:     PetscFree(values);
1071:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1072:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1073:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1074:   }
1075:   PetscFree2(rows,cols);
1076:   return(0);
1077: }

1079: PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da,Mat J)
1080: {
1081:   PetscErrorCode         ierr;
1082:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1083:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1084:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1085:   MPI_Comm               comm;
1086:   PetscScalar            *values;
1087:   DMBoundaryType         bx,by,bz;
1088:   ISLocalToGlobalMapping ltog;
1089:   DMDAStencilType        st;

1092:   /*
1093:          nc - number of components per grid point
1094:          col - number of colors needed in one direction for single component problem

1096:   */
1097:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1098:   col  = 2*s + 1;
1099:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1100:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1101:   PetscObjectGetComm((PetscObject)da,&comm);

1103:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1104:   DMGetLocalToGlobalMapping(da,&ltog);

1106:   MatSetBlockSize(J,nc);
1107:   /* determine the matrix preallocation information */
1108:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1109:   for (i=xs; i<xs+nx; i++) {
1110:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1111:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1112:     for (j=ys; j<ys+ny; j++) {
1113:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1114:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1115:       for (k=zs; k<zs+nz; k++) {
1116:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1117:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1119:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

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:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1135:       }
1136:     }
1137:   }
1138:   MatSetBlockSize(J,nc);
1139:   MatSeqSELLSetPreallocation(J,0,dnz);
1140:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1141:   MatPreallocateFinalize(dnz,onz);
1142:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1144:   /*
1145:     For each node in the grid: we get the neighbors in the local (on processor ordering
1146:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1147:     PETSc ordering.
1148:   */
1149:   if (!da->prealloc_only) {
1150:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1151:     for (i=xs; i<xs+nx; i++) {
1152:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1153:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1154:       for (j=ys; j<ys+ny; j++) {
1155:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1156:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1157:         for (k=zs; k<zs+nz; k++) {
1158:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1159:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1161:           slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

1163:           cnt = 0;
1164:           for (l=0; l<nc; l++) {
1165:             for (ii=istart; ii<iend+1; ii++) {
1166:               for (jj=jstart; jj<jend+1; jj++) {
1167:                 for (kk=kstart; kk<kend+1; kk++) {
1168:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1169:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1170:                   }
1171:                 }
1172:               }
1173:             }
1174:             rows[l] = l + nc*(slot);
1175:           }
1176:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1177:         }
1178:       }
1179:     }
1180:     PetscFree(values);
1181:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1182:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1183:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1184:   }
1185:   PetscFree2(rows,cols);
1186:   return(0);
1187: }

1189: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
1190: {
1191:   PetscErrorCode         ierr;
1192:   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,M,N;
1193:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1194:   MPI_Comm               comm;
1195:   PetscScalar            *values;
1196:   DMBoundaryType         bx,by;
1197:   ISLocalToGlobalMapping ltog,mltog;
1198:   DMDAStencilType        st;
1199:   PetscBool              removedups = PETSC_FALSE;

1202:   /*
1203:          nc - number of components per grid point
1204:          col - number of colors needed in one direction for single component problem

1206:   */
1207:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1208:   col  = 2*s + 1;
1209:   /*
1210:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1211:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1212:   */
1213:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1214:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1215:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1216:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1217:   PetscObjectGetComm((PetscObject)da,&comm);

1219:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1220:   DMGetLocalToGlobalMapping(da,&ltog);

1222:   MatSetBlockSize(J,nc);
1223:   /* determine the matrix preallocation information */
1224:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1225:   for (i=xs; i<xs+nx; i++) {

1227:     pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1228:     pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1230:     for (j=ys; j<ys+ny; j++) {
1231:       slot = i - gxs + gnx*(j - gys);

1233:       lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1234:       lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1236:       cnt = 0;
1237:       for (k=0; k<nc; k++) {
1238:         for (l=lstart; l<lend+1; l++) {
1239:           for (p=pstart; p<pend+1; p++) {
1240:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1241:               cols[cnt++] = k + nc*(slot + gnx*l + p);
1242:             }
1243:           }
1244:         }
1245:         rows[k] = k + nc*(slot);
1246:       }
1247:       if (removedups) {
1248:         MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1249:       } else {
1250:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1251:       }
1252:     }
1253:   }
1254:   MatSetBlockSize(J,nc);
1255:   MatSeqAIJSetPreallocation(J,0,dnz);
1256:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1257:   MatPreallocateFinalize(dnz,onz);
1258:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1259:   if (!mltog) {
1260:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1261:   }

1263:   /*
1264:     For each node in the grid: we get the neighbors in the local (on processor ordering
1265:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1266:     PETSc ordering.
1267:   */
1268:   if (!da->prealloc_only) {
1269:     PetscCalloc1(col*col*nc*nc,&values);
1270:     for (i=xs; i<xs+nx; i++) {

1272:       pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1273:       pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1275:       for (j=ys; j<ys+ny; j++) {
1276:         slot = i - gxs + gnx*(j - gys);

1278:         lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1279:         lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1281:         cnt = 0;
1282:         for (k=0; k<nc; k++) {
1283:           for (l=lstart; l<lend+1; l++) {
1284:             for (p=pstart; p<pend+1; p++) {
1285:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1286:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
1287:               }
1288:             }
1289:           }
1290:           rows[k] = k + nc*(slot);
1291:         }
1292:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1293:       }
1294:     }
1295:     PetscFree(values);
1296:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1297:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1298:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1299:   }
1300:   PetscFree2(rows,cols);
1301:   return(0);
1302: }

1304: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
1305: {
1306:   PetscErrorCode         ierr;
1307:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1308:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
1309:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1310:   DM_DA                  *dd = (DM_DA*)da->data;
1311:   PetscInt               ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
1312:   MPI_Comm               comm;
1313:   PetscScalar            *values;
1314:   DMBoundaryType         bx,by;
1315:   ISLocalToGlobalMapping ltog;
1316:   DMDAStencilType        st;
1317:   PetscBool              removedups = PETSC_FALSE;

1320:   /*
1321:          nc - number of components per grid point
1322:          col - number of colors needed in one direction for single component problem

1324:   */
1325:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1326:   col  = 2*s + 1;
1327:   /*
1328:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1329:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1330:   */
1331:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1332:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1333:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1334:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1335:   PetscObjectGetComm((PetscObject)da,&comm);

1337:   PetscMalloc1(col*col*nc,&cols);
1338:   DMGetLocalToGlobalMapping(da,&ltog);

1340:   MatSetBlockSize(J,nc);
1341:   /* determine the matrix preallocation information */
1342:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1343:   for (i=xs; i<xs+nx; i++) {

1345:     pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1346:     pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1348:     for (j=ys; j<ys+ny; j++) {
1349:       slot = i - gxs + gnx*(j - gys);

1351:       lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1352:       lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1354:       for (k=0; k<nc; k++) {
1355:         cnt = 0;
1356:         for (l=lstart; l<lend+1; l++) {
1357:           for (p=pstart; p<pend+1; p++) {
1358:             if (l || p) {
1359:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1360:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1361:               }
1362:             } else {
1363:               if (dfill) {
1364:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1365:               } else {
1366:                 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1367:               }
1368:             }
1369:           }
1370:         }
1371:         row    = k + nc*(slot);
1372:         maxcnt = PetscMax(maxcnt,cnt);
1373:         if (removedups) {
1374:           MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1375:         } else {
1376:           MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1377:         }
1378:       }
1379:     }
1380:   }
1381:   MatSeqAIJSetPreallocation(J,0,dnz);
1382:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1383:   MatPreallocateFinalize(dnz,onz);
1384:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1386:   /*
1387:     For each node in the grid: we get the neighbors in the local (on processor ordering
1388:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1389:     PETSc ordering.
1390:   */
1391:   if (!da->prealloc_only) {
1392:     PetscCalloc1(maxcnt,&values);
1393:     for (i=xs; i<xs+nx; i++) {

1395:       pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1396:       pend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));

1398:       for (j=ys; j<ys+ny; j++) {
1399:         slot = i - gxs + gnx*(j - gys);

1401:         lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1402:         lend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));

1404:         for (k=0; k<nc; k++) {
1405:           cnt = 0;
1406:           for (l=lstart; l<lend+1; l++) {
1407:             for (p=pstart; p<pend+1; p++) {
1408:               if (l || p) {
1409:                 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1410:                   for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1411:                 }
1412:               } else {
1413:                 if (dfill) {
1414:                   for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1415:                 } else {
1416:                   for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1417:                 }
1418:               }
1419:             }
1420:           }
1421:           row  = k + nc*(slot);
1422:           MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1423:         }
1424:       }
1425:     }
1426:     PetscFree(values);
1427:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1428:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1429:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1430:   }
1431:   PetscFree(cols);
1432:   return(0);
1433: }

1435: /* ---------------------------------------------------------------------------------*/

1437: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1438: {
1439:   PetscErrorCode         ierr;
1440:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1441:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1442:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1443:   MPI_Comm               comm;
1444:   PetscScalar            *values;
1445:   DMBoundaryType         bx,by,bz;
1446:   ISLocalToGlobalMapping ltog,mltog;
1447:   DMDAStencilType        st;
1448:   PetscBool              removedups = PETSC_FALSE;

1451:   /*
1452:          nc - number of components per grid point
1453:          col - number of colors needed in one direction for single component problem

1455:   */
1456:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1457:   col  = 2*s + 1;

1459:   /*
1460:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1461:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1462:   */
1463:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1464:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1465:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

1467:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1468:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1469:   PetscObjectGetComm((PetscObject)da,&comm);

1471:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1472:   DMGetLocalToGlobalMapping(da,&ltog);

1474:   MatSetBlockSize(J,nc);
1475:   /* determine the matrix preallocation information */
1476:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1477:   for (i=xs; i<xs+nx; i++) {
1478:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1479:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1480:     for (j=ys; j<ys+ny; j++) {
1481:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1482:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1483:       for (k=zs; k<zs+nz; k++) {
1484:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1485:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1487:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

1489:         cnt = 0;
1490:         for (l=0; l<nc; l++) {
1491:           for (ii=istart; ii<iend+1; ii++) {
1492:             for (jj=jstart; jj<jend+1; jj++) {
1493:               for (kk=kstart; kk<kend+1; kk++) {
1494:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1495:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1496:                 }
1497:               }
1498:             }
1499:           }
1500:           rows[l] = l + nc*(slot);
1501:         }
1502:         if (removedups) {
1503:           MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1504:         } else {
1505:           MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1506:         }
1507:       }
1508:     }
1509:   }
1510:   MatSetBlockSize(J,nc);
1511:   MatSeqAIJSetPreallocation(J,0,dnz);
1512:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1513:   MatPreallocateFinalize(dnz,onz);
1514:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1515:   if (!mltog) {
1516:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1517:   }

1519:   /*
1520:     For each node in the grid: we get the neighbors in the local (on processor ordering
1521:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1522:     PETSc ordering.
1523:   */
1524:   if (!da->prealloc_only) {
1525:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1526:     for (i=xs; i<xs+nx; i++) {
1527:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1528:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1529:       for (j=ys; j<ys+ny; j++) {
1530:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1531:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1532:         for (k=zs; k<zs+nz; k++) {
1533:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1534:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1536:           slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

1538:           cnt = 0;
1539:           for (l=0; l<nc; l++) {
1540:             for (ii=istart; ii<iend+1; ii++) {
1541:               for (jj=jstart; jj<jend+1; jj++) {
1542:                 for (kk=kstart; kk<kend+1; kk++) {
1543:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1544:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1545:                   }
1546:                 }
1547:               }
1548:             }
1549:             rows[l] = l + nc*(slot);
1550:           }
1551:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1552:         }
1553:       }
1554:     }
1555:     PetscFree(values);
1556:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1557:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1558:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1559:   }
1560:   PetscFree2(rows,cols);
1561:   return(0);
1562: }

1564: /* ---------------------------------------------------------------------------------*/

1566: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1567: {
1568:   PetscErrorCode         ierr;
1569:   DM_DA                  *dd = (DM_DA*)da->data;
1570:   PetscInt               xs,nx,i,j,gxs,gnx,row,k,l;
1571:   PetscInt               m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1572:   PetscInt               *ofill = dd->ofill,*dfill = dd->dfill;
1573:   PetscScalar            *values;
1574:   DMBoundaryType         bx;
1575:   ISLocalToGlobalMapping ltog;
1576:   PetscMPIInt            rank,size;

1579:   if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1580:   MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1581:   MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);

1583:   /*
1584:          nc - number of components per grid point

1586:   */
1587:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1588:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1589:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1591:   MatSetBlockSize(J,nc);
1592:   PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);

1594:   /*
1595:         note should be smaller for first and last process with no periodic
1596:         does not handle dfill
1597:   */
1598:   cnt = 0;
1599:   /* coupling with process to the left */
1600:   for (i=0; i<s; i++) {
1601:     for (j=0; j<nc; j++) {
1602:       ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1603:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1604:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1605:       cnt++;
1606:     }
1607:   }
1608:   for (i=s; i<nx-s; i++) {
1609:     for (j=0; j<nc; j++) {
1610:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1611:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1612:       cnt++;
1613:     }
1614:   }
1615:   /* coupling with process to the right */
1616:   for (i=nx-s; i<nx; i++) {
1617:     for (j=0; j<nc; j++) {
1618:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1619:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1620:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1621:       cnt++;
1622:     }
1623:   }

1625:   MatSeqAIJSetPreallocation(J,0,cols);
1626:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1627:   PetscFree2(cols,ocols);

1629:   DMGetLocalToGlobalMapping(da,&ltog);
1630:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1632:   /*
1633:     For each node in the grid: we get the neighbors in the local (on processor ordering
1634:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1635:     PETSc ordering.
1636:   */
1637:   if (!da->prealloc_only) {
1638:     PetscCalloc2(maxcnt,&values,maxcnt,&cols);

1640:     row = xs*nc;
1641:     /* coupling with process to the left */
1642:     for (i=xs; i<xs+s; i++) {
1643:       for (j=0; j<nc; j++) {
1644:         cnt = 0;
1645:         if (rank) {
1646:           for (l=0; l<s; l++) {
1647:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1648:           }
1649:         }
1650:         if (dfill) {
1651:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1652:             cols[cnt++] = i*nc + dfill[k];
1653:           }
1654:         } else {
1655:           for (k=0; k<nc; k++) {
1656:             cols[cnt++] = i*nc + k;
1657:           }
1658:         }
1659:         for (l=0; l<s; l++) {
1660:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1661:         }
1662:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1663:         row++;
1664:       }
1665:     }
1666:     for (i=xs+s; i<xs+nx-s; i++) {
1667:       for (j=0; j<nc; j++) {
1668:         cnt = 0;
1669:         for (l=0; l<s; l++) {
1670:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1671:         }
1672:         if (dfill) {
1673:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1674:             cols[cnt++] = i*nc + dfill[k];
1675:           }
1676:         } else {
1677:           for (k=0; k<nc; k++) {
1678:             cols[cnt++] = i*nc + k;
1679:           }
1680:         }
1681:         for (l=0; l<s; l++) {
1682:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1683:         }
1684:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1685:         row++;
1686:       }
1687:     }
1688:     /* coupling with process to the right */
1689:     for (i=xs+nx-s; i<xs+nx; i++) {
1690:       for (j=0; j<nc; j++) {
1691:         cnt = 0;
1692:         for (l=0; l<s; l++) {
1693:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1694:         }
1695:         if (dfill) {
1696:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1697:             cols[cnt++] = i*nc + dfill[k];
1698:           }
1699:         } else {
1700:           for (k=0; k<nc; k++) {
1701:             cols[cnt++] = i*nc + k;
1702:           }
1703:         }
1704:         if (rank < size-1) {
1705:           for (l=0; l<s; l++) {
1706:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1707:           }
1708:         }
1709:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1710:         row++;
1711:       }
1712:     }
1713:     PetscFree2(values,cols);
1714:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1715:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1716:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1717:   }
1718:   return(0);
1719: }

1721: /* ---------------------------------------------------------------------------------*/

1723: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1724: {
1725:   PetscErrorCode         ierr;
1726:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1727:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1728:   PetscInt               istart,iend;
1729:   PetscScalar            *values;
1730:   DMBoundaryType         bx;
1731:   ISLocalToGlobalMapping ltog,mltog;

1734:   /*
1735:          nc - number of components per grid point
1736:          col - number of colors needed in one direction for single component problem

1738:   */
1739:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1740:   col  = 2*s + 1;

1742:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1743:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1745:   MatSetBlockSize(J,nc);
1746:   MatSeqAIJSetPreallocation(J,col*nc,0);
1747:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1749:   DMGetLocalToGlobalMapping(da,&ltog);
1750:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1751:   if (!mltog) {
1752:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1753:   }

1755:   /*
1756:     For each node in the grid: we get the neighbors in the local (on processor ordering
1757:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1758:     PETSc ordering.
1759:   */
1760:   if (!da->prealloc_only) {
1761:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1762:     PetscCalloc1(col*nc*nc,&values);
1763:     for (i=xs; i<xs+nx; i++) {
1764:       istart = PetscMax(-s,gxs - i);
1765:       iend   = PetscMin(s,gxs + gnx - i - 1);
1766:       slot   = i - gxs;

1768:       cnt = 0;
1769:       for (l=0; l<nc; l++) {
1770:         for (i1=istart; i1<iend+1; i1++) {
1771:           cols[cnt++] = l + nc*(slot + i1);
1772:         }
1773:         rows[l] = l + nc*(slot);
1774:       }
1775:       MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1776:     }
1777:     PetscFree(values);
1778:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1779:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1780:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1781:     PetscFree2(rows,cols);
1782:   }
1783:   return(0);
1784: }

1786: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1787: {
1788:   PetscErrorCode         ierr;
1789:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1790:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1791:   PetscInt               istart,iend,jstart,jend,ii,jj;
1792:   MPI_Comm               comm;
1793:   PetscScalar            *values;
1794:   DMBoundaryType         bx,by;
1795:   DMDAStencilType        st;
1796:   ISLocalToGlobalMapping ltog;

1799:   /*
1800:      nc - number of components per grid point
1801:      col - number of colors needed in one direction for single component problem
1802:   */
1803:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1804:   col  = 2*s + 1;

1806:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1807:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1808:   PetscObjectGetComm((PetscObject)da,&comm);

1810:   PetscMalloc1(col*col*nc*nc,&cols);

1812:   DMGetLocalToGlobalMapping(da,&ltog);

1814:   /* determine the matrix preallocation information */
1815:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1816:   for (i=xs; i<xs+nx; i++) {
1817:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1818:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1819:     for (j=ys; j<ys+ny; j++) {
1820:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1821:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1822:       slot   = i - gxs + gnx*(j - gys);

1824:       /* Find block columns in block row */
1825:       cnt = 0;
1826:       for (ii=istart; ii<iend+1; ii++) {
1827:         for (jj=jstart; jj<jend+1; jj++) {
1828:           if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1829:             cols[cnt++] = slot + ii + gnx*jj;
1830:           }
1831:         }
1832:       }
1833:       MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1834:     }
1835:   }
1836:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1837:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1838:   MatPreallocateFinalize(dnz,onz);

1840:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1842:   /*
1843:     For each node in the grid: we get the neighbors in the local (on processor ordering
1844:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1845:     PETSc ordering.
1846:   */
1847:   if (!da->prealloc_only) {
1848:     PetscCalloc1(col*col*nc*nc,&values);
1849:     for (i=xs; i<xs+nx; i++) {
1850:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1851:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1852:       for (j=ys; j<ys+ny; j++) {
1853:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1854:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1855:         slot = i - gxs + gnx*(j - gys);
1856:         cnt  = 0;
1857:         for (ii=istart; ii<iend+1; ii++) {
1858:           for (jj=jstart; jj<jend+1; jj++) {
1859:             if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1860:               cols[cnt++] = slot + ii + gnx*jj;
1861:             }
1862:           }
1863:         }
1864:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1865:       }
1866:     }
1867:     PetscFree(values);
1868:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1869:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1870:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1871:   }
1872:   PetscFree(cols);
1873:   return(0);
1874: }

1876: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1877: {
1878:   PetscErrorCode         ierr;
1879:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1880:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1881:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1882:   MPI_Comm               comm;
1883:   PetscScalar            *values;
1884:   DMBoundaryType         bx,by,bz;
1885:   DMDAStencilType        st;
1886:   ISLocalToGlobalMapping ltog;

1889:   /*
1890:          nc - number of components per grid point
1891:          col - number of colors needed in one direction for single component problem

1893:   */
1894:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1895:   col  = 2*s + 1;

1897:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1898:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1899:   PetscObjectGetComm((PetscObject)da,&comm);

1901:   PetscMalloc1(col*col*col,&cols);

1903:   DMGetLocalToGlobalMapping(da,&ltog);

1905:   /* determine the matrix preallocation information */
1906:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1907:   for (i=xs; i<xs+nx; i++) {
1908:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1909:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1910:     for (j=ys; j<ys+ny; j++) {
1911:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1912:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1913:       for (k=zs; k<zs+nz; k++) {
1914:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1915:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1917:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

1919:         /* Find block columns in block row */
1920:         cnt = 0;
1921:         for (ii=istart; ii<iend+1; ii++) {
1922:           for (jj=jstart; jj<jend+1; jj++) {
1923:             for (kk=kstart; kk<kend+1; kk++) {
1924:               if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1925:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1926:               }
1927:             }
1928:           }
1929:         }
1930:         MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1931:       }
1932:     }
1933:   }
1934:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1935:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1936:   MatPreallocateFinalize(dnz,onz);

1938:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1940:   /*
1941:     For each node in the grid: we get the neighbors in the local (on processor ordering
1942:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1943:     PETSc ordering.
1944:   */
1945:   if (!da->prealloc_only) {
1946:     PetscCalloc1(col*col*col*nc*nc,&values);
1947:     for (i=xs; i<xs+nx; i++) {
1948:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1949:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1950:       for (j=ys; j<ys+ny; j++) {
1951:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1952:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1953:         for (k=zs; k<zs+nz; k++) {
1954:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1955:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

1957:           slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

1959:           cnt = 0;
1960:           for (ii=istart; ii<iend+1; ii++) {
1961:             for (jj=jstart; jj<jend+1; jj++) {
1962:               for (kk=kstart; kk<kend+1; kk++) {
1963:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1964:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1965:                 }
1966:               }
1967:             }
1968:           }
1969:           MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1970:         }
1971:       }
1972:     }
1973:     PetscFree(values);
1974:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1975:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1976:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1977:   }
1978:   PetscFree(cols);
1979:   return(0);
1980: }

1982: /*
1983:   This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1984:   identify in the local ordering with periodic domain.
1985: */
1986: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1987: {
1989:   PetscInt       i,n;

1992:   ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1993:   ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1994:   for (i=0,n=0; i<*cnt; i++) {
1995:     if (col[i] >= *row) col[n++] = col[i];
1996:   }
1997:   *cnt = n;
1998:   return(0);
1999: }

2001: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
2002: {
2003:   PetscErrorCode         ierr;
2004:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2005:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
2006:   PetscInt               istart,iend,jstart,jend,ii,jj;
2007:   MPI_Comm               comm;
2008:   PetscScalar            *values;
2009:   DMBoundaryType         bx,by;
2010:   DMDAStencilType        st;
2011:   ISLocalToGlobalMapping ltog;

2014:   /*
2015:      nc - number of components per grid point
2016:      col - number of colors needed in one direction for single component problem
2017:   */
2018:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
2019:   col  = 2*s + 1;

2021:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
2022:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
2023:   PetscObjectGetComm((PetscObject)da,&comm);

2025:   PetscMalloc1(col*col*nc*nc,&cols);

2027:   DMGetLocalToGlobalMapping(da,&ltog);

2029:   /* determine the matrix preallocation information */
2030:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
2031:   for (i=xs; i<xs+nx; i++) {
2032:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2033:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2034:     for (j=ys; j<ys+ny; j++) {
2035:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2036:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2037:       slot   = i - gxs + gnx*(j - gys);

2039:       /* Find block columns in block row */
2040:       cnt = 0;
2041:       for (ii=istart; ii<iend+1; ii++) {
2042:         for (jj=jstart; jj<jend+1; jj++) {
2043:           if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2044:             cols[cnt++] = slot + ii + gnx*jj;
2045:           }
2046:         }
2047:       }
2048:       L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2049:       MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2050:     }
2051:   }
2052:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2053:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2054:   MatPreallocateFinalize(dnz,onz);

2056:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2058:   /*
2059:     For each node in the grid: we get the neighbors in the local (on processor ordering
2060:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2061:     PETSc ordering.
2062:   */
2063:   if (!da->prealloc_only) {
2064:     PetscCalloc1(col*col*nc*nc,&values);
2065:     for (i=xs; i<xs+nx; i++) {
2066:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2067:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2068:       for (j=ys; j<ys+ny; j++) {
2069:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2070:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2071:         slot   = i - gxs + gnx*(j - gys);

2073:         /* Find block columns in block row */
2074:         cnt = 0;
2075:         for (ii=istart; ii<iend+1; ii++) {
2076:           for (jj=jstart; jj<jend+1; jj++) {
2077:             if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2078:               cols[cnt++] = slot + ii + gnx*jj;
2079:             }
2080:           }
2081:         }
2082:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2083:         MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2084:       }
2085:     }
2086:     PetscFree(values);
2087:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2088:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2089:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2090:   }
2091:   PetscFree(cols);
2092:   return(0);
2093: }

2095: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
2096: {
2097:   PetscErrorCode         ierr;
2098:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2099:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
2100:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
2101:   MPI_Comm               comm;
2102:   PetscScalar            *values;
2103:   DMBoundaryType         bx,by,bz;
2104:   DMDAStencilType        st;
2105:   ISLocalToGlobalMapping ltog;

2108:   /*
2109:      nc - number of components per grid point
2110:      col - number of colors needed in one direction for single component problem
2111:   */
2112:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2113:   col  = 2*s + 1;

2115:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2116:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2117:   PetscObjectGetComm((PetscObject)da,&comm);

2119:   /* create the matrix */
2120:   PetscMalloc1(col*col*col,&cols);

2122:   DMGetLocalToGlobalMapping(da,&ltog);

2124:   /* determine the matrix preallocation information */
2125:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2126:   for (i=xs; i<xs+nx; i++) {
2127:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2128:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2129:     for (j=ys; j<ys+ny; j++) {
2130:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2131:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2132:       for (k=zs; k<zs+nz; k++) {
2133:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2134:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

2136:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

2138:         /* Find block columns in block row */
2139:         cnt = 0;
2140:         for (ii=istart; ii<iend+1; ii++) {
2141:           for (jj=jstart; jj<jend+1; jj++) {
2142:             for (kk=kstart; kk<kend+1; kk++) {
2143:               if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2144:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2145:               }
2146:             }
2147:           }
2148:         }
2149:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2150:         MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2151:       }
2152:     }
2153:   }
2154:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2155:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2156:   MatPreallocateFinalize(dnz,onz);

2158:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2160:   /*
2161:     For each node in the grid: we get the neighbors in the local (on processor ordering
2162:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2163:     PETSc ordering.
2164:   */
2165:   if (!da->prealloc_only) {
2166:     PetscCalloc1(col*col*col*nc*nc,&values);
2167:     for (i=xs; i<xs+nx; i++) {
2168:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2169:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2170:       for (j=ys; j<ys+ny; j++) {
2171:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2172:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2173:         for (k=zs; k<zs+nz; k++) {
2174:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2175:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

2177:           slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

2179:           cnt = 0;
2180:           for (ii=istart; ii<iend+1; ii++) {
2181:             for (jj=jstart; jj<jend+1; jj++) {
2182:               for (kk=kstart; kk<kend+1; kk++) {
2183:                 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2184:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2185:                 }
2186:               }
2187:             }
2188:           }
2189:           L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2190:           MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2191:         }
2192:       }
2193:     }
2194:     PetscFree(values);
2195:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2196:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2197:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2198:   }
2199:   PetscFree(cols);
2200:   return(0);
2201: }

2203: /* ---------------------------------------------------------------------------------*/

2205: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2206: {
2207:   PetscErrorCode         ierr;
2208:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2209:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2210:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2211:   DM_DA                  *dd = (DM_DA*)da->data;
2212:   PetscInt               ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2213:   MPI_Comm               comm;
2214:   PetscScalar            *values;
2215:   DMBoundaryType         bx,by,bz;
2216:   ISLocalToGlobalMapping ltog;
2217:   DMDAStencilType        st;
2218:   PetscBool              removedups = PETSC_FALSE;

2221:   /*
2222:          nc - number of components per grid point
2223:          col - number of colors needed in one direction for single component problem

2225:   */
2226:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2227:   col  = 2*s + 1;
2228:   if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
2229:                  by 2*stencil_width + 1\n");
2230:   if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
2231:                  by 2*stencil_width + 1\n");
2232:   if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
2233:                  by 2*stencil_width + 1\n");

2235:   /*
2236:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2237:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2238:   */
2239:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2240:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2241:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

2243:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2244:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2245:   PetscObjectGetComm((PetscObject)da,&comm);

2247:   PetscMalloc1(col*col*col*nc,&cols);
2248:   DMGetLocalToGlobalMapping(da,&ltog);

2250:   /* determine the matrix preallocation information */
2251:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);

2253:   MatSetBlockSize(J,nc);
2254:   for (i=xs; i<xs+nx; i++) {
2255:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2256:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2257:     for (j=ys; j<ys+ny; j++) {
2258:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2259:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2260:       for (k=zs; k<zs+nz; k++) {
2261:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2262:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

2264:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

2266:         for (l=0; l<nc; l++) {
2267:           cnt = 0;
2268:           for (ii=istart; ii<iend+1; ii++) {
2269:             for (jj=jstart; jj<jend+1; jj++) {
2270:               for (kk=kstart; kk<kend+1; kk++) {
2271:                 if (ii || jj || kk) {
2272:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2273:                     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);
2274:                   }
2275:                 } else {
2276:                   if (dfill) {
2277:                     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);
2278:                   } else {
2279:                     for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2280:                   }
2281:                 }
2282:               }
2283:             }
2284:           }
2285:           row  = l + nc*(slot);
2286:           maxcnt = PetscMax(maxcnt,cnt);
2287:           if (removedups) {
2288:             MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2289:           } else {
2290:             MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2291:           }
2292:         }
2293:       }
2294:     }
2295:   }
2296:   MatSeqAIJSetPreallocation(J,0,dnz);
2297:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2298:   MatPreallocateFinalize(dnz,onz);
2299:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2301:   /*
2302:     For each node in the grid: we get the neighbors in the local (on processor ordering
2303:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2304:     PETSc ordering.
2305:   */
2306:   if (!da->prealloc_only) {
2307:     PetscCalloc1(maxcnt,&values);
2308:     for (i=xs; i<xs+nx; i++) {
2309:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2310:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2311:       for (j=ys; j<ys+ny; j++) {
2312:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2313:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2314:         for (k=zs; k<zs+nz; k++) {
2315:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2316:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

2318:           slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);

2320:           for (l=0; l<nc; l++) {
2321:             cnt = 0;
2322:             for (ii=istart; ii<iend+1; ii++) {
2323:               for (jj=jstart; jj<jend+1; jj++) {
2324:                 for (kk=kstart; kk<kend+1; kk++) {
2325:                   if (ii || jj || kk) {
2326:                     if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2327:                       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);
2328:                     }
2329:                   } else {
2330:                     if (dfill) {
2331:                       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);
2332:                     } else {
2333:                       for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2334:                     }
2335:                   }
2336:                 }
2337:               }
2338:             }
2339:             row  = l + nc*(slot);
2340:             MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2341:           }
2342:         }
2343:       }
2344:     }
2345:     PetscFree(values);
2346:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2347:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2348:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2349:   }
2350:   PetscFree(cols);
2351:   return(0);
2352: }