Actual source code: fdda.c

petsc-3.8.0 2017-09-26
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: }

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

 55:     Logically Collective on DMDA

 57:     Input Parameter:
 58: +   da - the distributed array
 59: .   dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
 60: -   ofill - the fill pattern in the off-diagonal blocks


 63:     Level: developer

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

 68:            The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
 69:        representing coupling and 0 entries for missing coupling. For example
 70: $             dfill[9] = {1, 0, 0,
 71: $                         1, 1, 0,
 72: $                         0, 1, 1}
 73:        means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
 74:        itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
 75:        diagonal block).

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

 80:    Contributed by Glenn Hammond

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

 84: @*/
 85: PetscErrorCode  DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
 86: {
 87:   DM_DA          *dd = (DM_DA*)da->data;
 89:   PetscInt       i,k,cnt = 1;

 92:   DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
 93:   DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);

 95:   /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of
 96:    columns to the left with any nonzeros in them plus 1 */
 97:   PetscCalloc1(dd->w,&dd->ofillcols);
 98:   for (i=0; i<dd->w; i++) {
 99:     for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
100:   }
101:   for (i=0; i<dd->w; i++) {
102:     if (dd->ofillcols[i]) {
103:       dd->ofillcols[i] = cnt++;
104:     }
105:   }
106:   return(0);
107: }


110: PetscErrorCode  DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
111: {
112:   PetscErrorCode   ierr;
113:   PetscInt         dim,m,n,p,nc;
114:   DMBoundaryType bx,by,bz;
115:   MPI_Comm         comm;
116:   PetscMPIInt      size;
117:   PetscBool        isBAIJ;
118:   DM_DA            *dd = (DM_DA*)da->data;

121:   /*
122:                                   m
123:           ------------------------------------------------------
124:          |                                                     |
125:          |                                                     |
126:          |               ----------------------                |
127:          |               |                    |                |
128:       n  |           yn  |                    |                |
129:          |               |                    |                |
130:          |               .---------------------                |
131:          |             (xs,ys)     xn                          |
132:          |            .                                        |
133:          |         (gxs,gys)                                   |
134:          |                                                     |
135:           -----------------------------------------------------
136:   */

138:   /*
139:          nc - number of components per grid point
140:          col - number of colors needed in one direction for single component problem

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

145:   PetscObjectGetComm((PetscObject)da,&comm);
146:   MPI_Comm_size(comm,&size);
147:   if (ctype == IS_COLORING_LOCAL) {
148:     if (size == 1) {
149:       ctype = IS_COLORING_GLOBAL;
150:     } else if (dim > 1) {
151:       if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
152:         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");
153:       }
154:     }
155:   }

157:   /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
158:      matrices is for the blocks, not the individual matrix elements  */
159:   PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);
160:   if (!isBAIJ) {PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);}
161:   if (!isBAIJ) {PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);}
162:   if (isBAIJ) {
163:     dd->w  = 1;
164:     dd->xs = dd->xs/nc;
165:     dd->xe = dd->xe/nc;
166:     dd->Xs = dd->Xs/nc;
167:     dd->Xe = dd->Xe/nc;
168:   }

170:   /*
171:      We do not provide a getcoloring function in the DMDA operations because
172:    the basic DMDA does not know about matrices. We think of DMDA as being more
173:    more low-level then matrices.
174:   */
175:   if (dim == 1) {
176:     DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
177:   } else if (dim == 2) {
178:      DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
179:   } else if (dim == 3) {
180:      DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
181:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
182:   if (isBAIJ) {
183:     dd->w  = nc;
184:     dd->xs = dd->xs*nc;
185:     dd->xe = dd->xe*nc;
186:     dd->Xs = dd->Xs*nc;
187:     dd->Xe = dd->Xe*nc;
188:   }
189:   return(0);
190: }

192: /* ---------------------------------------------------------------------------------*/

