Actual source code: fdda.c

petsc-3.11.4 2019-09-28
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:   MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1580:   MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);

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

1585:   */
1586:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1587:   if (s > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Matrix creation for 1d not implemented correctly for stencil width larger than 1");
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:       if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1605:         if (size > 1) ocols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1606:         else cols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1607:       }
1608:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1609:       cnt++;
1610:     }
1611:   }
1612:   for (i=s; i<nx-s; i++) {
1613:     for (j=0; j<nc; j++) {
1614:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1615:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1616:       cnt++;
1617:     }
1618:   }
1619:   /* coupling with process to the right */
1620:   for (i=nx-s; i<nx; i++) {
1621:     for (j=0; j<nc; j++) {
1622:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1623:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1624:       if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1625:         if (size > 1) ocols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1626:         else cols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1627:       }
1628:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1629:       cnt++;
1630:     }
1631:   }

1633:   MatSeqAIJSetPreallocation(J,0,cols);
1634:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1635:   PetscFree2(cols,ocols);

1637:   DMGetLocalToGlobalMapping(da,&ltog);
1638:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1640:   /*
1641:     For each node in the grid: we get the neighbors in the local (on processor ordering
1642:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1643:     PETSc ordering.
1644:   */
1645:   if (!da->prealloc_only) {
1646:     PetscCalloc2(maxcnt,&values,maxcnt,&cols);

1648:     row = xs*nc;
1649:     /* coupling with process to the left */
1650:     for (i=xs; i<xs+s; i++) {
1651:       for (j=0; j<nc; j++) {
1652:         cnt = 0;
1653:         if (rank) {
1654:           for (l=0; l<s; l++) {
1655:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1656:           }
1657:         }
1658:         if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1659:           for (l=0; l<s; l++) {
1660:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (m + i - s - l)*nc + ofill[k];
1661:           }
1662:         }
1663:         if (dfill) {
1664:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1665:             cols[cnt++] = i*nc + dfill[k];
1666:           }
1667:         } else {
1668:           for (k=0; k<nc; k++) {
1669:             cols[cnt++] = i*nc + k;
1670:           }
1671:         }
1672:         for (l=0; l<s; l++) {
1673:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1674:         }
1675:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1676:         row++;
1677:       }
1678:     }
1679:     for (i=xs+s; i<xs+nx-s; i++) {
1680:       for (j=0; j<nc; j++) {
1681:         cnt = 0;
1682:         for (l=0; l<s; l++) {
1683:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1684:         }
1685:         if (dfill) {
1686:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1687:             cols[cnt++] = i*nc + dfill[k];
1688:           }
1689:         } else {
1690:           for (k=0; k<nc; k++) {
1691:             cols[cnt++] = i*nc + k;
1692:           }
1693:         }
1694:         for (l=0; l<s; l++) {
1695:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1696:         }
1697:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1698:         row++;
1699:       }
1700:     }
1701:     /* coupling with process to the right */
1702:     for (i=xs+nx-s; i<xs+nx; i++) {
1703:       for (j=0; j<nc; j++) {
1704:         cnt = 0;
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:         if (dfill) {
1709:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1710:             cols[cnt++] = i*nc + dfill[k];
1711:           }
1712:         } else {
1713:           for (k=0; k<nc; k++) {
1714:             cols[cnt++] = i*nc + k;
1715:           }
1716:         }
1717:         if (rank < size-1) {
1718:           for (l=0; l<s; l++) {
1719:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1720:           }
1721:         }
1722:         if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1723:           for (l=0; l<s; l++) {
1724:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s - l - m + 2)*nc + ofill[k];
1725:           }
1726:         }
1727:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1728:         row++;
1729:       }
1730:     }
1731:     PetscFree2(values,cols);
1732:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1733:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1734:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1735:   }
1736:   return(0);
1737: }

1739: /* ---------------------------------------------------------------------------------*/