194: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
195: {
196:   PetscErrorCode   ierr;
197:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
198:   PetscInt         ncolors;
199:   MPI_Comm         comm;
200:   DMBoundaryType bx,by;
201:   DMDAStencilType  st;
202:   ISColoringValue  *colors;
203:   DM_DA            *dd = (DM_DA*)da->data;

206:   /*
207:          nc - number of components per grid point
208:          col - number of colors needed in one direction for single component problem

210:   */
211:   DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
212:   col  = 2*s + 1;
213:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
214:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
215:   PetscObjectGetComm((PetscObject)da,&comm);

217:   /* special case as taught to us by Paul Hovland */
218:   if (st == DMDA_STENCIL_STAR && s == 1) {
219:     DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
220:   } else {

222:     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\
223:                                                             by 2*stencil_width + 1 (%d)\n", m, col);
224:     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\
225:                                                             by 2*stencil_width + 1 (%d)\n", n, col);
226:     if (ctype == IS_COLORING_GLOBAL) {
227:       if (!dd->localcoloring) {
228:         PetscMalloc1(nc*nx*ny,&colors);
229:         ii   = 0;
230:         for (j=ys; j<ys+ny; j++) {
231:           for (i=xs; i<xs+nx; i++) {
232:             for (k=0; k<nc; k++) {
233:               colors[ii++] = k + nc*((i % col) + col*(j % col));
234:             }
235:           }
236:         }
237:         ncolors = nc + nc*(col-1 + col*(col-1));
238:         ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
239:       }
240:       *coloring = dd->localcoloring;
241:     } else if (ctype == IS_COLORING_LOCAL) {
242:       if (!dd->ghostedcoloring) {
243:         PetscMalloc1(nc*gnx*gny,&colors);
244:         ii   = 0;
245:         for (j=gys; j<gys+gny; j++) {
246:           for (i=gxs; i<gxs+gnx; i++) {
247:             for (k=0; k<nc; k++) {
248:               /* the complicated stuff is to handle periodic boundaries */
249:               colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
250:             }
251:           }
252:         }
253:         ncolors = nc + nc*(col - 1 + col*(col-1));
254:         ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
255:         /* PetscIntView(ncolors,(PetscInt*)colors,0); */

257:         ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
258:       }
259:       *coloring = dd->ghostedcoloring;
260:     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
261:   }
262:   ISColoringReference(*coloring);
263:   return(0);
264: }

266: /* ---------------------------------------------------------------------------------*/

268: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
269: {
270:   PetscErrorCode   ierr;
271:   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;
272:   PetscInt         ncolors;
273:   MPI_Comm         comm;
274:   DMBoundaryType bx,by,bz;
275:   DMDAStencilType  st;
276:   ISColoringValue  *colors;
277:   DM_DA            *dd = (DM_DA*)da->data;

280:   /*
281:          nc - number of components per grid point
282:          col - number of colors needed in one direction for single component problem

284:   */
285:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
286:   col  = 2*s + 1;
287:   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\
288:                                                          by 2*stencil_width + 1\n");
289:   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\
290:                                                          by 2*stencil_width + 1\n");
291:   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\
292:                                                          by 2*stencil_width + 1\n");

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

298:   /* create the coloring */
299:   if (ctype == IS_COLORING_GLOBAL) {
300:     if (!dd->localcoloring) {
301:       PetscMalloc1(nc*nx*ny*nz,&colors);
302:       ii   = 0;
303:       for (k=zs; k<zs+nz; k++) {
304:         for (j=ys; j<ys+ny; j++) {
305:           for (i=xs; i<xs+nx; i++) {
306:             for (l=0; l<nc; l++) {
307:               colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
308:             }
309:           }
310:         }
311:       }
312:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
313:       ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
314:     }
315:     *coloring = dd->localcoloring;
316:   } else if (ctype == IS_COLORING_LOCAL) {
317:     if (!dd->ghostedcoloring) {
318:       PetscMalloc1(nc*gnx*gny*gnz,&colors);
319:       ii   = 0;
320:       for (k=gzs; k<gzs+gnz; k++) {
321:         for (j=gys; j<gys+gny; j++) {
322:           for (i=gxs; i<gxs+gnx; i++) {
323:             for (l=0; l<nc; l++) {
324:               /* the complicated stuff is to handle periodic boundaries */
325:               colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
326:             }
327:           }
328:         }
329:       }
330:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
331:       ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
332:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
333:     }
334:     *coloring = dd->ghostedcoloring;
335:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
336:   ISColoringReference(*coloring);
337:   return(0);
338: }

340: /* ---------------------------------------------------------------------------------*/

342: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
343: {
344:   PetscErrorCode   ierr;
345:   PetscInt         xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
346:   PetscInt         ncolors;
347:   MPI_Comm         comm;
348:   DMBoundaryType bx;
349:   ISColoringValue  *colors;
350:   DM_DA            *dd = (DM_DA*)da->data;

353:   /*
354:          nc - number of components per grid point
355:          col - number of colors needed in one direction for single component problem

357:   */
358:   DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
359:   col  = 2*s + 1;

361:   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\
362:                                                           by 2*stencil_width + 1 %d\n",(int)m,(int)col);

364:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
365:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
366:   PetscObjectGetComm((PetscObject)da,&comm);

368:   /* create the coloring */
369:   if (ctype == IS_COLORING_GLOBAL) {
370:     if (!dd->localcoloring) {
371:       PetscMalloc1(nc*nx,&colors);
372:       if (dd->ofillcols) {
373:         PetscInt tc = 0;
374:         for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
375:         i1 = 0;
376:         for (i=xs; i<xs+nx; i++) {
377:           for (l=0; l<nc; l++) {
378:             if (dd->ofillcols[l] && (i % col)) {
379:               colors[i1++] =  nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
380:             } else {
381:               colors[i1++] = l;
382:             }
383:           }
384:         }
385:         ncolors = nc + 2*s*tc;
386:       } else {
387:         i1 = 0;
388:         for (i=xs; i<xs+nx; i++) {
389:           for (l=0; l<nc; l++) {
390:             colors[i1++] = l + nc*(i % col);
391:           }
392:         }
393:         ncolors = nc + nc*(col-1);
394:       }
395:       ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
396:     }
397:     *coloring = dd->localcoloring;
398:   } else if (ctype == IS_COLORING_LOCAL) {
399:     if (!dd->ghostedcoloring) {
400:       PetscMalloc1(nc*gnx,&colors);
401:       i1   = 0;
402:       for (i=gxs; i<gxs+gnx; i++) {
403:         for (l=0; l<nc; l++) {
404:           /* the complicated stuff is to handle periodic boundaries */
405:           colors[i1++] = l + nc*(SetInRange(i,m) % col);
406:         }
407:       }
408:       ncolors = nc + nc*(col-1);
409:       ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
410:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
411:     }
412:     *coloring = dd->ghostedcoloring;
413:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
414:   ISColoringReference(*coloring);
415:   return(0);
416: }

418: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
419: {
420:   PetscErrorCode   ierr;
421:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
422:   PetscInt         ncolors;
423:   MPI_Comm         comm;
424:   DMBoundaryType bx,by;
425:   ISColoringValue  *colors;
426:   DM_DA            *dd = (DM_DA*)da->data;

429:   /*
430:          nc - number of components per grid point
431:          col - number of colors needed in one direction for single component problem

433:   */
434:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
435:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
436:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
437:   PetscObjectGetComm((PetscObject)da,&comm);

439:   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");
440:   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");

442:   /* create the coloring */
443:   if (ctype == IS_COLORING_GLOBAL) {
444:     if (!dd->localcoloring) {
445:       PetscMalloc1(nc*nx*ny,&colors);
446:       ii   = 0;
447:       for (j=ys; j<ys+ny; j++) {
448:         for (i=xs; i<xs+nx; i++) {
449:           for (k=0; k<nc; k++) {
450:             colors[ii++] = k + nc*((3*j+i) % 5);
451:           }
452:         }
453:       }
454:       ncolors = 5*nc;
455:       ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
456:     }
457:     *coloring = dd->localcoloring;
458:   } else if (ctype == IS_COLORING_LOCAL) {
459:     if (!dd->ghostedcoloring) {
460:       PetscMalloc1(nc*gnx*gny,&colors);
461:       ii = 0;
462:       for (j=gys; j<gys+gny; j++) {
463:         for (i=gxs; i<gxs+gnx; i++) {
464:           for (k=0; k<nc; k++) {
465:             colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
466:           }
467:         }
468:       }
469:       ncolors = 5*nc;
470:       ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
471:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
472:     }
473:     *coloring = dd->ghostedcoloring;
474:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
475:   return(0);
476: }

478: /* =========================================================================== */
479: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
480: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
481: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
482: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
483: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
484: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
485: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
486: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
487: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
488: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);

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

493:    Logically Collective on Mat

495:    Input Parameters:
496: +  mat - the matrix
497: -  da - the da

499:    Level: intermediate

501: @*/
502: PetscErrorCode MatSetupDM(Mat mat,DM da)
503: {

509:   PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
510:   return(0);
511: }

513: PetscErrorCode  MatView_MPI_DA(Mat A,PetscViewer viewer)
514: {
515:   DM                da;
516:   PetscErrorCode    ierr;
517:   const char        *prefix;
518:   Mat               Anatural;
519:   AO                ao;
520:   PetscInt          rstart,rend,*petsc,i;
521:   IS                is;
522:   MPI_Comm          comm;
523:   PetscViewerFormat format;

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

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

534:   DMDAGetAO(da,&ao);
535:   MatGetOwnershipRange(A,&rstart,&rend);
536:   PetscMalloc1(rend-rstart,&petsc);
537:   for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
538:   AOApplicationToPetsc(ao,rend-rstart,petsc);
539:   ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);

541:   /* call viewer on natural ordering */
542:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
543:   ISDestroy(&is);
544:   PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
545:   PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
546:   PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
547:   MatView(Anatural,viewer);
548:   MatDestroy(&Anatural);
549:   return(0);
550: }

552: PetscErrorCode  MatLoad_MPI_DA(Mat A,PetscViewer viewer)
553: {
554:   DM             da;
556:   Mat            Anatural,Aapp;
557:   AO             ao;
558:   PetscInt       rstart,rend,*app,i,m,n,M,N;
559:   IS             is;
560:   MPI_Comm       comm;

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

567:   /* Load the matrix in natural ordering */
568:   MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
569:   MatSetType(Anatural,((PetscObject)A)->type_name);
570:   MatGetSize(A,&M,&N);
571:   MatGetLocalSize(A,&m,&n);
572:   MatSetSizes(Anatural,m,n,M,N);
573:   MatLoad(Anatural,viewer);

575:   /* Map natural ordering to application ordering and create IS */
576:   DMDAGetAO(da,&ao);
577:   MatGetOwnershipRange(Anatural,&rstart,&rend);
578:   PetscMalloc1(rend-rstart,&app);
579:   for (i=rstart; i<rend; i++) app[i-rstart] = i;
580:   AOPetscToApplication(ao,rend-rstart,app);
581:   ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);

583:   /* Do permutation and replace header */
584:   MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
585:   MatHeaderReplace(A,&Aapp);
586:   ISDestroy(&is);
587:   MatDestroy(&Anatural);
588:   return(0);
589: }

591: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
592: {
594:   PetscInt       dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
595:   Mat            A;
596:   MPI_Comm       comm;
597:   MatType        Atype;
598:   PetscSection   section, sectionGlobal;
599:   void           (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL;
600:   MatType        mtype;
601:   PetscMPIInt    size;
602:   DM_DA          *dd = (DM_DA*)da->data;

605:   MatInitializePackage();
606:   mtype = da->mattype;

608:   DMGetDefaultSection(da, &section);
609:   if (section) {
610:     PetscInt  bs = -1;
611:     PetscInt  localSize;
612:     PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;

614:     DMGetDefaultGlobalSection(da, &sectionGlobal);
615:     PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
616:     MatCreate(PetscObjectComm((PetscObject)da),&A);
617:     MatSetSizes(A,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
618:     MatSetType(A,mtype);
619:     PetscStrcmp(mtype,MATSHELL,&isShell);
620:     PetscStrcmp(mtype,MATBAIJ,&isBlock);
621:     PetscStrcmp(mtype,MATSEQBAIJ,&isSeqBlock);
622:     PetscStrcmp(mtype,MATMPIBAIJ,&isMPIBlock);
623:     PetscStrcmp(mtype,MATSBAIJ,&isSymBlock);
624:     PetscStrcmp(mtype,MATSEQSBAIJ,&isSymSeqBlock);
625:     PetscStrcmp(mtype,MATMPISBAIJ,&isSymMPIBlock);
626:     /* Check for symmetric storage */
627:     isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
628:     if (isSymmetric) {
629:       MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
630:     }
631:     if (!isShell) {
632:       PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;

634:       if (bs < 0) {
635:         if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
636:           PetscInt pStart, pEnd, p, dof;

638:           PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
639:           for (p = pStart; p < pEnd; ++p) {
640:             PetscSectionGetDof(sectionGlobal, p, &dof);
641:             if (dof) {
642:               bs = dof;
643:               break;
644:             }
645:           }
646:         } else {
647:           bs = 1;
648:         }
649:         /* Must have same blocksize on all procs (some might have no points) */
650:         bsLocal = bs;
651:         MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
652:       }
653:       PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
654:       /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
655:       PetscFree4(dnz, onz, dnzu, onzu);
656:     }
657:   }
658:   /*
659:                                   m
660:           ------------------------------------------------------
661:          |                                                     |
662:          |                                                     |
663:          |               ----------------------                |
664:          |               |                    |                |
665:       n  |           ny  |                    |                |
666:          |               |                    |                |
667:          |               .---------------------                |
668:          |             (xs,ys)     nx                          |
669:          |            .                                        |
670:          |         (gxs,gys)                                   |
671:          |                                                     |
672:           -----------------------------------------------------
673:   */

675:   /*
676:          nc - number of components per grid point
677:          col - number of colors needed in one direction for single component problem

679:   */
680:   M   = dd->M;
681:   N   = dd->N;
682:   P   = dd->P;
683:   dim = da->dim;
684:   dof = dd->w;
685:   /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
686:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
687:   PetscObjectGetComm((PetscObject)da,&comm);
688:   MatCreate(comm,&A);
689:   MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
690:   MatSetType(A,mtype);
691:   MatSetDM(A,da);
692:   if (da->structure_only) {
693:     MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);
694:   }
695:   MatGetType(A,&Atype);
696:   /*
697:      We do not provide a getmatrix function in the DMDA operations because
698:    the basic DMDA does not know about matrices. We think of DMDA as being more
699:    more low-level than matrices. This is kind of cheating but, cause sometimes
700:    we think of DMDA has higher level than matrices.

702:      We could switch based on Atype (or mtype), but we do not since the
703:    specialized setting routines depend only the particular preallocation
704:    details of the matrix, not the type itself.
705:   */
706:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
707:   if (!aij) {
708:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
709:   }
710:   if (!aij) {
711:     PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
712:     if (!baij) {
713:       PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
714:     }
715:     if (!baij) {
716:       PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
717:       if (!sbaij) {
718:         PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
719:       }
720:     }
721:   }
722:   if (aij) {
723:     if (dim == 1) {
724:       if (dd->ofill) {
725:         DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
726:       } else {
727:         DMCreateMatrix_DA_1d_MPIAIJ(da,A);
728:       }
729:     } else if (dim == 2) {
730:       if (dd->ofill) {
731:         DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
732:       } else {
733:         DMCreateMatrix_DA_2d_MPIAIJ(da,A);
734:       }
735:     } else if (dim == 3) {
736:       if (dd->ofill) {
737:         DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
738:       } else {
739:         DMCreateMatrix_DA_3d_MPIAIJ(da,A);
740:       }
741:     }
742:   } else if (baij) {
743:     if (dim == 2) {
744:       DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
745:     } else if (dim == 3) {
746:       DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
747:     } 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);
748:   } else if (sbaij) {
749:     if (dim == 2) {
750:       DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
751:     } else if (dim == 3) {
752:       DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
753:     } 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);
754:   } else {
755:     ISLocalToGlobalMapping ltog;
756:     MatSetBlockSize(A,dof);
757:     MatSetUp(A);
758:     DMGetLocalToGlobalMapping(da,&ltog);
759:     MatSetLocalToGlobalMapping(A,ltog,ltog);
760:   }
761:   DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
762:   MatSetStencil(A,dim,dims,starts,dof);
763:   MatSetDM(A,da);
764:   MPI_Comm_size(comm,&size);
765:   if (size > 1) {
766:     /* change viewer to display matrix in natural ordering */
767:     MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
768:     MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
769:   }
770:   MatSetFromOptions(A);
771:   *J = A;
772:   return(0);
773: }

775: /* ---------------------------------------------------------------------------------*/
776: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
777: {
778:   PetscErrorCode         ierr;
779:   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;
780:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
781:   MPI_Comm               comm;
782:   PetscScalar            *values;
783:   DMBoundaryType         bx,by;
784:   ISLocalToGlobalMapping ltog;
785:   DMDAStencilType        st;
786:   PetscBool              removedups = PETSC_FALSE;

789:   /*
790:          nc - number of components per grid point
791:          col - number of colors needed in one direction for single component problem

793:   */
794:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
795:   col  = 2*s + 1;
796:   /*
797:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
798:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
799:   */
800:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
801:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
802:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
803:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
804:   PetscObjectGetComm((PetscObject)da,&comm);

806:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
807:   DMGetLocalToGlobalMapping(da,&ltog);

809:   MatSetBlockSize(J,nc);
810:   /* determine the matrix preallocation information */
811:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
812:   for (i=xs; i<xs+nx; i++) {

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

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

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

823:       cnt = 0;
824:       for (k=0; k<nc; k++) {
825:         for (l=lstart; l<lend+1; l++) {
826:           for (p=pstart; p<pend+1; p++) {
827:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
828:               cols[cnt++] = k + nc*(slot + gnx*l + p);
829:             }
830:           }
831:         }
832:         rows[k] = k + nc*(slot);
833:       }
834:       if (removedups) {
835:         MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
836:       } else {
837:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
838:       }
839:     }
840:   }
841:   MatSetBlockSize(J,nc);
842:   MatSeqAIJSetPreallocation(J,0,dnz);
843:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
844:   MatPreallocateFinalize(dnz,onz);