1741: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1742: {
1743:   PetscErrorCode         ierr;
1744:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1745:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1746:   PetscInt               istart,iend;
1747:   PetscScalar            *values;
1748:   DMBoundaryType         bx;
1749:   ISLocalToGlobalMapping ltog,mltog;

1752:   /*
1753:          nc - number of components per grid point
1754:          col - number of colors needed in one direction for single component problem

1756:   */
1757:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1758:   col  = 2*s + 1;

1760:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1761:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1763:   MatSetBlockSize(J,nc);
1764:   MatSeqAIJSetPreallocation(J,col*nc,0);
1765:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1767:   DMGetLocalToGlobalMapping(da,&ltog);
1768:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1769:   if (!mltog) {
1770:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1771:   }

1773:   /*
1774:     For each node in the grid: we get the neighbors in the local (on processor ordering
1775:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1776:     PETSc ordering.
1777:   */
1778:   if (!da->prealloc_only) {
1779:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1780:     PetscCalloc1(col*nc*nc,&values);
1781:     for (i=xs; i<xs+nx; i++) {
1782:       istart = PetscMax(-s,gxs - i);
1783:       iend   = PetscMin(s,gxs + gnx - i - 1);
1784:       slot   = i - gxs;

1786:       cnt = 0;
1787:       for (l=0; l<nc; l++) {
1788:         for (i1=istart; i1<iend+1; i1++) {
1789:           cols[cnt++] = l + nc*(slot + i1);
1790:         }
1791:         rows[l] = l + nc*(slot);
1792:       }
1793:       MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1794:     }
1795:     PetscFree(values);
1796:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1797:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1798:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1799:     PetscFree2(rows,cols);
1800:   }
1801:   return(0);
1802: }

1804: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1805: {
1806:   PetscErrorCode         ierr;
1807:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1808:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1809:   PetscInt               istart,iend,jstart,jend,ii,jj;
1810:   MPI_Comm               comm;
1811:   PetscScalar            *values;
1812:   DMBoundaryType         bx,by;
1813:   DMDAStencilType        st;
1814:   ISLocalToGlobalMapping ltog;

1817:   /*
1818:      nc - number of components per grid point
1819:      col - number of colors needed in one direction for single component problem
1820:   */
1821:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1822:   col  = 2*s + 1;

1824:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1825:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1826:   PetscObjectGetComm((PetscObject)da,&comm);

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

1830:   DMGetLocalToGlobalMapping(da,&ltog);

1832:   /* determine the matrix preallocation information */
1833:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1834:   for (i=xs; i<xs+nx; i++) {
1835:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1836:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1837:     for (j=ys; j<ys+ny; j++) {
1838:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1839:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1840:       slot   = i - gxs + gnx*(j - gys);

1842:       /* Find block columns in block row */
1843:       cnt = 0;
1844:       for (ii=istart; ii<iend+1; ii++) {
1845:         for (jj=jstart; jj<jend+1; jj++) {
1846:           if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1847:             cols[cnt++] = slot + ii + gnx*jj;
1848:           }
1849:         }
1850:       }
1851:       MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1852:     }
1853:   }
1854:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1855:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1856:   MatPreallocateFinalize(dnz,onz);

1858:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1860:   /*
1861:     For each node in the grid: we get the neighbors in the local (on processor ordering
1862:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1863:     PETSc ordering.
1864:   */
1865:   if (!da->prealloc_only) {
1866:     PetscCalloc1(col*col*nc*nc,&values);
1867:     for (i=xs; i<xs+nx; i++) {
1868:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1869:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1870:       for (j=ys; j<ys+ny; j++) {
1871:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1872:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1873:         slot = i - gxs + gnx*(j - gys);
1874:         cnt  = 0;
1875:         for (ii=istart; ii<iend+1; ii++) {
1876:           for (jj=jstart; jj<jend+1; jj++) {
1877:             if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1878:               cols[cnt++] = slot + ii + gnx*jj;
1879:             }
1880:           }
1881:         }
1882:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1883:       }
1884:     }
1885:     PetscFree(values);
1886:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1887:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1888:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1889:   }
1890:   PetscFree(cols);
1891:   return(0);
1892: }