846:   MatSetLocalToGlobalMapping(J,ltog,ltog);

848:   /*
849:     For each node in the grid: we get the neighbors in the local (on processor ordering
850:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
851:     PETSc ordering.
852:   */
853:   if (!da->prealloc_only) {
854:     PetscCalloc1(col*col*nc*nc,&values);
855:     for (i=xs; i<xs+nx; i++) {

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

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

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

866:         cnt = 0;
867:         for (k=0; k<nc; k++) {
868:           for (l=lstart; l<lend+1; l++) {
869:             for (p=pstart; p<pend+1; p++) {
870:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
871:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
872:               }
873:             }
874:           }
875:           rows[k] = k + nc*(slot);
876:         }
877:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
878:       }
879:     }
880:     PetscFree(values);
881:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
882:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
883:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
884:   }
885:   PetscFree2(rows,cols);
886:   return(0);
887: }

889: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
890: {
891:   PetscErrorCode         ierr;
892:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
893:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
894:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
895:   DM_DA                  *dd = (DM_DA*)da->data;
896:   PetscInt               ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
897:   MPI_Comm               comm;
898:   PetscScalar            *values;
899:   DMBoundaryType         bx,by;
900:   ISLocalToGlobalMapping ltog;
901:   DMDAStencilType        st;
902:   PetscBool              removedups = PETSC_FALSE;

905:   /*
906:          nc - number of components per grid point
907:          col - number of colors needed in one direction for single component problem

909:   */
910:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
911:   col  = 2*s + 1;
912:   /*
913:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
914:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
915:   */
916:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
917:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
918:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
919:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
920:   PetscObjectGetComm((PetscObject)da,&comm);

922:   PetscMalloc1(col*col*nc,&cols);
923:   DMGetLocalToGlobalMapping(da,&ltog);

925:   MatSetBlockSize(J,nc);
926:   /* determine the matrix preallocation information */
927:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
928:   for (i=xs; i<xs+nx; i++) {

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

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

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

939:       for (k=0; k<nc; k++) {
940:         cnt = 0;
941:         for (l=lstart; l<lend+1; l++) {
942:           for (p=pstart; p<pend+1; p++) {
943:             if (l || p) {
944:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
945:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
946:               }
947:             } else {
948:               if (dfill) {
949:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
950:               } else {
951:                 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
952:               }
953:             }
954:           }
955:         }
956:         row    = k + nc*(slot);
957:         maxcnt = PetscMax(maxcnt,cnt);
958:         if (removedups) {
959:           MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
960:         } else {
961:           MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
962:         }
963:       }
964:     }
965:   }
966:   MatSeqAIJSetPreallocation(J,0,dnz);
967:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
968:   MatPreallocateFinalize(dnz,onz);
969:   MatSetLocalToGlobalMapping(J,ltog,ltog);