1894: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1895: {
1896:   PetscErrorCode         ierr;
1897:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1898:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1899:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1900:   MPI_Comm               comm;
1901:   PetscScalar            *values;
1902:   DMBoundaryType         bx,by,bz;
1903:   DMDAStencilType        st;
1904:   ISLocalToGlobalMapping ltog;

1907:   /*
1908:          nc - number of components per grid point
1909:          col - number of colors needed in one direction for single component problem

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

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

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

1921:   DMGetLocalToGlobalMapping(da,&ltog);

1923:   /* determine the matrix preallocation information */
1924:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1925:   for (i=xs; i<xs+nx; i++) {
1926:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1927:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1928:     for (j=ys; j<ys+ny; j++) {
1929:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1930:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1931:       for (k=zs; k<zs+nz; k++) {
1932:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1933:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1937:         /* Find block columns in block row */
1938:         cnt = 0;
1939:         for (ii=istart; ii<iend+1; ii++) {
1940:           for (jj=jstart; jj<jend+1; jj++) {
1941:             for (kk=kstart; kk<kend+1; kk++) {
1942:               if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1943:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1944:               }
1945:             }
1946:           }
1947:         }
1948:         MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1949:       }
1950:     }
1951:   }
1952:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1953:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1954:   MatPreallocateFinalize(dnz,onz);

1956:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1958:   /*
1959:     For each node in the grid: we get the neighbors in the local (on processor ordering
1960:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1961:     PETSc ordering.
1962:   */
1963:   if (!da->prealloc_only) {
1964:     PetscCalloc1(col*col*col*nc*nc,&values);
1965:     for (i=xs; i<xs+nx; i++) {
1966:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1967:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1968:       for (j=ys; j<ys+ny; j++) {
1969:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1970:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1971:         for (k=zs; k<zs+nz; k++) {
1972:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1973:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1977:           cnt = 0;
1978:           for (ii=istart; ii<iend+1; ii++) {
1979:             for (jj=jstart; jj<jend+1; jj++) {
1980:               for (kk=kstart; kk<kend+1; kk++) {
1981:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1982:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1983:                 }
1984:               }
1985:             }
1986:           }
1987:           MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1988:         }
1989:       }
1990:     }
1991:     PetscFree(values);
1992:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1993:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1994:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1995:   }
1996:   PetscFree(cols);
1997:   return(0);
1998: }

2000: /*
2001:   This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
2002:   identify in the local ordering with periodic domain.
2003: */
2004: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
2005: {
2007:   PetscInt       i,n;

2010:   ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
2011:   ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
2012:   for (i=0,n=0; i<*cnt; i++) {
2013:     if (col[i] >= *row) col[n++] = col[i];
2014:   }
2015:   *cnt = n;
2016:   return(0);
2017: }

2019: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
2020: {
2021:   PetscErrorCode         ierr;
2022:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2023:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
2024:   PetscInt               istart,iend,jstart,jend,ii,jj;
2025:   MPI_Comm               comm;
2026:   PetscScalar            *values;
2027:   DMBoundaryType         bx,by;
2028:   DMDAStencilType        st;
2029:   ISLocalToGlobalMapping ltog;

2032:   /*
2033:      nc - number of components per grid point
2034:      col - number of colors needed in one direction for single component problem
2035:   */
2036:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
2037:   col  = 2*s + 1;

2039:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
2040:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
2041:   PetscObjectGetComm((PetscObject)da,&comm);

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

2045:   DMGetLocalToGlobalMapping(da,&ltog);

2047:   /* determine the matrix preallocation information */
2048:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
2049:   for (i=xs; i<xs+nx; i++) {
2050:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2051:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2052:     for (j=ys; j<ys+ny; j++) {
2053:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2054:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2055:       slot   = i - gxs + gnx*(j - gys);

2057:       /* Find block columns in block row */
2058:       cnt = 0;
2059:       for (ii=istart; ii<iend+1; ii++) {
2060:         for (jj=jstart; jj<jend+1; jj++) {
2061:           if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2062:             cols[cnt++] = slot + ii + gnx*jj;
2063:           }
2064:         }
2065:       }
2066:       L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2067:       MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2068:     }
2069:   }
2070:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2071:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2072:   MatPreallocateFinalize(dnz,onz);