971:   /*
972:     For each node in the grid: we get the neighbors in the local (on processor ordering
973:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
974:     PETSc ordering.
975:   */
976:   if (!da->prealloc_only) {
977:     PetscCalloc1(maxcnt,&values);
978:     for (i=xs; i<xs+nx; i++) {

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

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

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

989:         for (k=0; k<nc; k++) {
990:           cnt = 0;
991:           for (l=lstart; l<lend+1; l++) {
992:             for (p=pstart; p<pend+1; p++) {
993:               if (l || p) {
994:                 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
995:                   for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
996:                 }
997:               } else {
998:                 if (dfill) {
999:                   for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1000:                 } else {
1001:                   for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1002:                 }
1003:               }
1004:             }
1005:           }
1006:           row  = k + nc*(slot);
1007:           MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1008:         }
1009:       }
1010:     }
1011:     PetscFree(values);
1012:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1013:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1014:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1015:   }
1016:   PetscFree(cols);
1017:   return(0);
1018: }

1020: /* ---------------------------------------------------------------------------------*/

1022: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1023: {
1024:   PetscErrorCode         ierr;
1025:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1026:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1027:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1028:   MPI_Comm               comm;
1029:   PetscScalar            *values;
1030:   DMBoundaryType         bx,by,bz;
1031:   ISLocalToGlobalMapping ltog;
1032:   DMDAStencilType        st;
1033:   PetscBool              removedups = PETSC_FALSE;

1036:   /*
1037:          nc - number of components per grid point
1038:          col - number of colors needed in one direction for single component problem

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

1044:   /*
1045:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1046:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1047:   */
1048:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1049:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1050:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

1056:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1057:   DMGetLocalToGlobalMapping(da,&ltog);

1059:   MatSetBlockSize(J,nc);
1060:   /* determine the matrix preallocation information */
1061:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1062:   for (i=xs; i<xs+nx; i++) {
1063:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1064:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1065:     for (j=ys; j<ys+ny; j++) {
1066:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1067:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1068:       for (k=zs; k<zs+nz; k++) {
1069:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1070:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1074:         cnt = 0;
1075:         for (l=0; l<nc; l++) {
1076:           for (ii=istart; ii<iend+1; ii++) {
1077:             for (jj=jstart; jj<jend+1; jj++) {
1078:               for (kk=kstart; kk<kend+1; kk++) {
1079:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1080:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1081:                 }
1082:               }
1083:             }
1084:           }
1085:           rows[l] = l + nc*(slot);
1086:         }
1087:         if (removedups) {
1088:           MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1089:         } else {
1090:           MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1091:         }
1092:       }
1093:     }
1094:   }
1095:   MatSetBlockSize(J,nc);
1096:   MatSeqAIJSetPreallocation(J,0,dnz);
1097:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1098:   MatPreallocateFinalize(dnz,onz);
1099:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1101:   /*
1102:     For each node in the grid: we get the neighbors in the local (on processor ordering
1103:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1104:     PETSc ordering.
1105:   */
1106:   if (!da->prealloc_only) {
1107:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1108:     for (i=xs; i<xs+nx; i++) {
1109:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1110:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1111:       for (j=ys; j<ys+ny; j++) {
1112:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1113:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1114:         for (k=zs; k<zs+nz; k++) {
1115:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1116:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1120:           cnt = 0;
1121:           for (l=0; l<nc; l++) {
1122:             for (ii=istart; ii<iend+1; ii++) {
1123:               for (jj=jstart; jj<jend+1; jj++) {
1124:                 for (kk=kstart; kk<kend+1; kk++) {
1125:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1126:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1127:                   }
1128:                 }
1129:               }
1130:             }
1131:             rows[l] = l + nc*(slot);
1132:           }
1133:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1134:         }
1135:       }
1136:     }
1137:     PetscFree(values);
1138:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1139:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1140:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1141:   }
1142:   PetscFree2(rows,cols);
1143:   return(0);
1144: }

1146: /* ---------------------------------------------------------------------------------*/

1148: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1149: {
1150:   PetscErrorCode         ierr;
1151:   DM_DA                  *dd = (DM_DA*)da->data;
1152:   PetscInt               xs,nx,i,j,gxs,gnx,row,k,l;
1153:   PetscInt               m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1154:   PetscInt               *ofill = dd->ofill,*dfill = dd->dfill;
1155:   PetscScalar            *values;
1156:   DMBoundaryType         bx;
1157:   ISLocalToGlobalMapping ltog;
1158:   PetscMPIInt            rank,size;

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

1165:   /*
1166:          nc - number of components per grid point

1168:   */
1169:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1170:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1171:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1173:   MatSetBlockSize(J,nc);
1174:   PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);