2074:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2076:   /*
2077:     For each node in the grid: we get the neighbors in the local (on processor ordering
2078:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2079:     PETSc ordering.
2080:   */
2081:   if (!da->prealloc_only) {
2082:     PetscCalloc1(col*col*nc*nc,&values);
2083:     for (i=xs; i<xs+nx; i++) {
2084:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2085:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2086:       for (j=ys; j<ys+ny; j++) {
2087:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2088:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2089:         slot   = i - gxs + gnx*(j - gys);

2091:         /* Find block columns in block row */
2092:         cnt = 0;
2093:         for (ii=istart; ii<iend+1; ii++) {
2094:           for (jj=jstart; jj<jend+1; jj++) {
2095:             if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2096:               cols[cnt++] = slot + ii + gnx*jj;
2097:             }
2098:           }
2099:         }
2100:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2101:         MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2102:       }
2103:     }
2104:     PetscFree(values);
2105:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2106:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2107:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2108:   }
2109:   PetscFree(cols);
2110:   return(0);
2111: }

2113: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
2114: {
2115:   PetscErrorCode         ierr;
2116:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2117:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
2118:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
2119:   MPI_Comm               comm;
2120:   PetscScalar            *values;
2121:   DMBoundaryType         bx,by,bz;
2122:   DMDAStencilType        st;
2123:   ISLocalToGlobalMapping ltog;

2126:   /*
2127:      nc - number of components per grid point
2128:      col - number of colors needed in one direction for single component problem
2129:   */
2130:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2131:   col  = 2*s + 1;

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

2137:   /* create the matrix */
2138:   PetscMalloc1(col*col*col,&cols);

2140:   DMGetLocalToGlobalMapping(da,&ltog);

2142:   /* determine the matrix preallocation information */
2143:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2144:   for (i=xs; i<xs+nx; i++) {
2145:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2146:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2147:     for (j=ys; j<ys+ny; j++) {
2148:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2149:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2150:       for (k=zs; k<zs+nz; k++) {
2151:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2152:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2156:         /* Find block columns in block row */
2157:         cnt = 0;
2158:         for (ii=istart; ii<iend+1; ii++) {
2159:           for (jj=jstart; jj<jend+1; jj++) {
2160:             for (kk=kstart; kk<kend+1; kk++) {
2161:               if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2162:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2163:               }
2164:             }
2165:           }
2166:         }
2167:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2168:         MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2169:       }
2170:     }
2171:   }
2172:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2173:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2174:   MatPreallocateFinalize(dnz,onz);

2176:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2178:   /*
2179:     For each node in the grid: we get the neighbors in the local (on processor ordering
2180:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2181:     PETSc ordering.
2182:   */
2183:   if (!da->prealloc_only) {
2184:     PetscCalloc1(col*col*col*nc*nc,&values);
2185:     for (i=xs; i<xs+nx; i++) {
2186:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2187:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2188:       for (j=ys; j<ys+ny; j++) {
2189:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2190:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2191:         for (k=zs; k<zs+nz; k++) {
2192:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2193:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2197:           cnt = 0;
2198:           for (ii=istart; ii<iend+1; ii++) {
2199:             for (jj=jstart; jj<jend+1; jj++) {
2200:               for (kk=kstart; kk<kend+1; kk++) {
2201:                 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2202:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2203:                 }
2204:               }
2205:             }
2206:           }
2207:           L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2208:           MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2209:         }
2210:       }
2211:     }
2212:     PetscFree(values);
2213:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2214:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2215:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2216:   }
2217:   PetscFree(cols);
2218:   return(0);
2219: }