1176:   /*
1177:         note should be smaller for first and last process with no periodic
1178:         does not handle dfill
1179:   */
1180:   cnt = 0;
1181:   /* coupling with process to the left */
1182:   for (i=0; i<s; i++) {
1183:     for (j=0; j<nc; j++) {
1184:       ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1185:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1186:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1187:       cnt++;
1188:     }
1189:   }
1190:   for (i=s; i<nx-s; i++) {
1191:     for (j=0; j<nc; j++) {
1192:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1193:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1194:       cnt++;
1195:     }
1196:   }
1197:   /* coupling with process to the right */
1198:   for (i=nx-s; i<nx; i++) {
1199:     for (j=0; j<nc; j++) {
1200:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1201:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1202:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1203:       cnt++;
1204:     }
1205:   }

1207:   MatSeqAIJSetPreallocation(J,0,cols);
1208:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1209:   PetscFree2(cols,ocols);

1211:   DMGetLocalToGlobalMapping(da,&ltog);
1212:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1214:   /*
1215:     For each node in the grid: we get the neighbors in the local (on processor ordering
1216:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1217:     PETSc ordering.
1218:   */
1219:   if (!da->prealloc_only) {
1220:     PetscCalloc2(maxcnt,&values,maxcnt,&cols);

1222:     row = xs*nc;
1223:     /* coupling with process to the left */
1224:     for (i=xs; i<xs+s; i++) {
1225:       for (j=0; j<nc; j++) {
1226:         cnt = 0;
1227:         if (rank) {
1228:           for (l=0; l<s; l++) {
1229:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1230:           }
1231:         }
1232:         if (dfill) {
1233:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1234:             cols[cnt++] = i*nc + dfill[k];
1235:           }
1236:         } else {
1237:           for (k=0; k<nc; k++) {
1238:             cols[cnt++] = i*nc + k;
1239:           }
1240:         }
1241:         for (l=0; l<s; l++) {
1242:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1243:         }
1244:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1245:         row++;
1246:       }
1247:     }
1248:     for (i=xs+s; i<xs+nx-s; i++) {
1249:       for (j=0; j<nc; j++) {
1250:         cnt = 0;
1251:         for (l=0; l<s; l++) {
1252:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1253:         }
1254:         if (dfill) {
1255:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1256:             cols[cnt++] = i*nc + dfill[k];
1257:           }
1258:         } else {
1259:           for (k=0; k<nc; k++) {
1260:             cols[cnt++] = i*nc + k;
1261:           }
1262:         }
1263:         for (l=0; l<s; l++) {
1264:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1265:         }
1266:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1267:         row++;
1268:       }
1269:     }
1270:     /* coupling with process to the right */
1271:     for (i=xs+nx-s; i<xs+nx; i++) {
1272:       for (j=0; j<nc; j++) {
1273:         cnt = 0;
1274:         for (l=0; l<s; l++) {
1275:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1276:         }
1277:         if (dfill) {
1278:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1279:             cols[cnt++] = i*nc + dfill[k];
1280:           }
1281:         } else {
1282:           for (k=0; k<nc; k++) {
1283:             cols[cnt++] = i*nc + k;
1284:           }
1285:         }
1286:         if (rank < size-1) {
1287:           for (l=0; l<s; l++) {
1288:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1289:           }
1290:         }
1291:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1292:         row++;
1293:       }
1294:     }
1295:     PetscFree2(values,cols);
1296:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1297:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1298:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1299:   }
1300:   return(0);
1301: }

1303: /* ---------------------------------------------------------------------------------*/

1305: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1306: {
1307:   PetscErrorCode         ierr;
1308:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1309:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1310:   PetscInt               istart,iend;
1311:   PetscScalar            *values;
1312:   DMBoundaryType         bx;
1313:   ISLocalToGlobalMapping ltog;

1316:   /*
1317:          nc - number of components per grid point
1318:          col - number of colors needed in one direction for single component problem

1320:   */
1321:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1322:   col  = 2*s + 1;

1324:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1325:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1327:   MatSetBlockSize(J,nc);
1328:   MatSeqAIJSetPreallocation(J,col*nc,0);
1329:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1331:   DMGetLocalToGlobalMapping(da,&ltog);
1332:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1334:   /*
1335:     For each node in the grid: we get the neighbors in the local (on processor ordering
1336:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1337:     PETSc ordering.
1338:   */
1339:   if (!da->prealloc_only) {
1340:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1341:     PetscCalloc1(col*nc*nc,&values);
1342:     for (i=xs; i<xs+nx; i++) {
1343:       istart = PetscMax(-s,gxs - i);
1344:       iend   = PetscMin(s,gxs + gnx - i - 1);
1345:       slot   = i - gxs;

1347:       cnt = 0;
1348:       for (l=0; l<nc; l++) {
1349:         for (i1=istart; i1<iend+1; i1++) {
1350:           cols[cnt++] = l + nc*(slot + i1);
1351:         }
1352:         rows[l] = l + nc*(slot);
1353:       }
1354:       MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1355:     }
1356:     PetscFree(values);
1357:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1358:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1359:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1360:     PetscFree2(rows,cols);
1361:   }
1362:   return(0);
1363: }

1365: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1366: {
1367:   PetscErrorCode         ierr;
1368:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1369:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1370:   PetscInt               istart,iend,jstart,jend,ii,jj;
1371:   MPI_Comm               comm;
1372:   PetscScalar            *values;
1373:   DMBoundaryType         bx,by;
1374:   DMDAStencilType        st;
1375:   ISLocalToGlobalMapping ltog;

1378:   /*
1379:      nc - number of components per grid point
1380:      col - number of colors needed in one direction for single component problem
1381:   */
1382:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1383:   col  = 2*s + 1;

1385:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1386:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1387:   PetscObjectGetComm((PetscObject)da,&comm);

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

1391:   DMGetLocalToGlobalMapping(da,&ltog);

1393:   /* determine the matrix preallocation information */
1394:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1395:   for (i=xs; i<xs+nx; i++) {
1396:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1397:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1398:     for (j=ys; j<ys+ny; j++) {
1399:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1400:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1401:       slot   = i - gxs + gnx*(j - gys);

1403:       /* Find block columns in block row */
1404:       cnt = 0;
1405:       for (ii=istart; ii<iend+1; ii++) {
1406:         for (jj=jstart; jj<jend+1; jj++) {
1407:           if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1408:             cols[cnt++] = slot + ii + gnx*jj;
1409:           }
1410:         }
1411:       }
1412:       MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1413:     }
1414:   }
1415:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1416:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1417:   MatPreallocateFinalize(dnz,onz);