2221: /* ---------------------------------------------------------------------------------*/

2223: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2224: {
2225:   PetscErrorCode         ierr;
2226:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2227:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2228:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2229:   DM_DA                  *dd = (DM_DA*)da->data;
2230:   PetscInt               ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2231:   MPI_Comm               comm;
2232:   PetscScalar            *values;
2233:   DMBoundaryType         bx,by,bz;
2234:   ISLocalToGlobalMapping ltog;
2235:   DMDAStencilType        st;
2236:   PetscBool              removedups = PETSC_FALSE;

2239:   /*
2240:          nc - number of components per grid point
2241:          col - number of colors needed in one direction for single component problem

2243:   */
2244:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2245:   col  = 2*s + 1;
2246:   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\
2247:                  by 2*stencil_width + 1\n");
2248:   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\
2249:                  by 2*stencil_width + 1\n");
2250:   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\
2251:                  by 2*stencil_width + 1\n");

2253:   /*
2254:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2255:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2256:   */
2257:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2258:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2259:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

2265:   PetscMalloc1(col*col*col*nc,&cols);
2266:   DMGetLocalToGlobalMapping(da,&ltog);

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

2271:   MatSetBlockSize(J,nc);
2272:   for (i=xs; i<xs+nx; i++) {
2273:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2274:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2275:     for (j=ys; j<ys+ny; j++) {
2276:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2277:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2278:       for (k=zs; k<zs+nz; k++) {
2279:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2280:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2284:         for (l=0; l<nc; l++) {
2285:           cnt = 0;
2286:           for (ii=istart; ii<iend+1; ii++) {
2287:             for (jj=jstart; jj<jend+1; jj++) {
2288:               for (kk=kstart; kk<kend+1; kk++) {
2289:                 if (ii || jj || kk) {
2290:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2291:                     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);
2292:                   }
2293:                 } else {
2294:                   if (dfill) {
2295:                     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);
2296:                   } else {
2297:                     for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2298:                   }
2299:                 }
2300:               }
2301:             }
2302:           }
2303:           row  = l + nc*(slot);
2304:           maxcnt = PetscMax(maxcnt,cnt);
2305:           if (removedups) {
2306:             MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2307:           } else {
2308:             MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2309:           }
2310:         }
2311:       }
2312:     }
2313:   }
2314:   MatSeqAIJSetPreallocation(J,0,dnz);
2315:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2316:   MatPreallocateFinalize(dnz,onz);
2317:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2319:   /*
2320:     For each node in the grid: we get the neighbors in the local (on processor ordering
2321:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2322:     PETSc ordering.
2323:   */
2324:   if (!da->prealloc_only) {
2325:     PetscCalloc1(maxcnt,&values);
2326:     for (i=xs; i<xs+nx; i++) {
2327:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2328:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2329:       for (j=ys; j<ys+ny; j++) {
2330:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2331:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2332:         for (k=zs; k<zs+nz; k++) {
2333:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2334:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2338:           for (l=0; l<nc; l++) {
2339:             cnt = 0;
2340:             for (ii=istart; ii<iend+1; ii++) {
2341:               for (jj=jstart; jj<jend+1; jj++) {
2342:                 for (kk=kstart; kk<kend+1; kk++) {
2343:                   if (ii || jj || kk) {
2344:                     if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2345:                       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);
2346:                     }
2347:                   } else {
2348:                     if (dfill) {
2349:                       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);
2350:                     } else {
2351:                       for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2352:                     }
2353:                   }
2354:                 }
2355:               }
2356:             }
2357:             row  = l + nc*(slot);
2358:             MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2359:           }
2360:         }
2361:       }
2362:     }
2363:     PetscFree(values);
2364:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2365:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2366:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2367:   }
2368:   PetscFree(cols);
2369:   return(0);
2370: }