1419:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1421:   /*
1422:     For each node in the grid: we get the neighbors in the local (on processor ordering
1423:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1424:     PETSc ordering.
1425:   */
1426:   if (!da->prealloc_only) {
1427:     PetscCalloc1(col*col*nc*nc,&values);
1428:     for (i=xs; i<xs+nx; i++) {
1429:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1430:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1431:       for (j=ys; j<ys+ny; j++) {
1432:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1433:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1434:         slot = i - gxs + gnx*(j - gys);
1435:         cnt  = 0;
1436:         for (ii=istart; ii<iend+1; ii++) {
1437:           for (jj=jstart; jj<jend+1; jj++) {
1438:             if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1439:               cols[cnt++] = slot + ii + gnx*jj;
1440:             }
1441:           }
1442:         }
1443:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1444:       }
1445:     }
1446:     PetscFree(values);
1447:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1448:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1449:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1450:   }
1451:   PetscFree(cols);
1452:   return(0);
1453: }

1455: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1456: {
1457:   PetscErrorCode         ierr;
1458:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1459:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1460:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1461:   MPI_Comm               comm;
1462:   PetscScalar            *values;
1463:   DMBoundaryType         bx,by,bz;
1464:   DMDAStencilType        st;
1465:   ISLocalToGlobalMapping ltog;

1468:   /*
1469:          nc - number of components per grid point
1470:          col - number of colors needed in one direction for single component problem

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

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

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

1482:   DMGetLocalToGlobalMapping(da,&ltog);

1484:   /* determine the matrix preallocation information */
1485:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1486:   for (i=xs; i<xs+nx; i++) {
1487:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1488:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1489:     for (j=ys; j<ys+ny; j++) {
1490:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1491:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1492:       for (k=zs; k<zs+nz; k++) {
1493:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1494:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1498:         /* Find block columns in block row */
1499:         cnt = 0;
1500:         for (ii=istart; ii<iend+1; ii++) {
1501:           for (jj=jstart; jj<jend+1; jj++) {
1502:             for (kk=kstart; kk<kend+1; kk++) {
1503:               if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1504:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1505:               }
1506:             }
1507:           }
1508:         }
1509:         MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1510:       }
1511:     }
1512:   }
1513:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1514:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1515:   MatPreallocateFinalize(dnz,onz);

1517:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

1561: /*
1562:   This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1563:   identify in the local ordering with periodic domain.
1564: */
1565: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1566: {
1568:   PetscInt       i,n;

1571:   ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1572:   ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1573:   for (i=0,n=0; i<*cnt; i++) {
1574:     if (col[i] >= *row) col[n++] = col[i];
1575:   }
1576:   *cnt = n;
1577:   return(0);
1578: }

1580: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1581: {
1582:   PetscErrorCode         ierr;
1583:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1584:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1585:   PetscInt               istart,iend,jstart,jend,ii,jj;
1586:   MPI_Comm               comm;
1587:   PetscScalar            *values;
1588:   DMBoundaryType         bx,by;
1589:   DMDAStencilType        st;
1590:   ISLocalToGlobalMapping ltog;

1593:   /*
1594:      nc - number of components per grid point
1595:      col - number of colors needed in one direction for single component problem
1596:   */
1597:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1598:   col  = 2*s + 1;

1600:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1601:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1602:   PetscObjectGetComm((PetscObject)da,&comm);

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

1606:   DMGetLocalToGlobalMapping(da,&ltog);

1608:   /* determine the matrix preallocation information */
1609:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1610:   for (i=xs; i<xs+nx; i++) {
1611:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1612:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1613:     for (j=ys; j<ys+ny; j++) {
1614:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1615:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1616:       slot   = i - gxs + gnx*(j - gys);

1618:       /* Find block columns in block row */
1619:       cnt = 0;
1620:       for (ii=istart; ii<iend+1; ii++) {
1621:         for (jj=jstart; jj<jend+1; jj++) {
1622:           if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1623:             cols[cnt++] = slot + ii + gnx*jj;
1624:           }
1625:         }
1626:       }
1627:       L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1628:       MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1629:     }
1630:   }
1631:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1632:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1633:   MatPreallocateFinalize(dnz,onz);

1635:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1637:   /*
1638:     For each node in the grid: we get the neighbors in the local (on processor ordering
1639:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1640:     PETSc ordering.
1641:   */
1642:   if (!da->prealloc_only) {
1643:     PetscCalloc1(col*col*nc*nc,&values);
1644:     for (i=xs; i<xs+nx; i++) {
1645:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1646:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1647:       for (j=ys; j<ys+ny; j++) {
1648:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1649:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1650:         slot   = i - gxs + gnx*(j - gys);

1652:         /* Find block columns in block row */
1653:         cnt = 0;
1654:         for (ii=istart; ii<iend+1; ii++) {
1655:           for (jj=jstart; jj<jend+1; jj++) {
1656:             if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1657:               cols[cnt++] = slot + ii + gnx*jj;
1658:             }
1659:           }
1660:         }
1661:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1662:         MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1663:       }
1664:     }
1665:     PetscFree(values);
1666:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1667:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1668:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1669:   }
1670:   PetscFree(cols);
1671:   return(0);
1672: }

1674: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1675: {
1676:   PetscErrorCode         ierr;
1677:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1678:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1679:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1680:   MPI_Comm               comm;
1681:   PetscScalar            *values;
1682:   DMBoundaryType         bx,by,bz;
1683:   DMDAStencilType        st;
1684:   ISLocalToGlobalMapping ltog;

1687:   /*
1688:      nc - number of components per grid point
1689:      col - number of colors needed in one direction for single component problem
1690:   */
1691:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1692:   col  = 2*s + 1;

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

1698:   /* create the matrix */
1699:   PetscMalloc1(col*col*col,&cols);

1701:   DMGetLocalToGlobalMapping(da,&ltog);

1703:   /* determine the matrix preallocation information */
1704:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1705:   for (i=xs; i<xs+nx; i++) {
1706:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1707:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1708:     for (j=ys; j<ys+ny; j++) {
1709:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1710:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1711:       for (k=zs; k<zs+nz; k++) {
1712:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1713:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1717:         /* Find block columns in block row */
1718:         cnt = 0;
1719:         for (ii=istart; ii<iend+1; ii++) {
1720:           for (jj=jstart; jj<jend+1; jj++) {
1721:             for (kk=kstart; kk<kend+1; kk++) {
1722:               if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1723:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1724:               }
1725:             }
1726:           }
1727:         }
1728:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1729:         MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1730:       }
1731:     }
1732:   }
1733:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1734:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1735:   MatPreallocateFinalize(dnz,onz);

1737:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1739:   /*
1740:     For each node in the grid: we get the neighbors in the local (on processor ordering
1741:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1742:     PETSc ordering.
1743:   */
1744:   if (!da->prealloc_only) {
1745:     PetscCalloc1(col*col*col*nc*nc,&values);
1746:     for (i=xs; i<xs+nx; i++) {
1747:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1748:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1749:       for (j=ys; j<ys+ny; j++) {
1750:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1751:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1752:         for (k=zs; k<zs+nz; k++) {
1753:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1754:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1758:           cnt = 0;
1759:           for (ii=istart; ii<iend+1; ii++) {
1760:             for (jj=jstart; jj<jend+1; jj++) {
1761:               for (kk=kstart; kk<kend+1; kk++) {
1762:                 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1763:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1764:                 }
1765:               }
1766:             }
1767:           }
1768:           L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1769:           MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1770:         }
1771:       }
1772:     }
1773:     PetscFree(values);
1774:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1775:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1776:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1777:   }
1778:   PetscFree(cols);
1779:   return(0);
1780: }

1782: /* ---------------------------------------------------------------------------------*/

1784: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
1785: {
1786:   PetscErrorCode         ierr;
1787:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1788:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
1789:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1790:   DM_DA                  *dd = (DM_DA*)da->data;
1791:   PetscInt               ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
1792:   MPI_Comm               comm;
1793:   PetscScalar            *values;
1794:   DMBoundaryType         bx,by,bz;
1795:   ISLocalToGlobalMapping ltog;
1796:   DMDAStencilType        st;
1797:   PetscBool              removedups = PETSC_FALSE;

1800:   /*
1801:          nc - number of components per grid point
1802:          col - number of colors needed in one direction for single component problem

1804:   */
1805:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1806:   col  = 2*s + 1;
1807:   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\
1808:                  by 2*stencil_width + 1\n");
1809:   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\
1810:                  by 2*stencil_width + 1\n");
1811:   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\
1812:                  by 2*stencil_width + 1\n");

1814:   /*
1815:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1816:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1817:   */
1818:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1819:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1820:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

1826:   PetscMalloc1(col*col*col*nc,&cols);
1827:   DMGetLocalToGlobalMapping(da,&ltog);

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

1832:   MatSetBlockSize(J,nc);
1833:   for (i=xs; i<xs+nx; i++) {
1834:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1835:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1836:     for (j=ys; j<ys+ny; j++) {
1837:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1838:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1839:       for (k=zs; k<zs+nz; k++) {
1840:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1841:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1845:         for (l=0; l<nc; l++) {
1846:           cnt = 0;
1847:           for (ii=istart; ii<iend+1; ii++) {
1848:             for (jj=jstart; jj<jend+1; jj++) {
1849:               for (kk=kstart; kk<kend+1; kk++) {
1850:                 if (ii || jj || kk) {
1851:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1852:                     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);
1853:                   }
1854:                 } else {
1855:                   if (dfill) {
1856:                     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);
1857:                   } else {
1858:                     for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1859:                   }
1860:                 }
1861:               }
1862:             }
1863:           }
1864:           row  = l + nc*(slot);
1865:           maxcnt = PetscMax(maxcnt,cnt);
1866:           if (removedups) {
1867:             MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1868:           } else {
1869:             MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1870:           }
1871:         }
1872:       }
1873:     }
1874:   }
1875:   MatSeqAIJSetPreallocation(J,0,dnz);
1876:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1877:   MatPreallocateFinalize(dnz,onz);
1878:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1880:   /*
1881:     For each node in the grid: we get the neighbors in the local (on processor ordering
1882:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1883:     PETSc ordering.
1884:   */
1885:   if (!da->prealloc_only) {
1886:     PetscCalloc1(maxcnt,&values);
1887:     for (i=xs; i<xs+nx; i++) {
1888:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1889:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1890:       for (j=ys; j<ys+ny; j++) {
1891:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1892:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1893:         for (k=zs; k<zs+nz; k++) {
1894:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1895:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1899:           for (l=0; l<nc; l++) {
1900:             cnt = 0;
1901:             for (ii=istart; ii<iend+1; ii++) {
1902:               for (jj=jstart; jj<jend+1; jj++) {
1903:                 for (kk=kstart; kk<kend+1; kk++) {
1904:                   if (ii || jj || kk) {
1905:                     if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1906:                       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);
1907:                     }
1908:                   } else {
1909:                     if (dfill) {
1910:                       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);
1911:                     } else {
1912:                       for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1913:                     }
1914:                   }
1915:                 }
1916:               }
1917:             }
1918:             row  = l + nc*(slot);
1919:             MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1920:           }
1921:         }
1922:       }
1923:     }
1924:     PetscFree(values);
1925:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1926:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1927:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1928:   }
1929:   PetscFree(cols);
1930:   return(0);
1931: }