Actual source code: fdda.c

petsc-3.14.6 2021-03-30
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:   PetscArraycpy(*rfill,dfillsparse,nz+w+1);
 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 da

 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 da

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,NULL,NULL,NULL,&m,&n,&p,&nc,NULL,&bx,&by,&bz,NULL);

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:   PetscStrbeginswith(da->mattype,MATBAIJ,&isBAIJ);
252:   if (!isBAIJ) {PetscStrbeginswith(da->mattype,MATMPIBAIJ,&isBAIJ);}
253:   if (!isBAIJ) {PetscStrbeginswith(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
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,NULL,&M,&N,NULL,&nc,&s,&bx,&by,NULL,&st);
304:   col  = 2*s + 1;
305:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
306:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
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 {
313:     if (ctype == IS_COLORING_GLOBAL) {
314:       if (!dd->localcoloring) {
315:         PetscMalloc1(nc*nx*ny,&colors);
316:         ii   = 0;
317:         for (j=ys; j<ys+ny; j++) {
318:           for (i=xs; i<xs+nx; i++) {
319:             for (k=0; k<nc; k++) {
320:               colors[ii++] = k + nc*((i % col) + col*(j % col));
321:             }
322:           }
323:         }
324:         ncolors = nc + nc*(col-1 + col*(col-1));
325:         ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
326:       }
327:       *coloring = dd->localcoloring;
328:     } else if (ctype == IS_COLORING_LOCAL) {
329:       if (!dd->ghostedcoloring) {
330:         PetscMalloc1(nc*gnx*gny,&colors);
331:         ii   = 0;
332:         for (j=gys; j<gys+gny; j++) {
333:           for (i=gxs; i<gxs+gnx; i++) {
334:             for (k=0; k<nc; k++) {
335:               /* the complicated stuff is to handle periodic boundaries */
336:               colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
337:             }
338:           }
339:         }
340:         ncolors = nc + nc*(col - 1 + col*(col-1));
341:         ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
342:         /* PetscIntView(ncolors,(PetscInt*)colors,0); */

344:         ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
345:       }
346:       *coloring = dd->ghostedcoloring;
347:     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
348:   }
349:   ISColoringReference(*coloring);
350:   return(0);
351: }

353: /* ---------------------------------------------------------------------------------*/

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

367:   /*
368:          nc - number of components per grid point
369:          col - number of colors needed in one direction for single component problem

371:   */
372:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
373:   col  = 2*s + 1;
374:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
375:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
376:   PetscObjectGetComm((PetscObject)da,&comm);

378:   /* create the coloring */
379:   if (ctype == IS_COLORING_GLOBAL) {
380:     if (!dd->localcoloring) {
381:       PetscMalloc1(nc*nx*ny*nz,&colors);
382:       ii   = 0;
383:       for (k=zs; k<zs+nz; k++) {
384:         for (j=ys; j<ys+ny; j++) {
385:           for (i=xs; i<xs+nx; i++) {
386:             for (l=0; l<nc; l++) {
387:               colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
388:             }
389:           }
390:         }
391:       }
392:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
393:       ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
394:     }
395:     *coloring = dd->localcoloring;
396:   } else if (ctype == IS_COLORING_LOCAL) {
397:     if (!dd->ghostedcoloring) {
398:       PetscMalloc1(nc*gnx*gny*gnz,&colors);
399:       ii   = 0;
400:       for (k=gzs; k<gzs+gnz; k++) {
401:         for (j=gys; j<gys+gny; j++) {
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[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
406:             }
407:           }
408:         }
409:       }
410:       ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
411:       ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
412:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
413:     }
414:     *coloring = dd->ghostedcoloring;
415:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
416:   ISColoringReference(*coloring);
417:   return(0);
418: }

420: /* ---------------------------------------------------------------------------------*/

422: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
423: {
424:   PetscErrorCode  ierr;
425:   PetscInt        xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
426:   PetscInt        ncolors;
427:   MPI_Comm        comm;
428:   DMBoundaryType  bx;
429:   ISColoringValue *colors;
430:   DM_DA           *dd = (DM_DA*)da->data;

433:   /*
434:          nc - number of components per grid point
435:          col - number of colors needed in one direction for single component problem

437:   */
438:   DMDAGetInfo(da,&dim,&m,NULL,NULL,&M,NULL,NULL,&nc,&s,&bx,NULL,NULL,NULL);
439:   col  = 2*s + 1;
440:   DMDAGetCorners(da,&xs,NULL,NULL,&nx,NULL,NULL);
441:   DMDAGetGhostCorners(da,&gxs,NULL,NULL,&gnx,NULL,NULL);
442:   PetscObjectGetComm((PetscObject)da,&comm);

444:   /* create the coloring */
445:   if (ctype == IS_COLORING_GLOBAL) {
446:     if (!dd->localcoloring) {
447:       PetscMalloc1(nc*nx,&colors);
448:       if (dd->ofillcols) {
449:         PetscInt tc = 0;
450:         for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
451:         i1 = 0;
452:         for (i=xs; i<xs+nx; i++) {
453:           for (l=0; l<nc; l++) {
454:             if (dd->ofillcols[l] && (i % col)) {
455:               colors[i1++] =  nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
456:             } else {
457:               colors[i1++] = l;
458:             }
459:           }
460:         }
461:         ncolors = nc + 2*s*tc;
462:       } else {
463:         i1 = 0;
464:         for (i=xs; i<xs+nx; i++) {
465:           for (l=0; l<nc; l++) {
466:             colors[i1++] = l + nc*(i % col);
467:           }
468:         }
469:         ncolors = nc + nc*(col-1);
470:       }
471:       ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
472:     }
473:     *coloring = dd->localcoloring;
474:   } else if (ctype == IS_COLORING_LOCAL) {
475:     if (!dd->ghostedcoloring) {
476:       PetscMalloc1(nc*gnx,&colors);
477:       i1   = 0;
478:       for (i=gxs; i<gxs+gnx; i++) {
479:         for (l=0; l<nc; l++) {
480:           /* the complicated stuff is to handle periodic boundaries */
481:           colors[i1++] = l + nc*(SetInRange(i,m) % col);
482:         }
483:       }
484:       ncolors = nc + nc*(col-1);
485:       ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
486:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
487:     }
488:     *coloring = dd->ghostedcoloring;
489:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
490:   ISColoringReference(*coloring);
491:   return(0);
492: }

494: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
495: {
496:   PetscErrorCode  ierr;
497:   PetscInt        xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
498:   PetscInt        ncolors;
499:   MPI_Comm        comm;
500:   DMBoundaryType  bx,by;
501:   ISColoringValue *colors;
502:   DM_DA           *dd = (DM_DA*)da->data;

505:   /*
506:          nc - number of components per grid point
507:          col - number of colors needed in one direction for single component problem

509:   */
510:   DMDAGetInfo(da,&dim,&m,&n,NULL,NULL,NULL,NULL,&nc,&s,&bx,&by,NULL,NULL);
511:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
512:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
513:   PetscObjectGetComm((PetscObject)da,&comm);
514:   /* create the coloring */
515:   if (ctype == IS_COLORING_GLOBAL) {
516:     if (!dd->localcoloring) {
517:       PetscMalloc1(nc*nx*ny,&colors);
518:       ii   = 0;
519:       for (j=ys; j<ys+ny; j++) {
520:         for (i=xs; i<xs+nx; i++) {
521:           for (k=0; k<nc; k++) {
522:             colors[ii++] = k + nc*((3*j+i) % 5);
523:           }
524:         }
525:       }
526:       ncolors = 5*nc;
527:       ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
528:     }
529:     *coloring = dd->localcoloring;
530:   } else if (ctype == IS_COLORING_LOCAL) {
531:     if (!dd->ghostedcoloring) {
532:       PetscMalloc1(nc*gnx*gny,&colors);
533:       ii = 0;
534:       for (j=gys; j<gys+gny; j++) {
535:         for (i=gxs; i<gxs+gnx; i++) {
536:           for (k=0; k<nc; k++) {
537:             colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
538:           }
539:         }
540:       }
541:       ncolors = 5*nc;
542:       ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
543:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
544:     }
545:     *coloring = dd->ghostedcoloring;
546:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
547:   return(0);
548: }

550: /* =========================================================================== */
551: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat,PetscBool);
552: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
553: extern PetscErrorCode DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(DM,Mat,PetscBool);
554: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat,PetscBool);
555: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
556: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat,PetscBool);
557: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
558: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
559: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
560: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
561: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
562: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM,Mat);
563: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM,Mat);
564: extern PetscErrorCode DMCreateMatrix_DA_IS(DM,Mat);

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

569:    Logically Collective on mat

571:    Input Parameters:
572: +  mat - the matrix
573: -  da - the da

575:    Level: intermediate

577: @*/
578: PetscErrorCode MatSetupDM(Mat mat,DM da)
579: {

585:   PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
586:   return(0);
587: }

589: PetscErrorCode  MatView_MPI_DA(Mat A,PetscViewer viewer)
590: {
591:   DM                da;
592:   PetscErrorCode    ierr;
593:   const char        *prefix;
594:   Mat               Anatural;
595:   AO                ao;
596:   PetscInt          rstart,rend,*petsc,i;
597:   IS                is;
598:   MPI_Comm          comm;
599:   PetscViewerFormat format;

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

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

610:   DMDAGetAO(da,&ao);
611:   MatGetOwnershipRange(A,&rstart,&rend);
612:   PetscMalloc1(rend-rstart,&petsc);
613:   for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
614:   AOApplicationToPetsc(ao,rend-rstart,petsc);
615:   ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);

617:   /* call viewer on natural ordering */
618:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
619:   ISDestroy(&is);
620:   PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
621:   PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
622:   PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
623:   ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
624:   MatView(Anatural,viewer);
625:   ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
626:   MatDestroy(&Anatural);
627:   return(0);
628: }

630: PetscErrorCode  MatLoad_MPI_DA(Mat A,PetscViewer viewer)
631: {
632:   DM             da;
634:   Mat            Anatural,Aapp;
635:   AO             ao;
636:   PetscInt       rstart,rend,*app,i,m,n,M,N;
637:   IS             is;
638:   MPI_Comm       comm;

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

645:   /* Load the matrix in natural ordering */
646:   MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
647:   MatSetType(Anatural,((PetscObject)A)->type_name);
648:   MatGetSize(A,&M,&N);
649:   MatGetLocalSize(A,&m,&n);
650:   MatSetSizes(Anatural,m,n,M,N);
651:   MatLoad(Anatural,viewer);

653:   /* Map natural ordering to application ordering and create IS */
654:   DMDAGetAO(da,&ao);
655:   MatGetOwnershipRange(Anatural,&rstart,&rend);
656:   PetscMalloc1(rend-rstart,&app);
657:   for (i=rstart; i<rend; i++) app[i-rstart] = i;
658:   AOPetscToApplication(ao,rend-rstart,app);
659:   ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);

661:   /* Do permutation and replace header */
662:   MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
663:   MatHeaderReplace(A,&Aapp);
664:   ISDestroy(&is);
665:   MatDestroy(&Anatural);
666:   return(0);
667: }

669: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
670: {
672:   PetscInt       dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
673:   Mat            A;
674:   MPI_Comm       comm;
675:   MatType        Atype;
676:   void           (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL,(*sell)(void)=NULL,(*is)(void)=NULL;
677:   MatType        mtype;
678:   PetscMPIInt    size;
679:   DM_DA          *dd = (DM_DA*)da->data;

682:   MatInitializePackage();
683:   mtype = da->mattype;

685:   /*
686:                                   m
687:           ------------------------------------------------------
688:          |                                                     |
689:          |                                                     |
690:          |               ----------------------                |
691:          |               |                    |                |
692:       n  |           ny  |                    |                |
693:          |               |                    |                |
694:          |               .---------------------                |
695:          |             (xs,ys)     nx                          |
696:          |            .                                        |
697:          |         (gxs,gys)                                   |
698:          |                                                     |
699:           -----------------------------------------------------
700:   */

702:   /*
703:          nc - number of components per grid point
704:          col - number of colors needed in one direction for single component problem

706:   */
707:   M   = dd->M;
708:   N   = dd->N;
709:   P   = dd->P;
710:   dim = da->dim;
711:   dof = dd->w;
712:   /* DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL); */
713:   DMDAGetCorners(da,NULL,NULL,NULL,&nx,&ny,&nz);
714:   PetscObjectGetComm((PetscObject)da,&comm);
715:   MatCreate(comm,&A);
716:   MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
717:   MatSetType(A,mtype);
718:   MatSetFromOptions(A);
719:   MatSetDM(A,da);
720:   if (da->structure_only) {
721:     MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);
722:   }
723:   MatGetType(A,&Atype);
724:   /*
725:      We do not provide a getmatrix function in the DMDA operations because
726:    the basic DMDA does not know about matrices. We think of DMDA as being more
727:    more low-level than matrices. This is kind of cheating but, cause sometimes
728:    we think of DMDA has higher level than matrices.

730:      We could switch based on Atype (or mtype), but we do not since the
731:    specialized setting routines depend only on the particular preallocation
732:    details of the matrix, not the type itself.
733:   */
734:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
735:   if (!aij) {
736:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
737:   }
738:   if (!aij) {
739:     PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
740:     if (!baij) {
741:       PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
742:     }
743:     if (!baij) {
744:       PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
745:       if (!sbaij) {
746:         PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
747:       }
748:       if (!sbaij) {
749:         PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);
750:         if (!sell) {
751:           PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);
752:         }
753:       }
754:       if (!sell) {
755:         PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
756:       }
757:     }
758:   }
759:   if (aij) {
760:     if (dim == 1) {
761:       if (dd->ofill) {
762:         DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
763:       } else {
764:         DMBoundaryType bx;
765:         PetscMPIInt  size;
766:         DMDAGetInfo(da,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&bx,NULL,NULL,NULL);
767:         MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
768:         if (size == 1 && bx == DM_BOUNDARY_NONE) {
769:           DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(da,A,PETSC_FALSE);
770:         } else {
771:           DMCreateMatrix_DA_1d_MPIAIJ(da,A,PETSC_FALSE);
772:         }
773:       }
774:     } else if (dim == 2) {
775:       if (dd->ofill) {
776:         DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
777:       } else {
778:         DMCreateMatrix_DA_2d_MPIAIJ(da,A,PETSC_FALSE);
779:       }
780:     } else if (dim == 3) {
781:       if (dd->ofill) {
782:         DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
783:       } else {
784:         DMCreateMatrix_DA_3d_MPIAIJ(da,A,PETSC_FALSE);
785:       }
786:     }
787:   } else if (baij) {
788:     if (dim == 2) {
789:       DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
790:     } else if (dim == 3) {
791:       DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
792:     } 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);
793:   } else if (sbaij) {
794:     if (dim == 2) {
795:       DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
796:     } else if (dim == 3) {
797:       DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
798:     } 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);
799:   } else if (sell) {
800:      if (dim == 2) {
801:        DMCreateMatrix_DA_2d_MPISELL(da,A);
802:      } else if (dim == 3) {
803:        DMCreateMatrix_DA_3d_MPISELL(da,A);
804:      } 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);
805:   } else if (is) {
806:     DMCreateMatrix_DA_IS(da,A);
807:   } else {
808:     ISLocalToGlobalMapping ltog;

810:     MatSetBlockSize(A,dof);
811:     MatSetUp(A);
812:     DMGetLocalToGlobalMapping(da,&ltog);
813:     MatSetLocalToGlobalMapping(A,ltog,ltog);
814:   }
815:   DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
816:   MatSetStencil(A,dim,dims,starts,dof);
817:   MatSetDM(A,da);
818:   MPI_Comm_size(comm,&size);
819:   if (size > 1) {
820:     /* change viewer to display matrix in natural ordering */
821:     MatSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
822:     MatSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
823:   }
824:   *J = A;
825:   return(0);
826: }

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

831: PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J)
832: {
833:   DM_DA                  *da = (DM_DA*)dm->data;
834:   Mat                    lJ;
835:   ISLocalToGlobalMapping ltog;
836:   IS                     is_loc_filt, is_glob;
837:   const PetscInt         *e_loc,*idx;
838:   PetscInt               nel,nen,nv,dof,dim,*gidx,nb;
839:   PetscBool              flg;
840:   PetscErrorCode         ierr;

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

849:   MatSetBlockSize(J,dof);

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

856:   /* obtain a consistent local ordering for MATIS */
857:   ISSortRemoveDups(is_loc_filt);
858:   ISBlockGetLocalSize(is_loc_filt,&nb);
859:   DMGetLocalToGlobalMapping(dm,&ltog);
860:   ISLocalToGlobalMappingGetSize(ltog,&nv);
861:   PetscMalloc1(PetscMax(nb,nv/dof),&gidx);
862:   ISBlockGetIndices(is_loc_filt,&idx);
863:   ISLocalToGlobalMappingApplyBlock(ltog,nb,idx,gidx);
864:   ISBlockRestoreIndices(is_loc_filt,&idx);
865:   ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nb,gidx,PETSC_USE_POINTER,&is_glob);
866:   ISLocalToGlobalMappingCreateIS(is_glob,&ltog);
867:   ISDestroy(&is_glob);
868:   MatSetLocalToGlobalMapping(J,ltog,ltog);
869:   ISLocalToGlobalMappingDestroy(&ltog);

871:   /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */
872:   MatISGetLocalMat(J,&lJ);
873:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
874:   ISDestroy(&is_loc_filt);
875:   ISCreateStride(PetscObjectComm((PetscObject)lJ),nv/dof,0,1,&is_glob);
876:   ISGetIndices(is_glob,&idx);
877:   ISGlobalToLocalMappingApplyBlock(ltog,IS_GTOLM_MASK,nv/dof,idx,&nb,gidx);
878:   ISRestoreIndices(is_glob,&idx);
879:   ISDestroy(&is_glob);
880:   ISLocalToGlobalMappingDestroy(&ltog);
881:   ISCreateBlock(PETSC_COMM_SELF,dof,nb,gidx,PETSC_USE_POINTER,&is_loc_filt);
882:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
883:   ISDestroy(&is_loc_filt);
884:   MatSetLocalToGlobalMapping(lJ,ltog,ltog);
885:   ISLocalToGlobalMappingDestroy(&ltog);
886:   PetscFree(gidx);

888:   /* Preallocation (not exact): we reuse the preallocation routines of the assembled version  */
889:   flg = dm->prealloc_only;
890:   dm->prealloc_only = PETSC_TRUE;
891:   switch (dim) {
892:   case 1:
893:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
894:     DMCreateMatrix_DA_1d_MPIAIJ(dm,J,PETSC_TRUE);
895:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
896:     break;
897:   case 2:
898:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
899:     DMCreateMatrix_DA_2d_MPIAIJ(dm,J,PETSC_TRUE);
900:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
901:     break;
902:   case 3:
903:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
904:     DMCreateMatrix_DA_3d_MPIAIJ(dm,J,PETSC_TRUE);
905:     PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
906:     break;
907:   default:
908:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dimension %d",dim);
909:   }
910:   dm->prealloc_only = flg;
911:   return(0);
912: }

914: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J)
915: {
916:   PetscErrorCode         ierr;
917:   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;
918:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
919:   MPI_Comm               comm;
920:   PetscScalar            *values;
921:   DMBoundaryType         bx,by;
922:   ISLocalToGlobalMapping ltog;
923:   DMDAStencilType        st;

926:   /*
927:          nc - number of components per grid point
928:          col - number of colors needed in one direction for single component problem

930:   */
931:   DMDAGetInfo(da,&dim,&m,&n,NULL,NULL,NULL,NULL,&nc,&s,&bx,&by,NULL,&st);
932:   col  = 2*s + 1;
933:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
934:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
935:   PetscObjectGetComm((PetscObject)da,&comm);

937:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
938:   DMGetLocalToGlobalMapping(da,&ltog);

940:   MatSetBlockSize(J,nc);
941:   /* determine the matrix preallocation information */
942:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
943:   for (i=xs; i<xs+nx; i++) {

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

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

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

954:       cnt = 0;
955:       for (k=0; k<nc; k++) {
956:         for (l=lstart; l<lend+1; l++) {
957:           for (p=pstart; p<pend+1; p++) {
958:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
959:               cols[cnt++] = k + nc*(slot + gnx*l + p);
960:             }
961:           }
962:         }
963:         rows[k] = k + nc*(slot);
964:       }
965:       MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
966:     }
967:   }
968:   MatSetBlockSize(J,nc);
969:   MatSeqSELLSetPreallocation(J,0,dnz);
970:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
971:   MatPreallocateFinalize(dnz,onz);

973:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

993:         cnt = 0;
994:         for (k=0; k<nc; k++) {
995:           for (l=lstart; l<lend+1; l++) {
996:             for (p=pstart; p<pend+1; p++) {
997:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
998:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
999:               }
1000:             }
1001:           }
1002:           rows[k] = k + nc*(slot);
1003:         }
1004:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1005:       }
1006:     }
1007:     PetscFree(values);
1008:     /* do not copy values to GPU since they are all zero and not yet needed there */
1009:     MatBindToCPU(J,PETSC_TRUE);
1010:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1011:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1012:     MatBindToCPU(J,PETSC_FALSE);
1013:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1014:   }
1015:   PetscFree2(rows,cols);
1016:   return(0);
1017: }

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

1032:   /*
1033:          nc - number of components per grid point
1034:          col - number of colors needed in one direction for single component problem

1036:   */
1037:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1038:   col  = 2*s + 1;
1039:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1040:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1041:   PetscObjectGetComm((PetscObject)da,&comm);

1043:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1044:   DMGetLocalToGlobalMapping(da,&ltog);

1046:   MatSetBlockSize(J,nc);
1047:   /* determine the matrix preallocation information */
1048:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1049:   for (i=xs; i<xs+nx; i++) {
1050:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1051:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1052:     for (j=ys; j<ys+ny; j++) {
1053:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1054:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1055:       for (k=zs; k<zs+nz; k++) {
1056:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1057:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1061:         cnt = 0;
1062:         for (l=0; l<nc; l++) {
1063:           for (ii=istart; ii<iend+1; ii++) {
1064:             for (jj=jstart; jj<jend+1; jj++) {
1065:               for (kk=kstart; kk<kend+1; kk++) {
1066:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1067:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1068:                 }
1069:               }
1070:             }
1071:           }
1072:           rows[l] = l + nc*(slot);
1073:         }
1074:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1075:       }
1076:     }
1077:   }
1078:   MatSetBlockSize(J,nc);
1079:   MatSeqSELLSetPreallocation(J,0,dnz);
1080:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1081:   MatPreallocateFinalize(dnz,onz);
1082:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1084:   /*
1085:     For each node in the grid: we get the neighbors in the local (on processor ordering
1086:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1087:     PETSc ordering.
1088:   */
1089:   if (!da->prealloc_only) {
1090:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1091:     for (i=xs; i<xs+nx; i++) {
1092:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1093:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1094:       for (j=ys; j<ys+ny; j++) {
1095:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1096:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1097:         for (k=zs; k<zs+nz; k++) {
1098:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1099:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1103:           cnt = 0;
1104:           for (l=0; l<nc; l++) {
1105:             for (ii=istart; ii<iend+1; ii++) {
1106:               for (jj=jstart; jj<jend+1; jj++) {
1107:                 for (kk=kstart; kk<kend+1; kk++) {
1108:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1109:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1110:                   }
1111:                 }
1112:               }
1113:             }
1114:             rows[l] = l + nc*(slot);
1115:           }
1116:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1117:         }
1118:       }
1119:     }
1120:     PetscFree(values);
1121:     /* do not copy values to GPU since they are all zero and not yet needed there */
1122:     MatBindToCPU(J,PETSC_TRUE);
1123:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1124:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1125:     MatBindToCPU(J,PETSC_FALSE);
1126:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1127:   }
1128:   PetscFree2(rows,cols);
1129:   return(0);
1130: }

1132: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J,PetscBool isIS)
1133: {
1134:   PetscErrorCode         ierr;
1135:   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;
1136:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1137:   MPI_Comm               comm;
1138:   DMBoundaryType         bx,by;
1139:   ISLocalToGlobalMapping ltog,mltog;
1140:   DMDAStencilType        st;
1141:   PetscBool              removedups = PETSC_FALSE;

1144:   /*
1145:          nc - number of components per grid point
1146:          col - number of colors needed in one direction for single component problem

1148:   */
1149:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,NULL,NULL,&nc,&s,&bx,&by,NULL,&st);
1150:   if (!isIS && bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE) {
1151:     MatSetOption(J,MAT_SORTED_FULL,PETSC_TRUE);
1152:   }
1153:   col  = 2*s + 1;
1154:   /*
1155:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1156:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1157:   */
1158:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1159:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1160:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
1161:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
1162:   PetscObjectGetComm((PetscObject)da,&comm);

1164:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1165:   DMGetLocalToGlobalMapping(da,&ltog);

1167:   MatSetBlockSize(J,nc);
1168:   /* determine the matrix preallocation information */
1169:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1170:   for (i=xs; i<xs+nx; i++) {

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

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

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

1181:       cnt = 0;
1182:       for (k=0; k<nc; k++) {
1183:         for (l=lstart; l<lend+1; l++) {
1184:           for (p=pstart; p<pend+1; p++) {
1185:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1186:               cols[cnt++] = k + nc*(slot + gnx*l + p);
1187:             }
1188:           }
1189:         }
1190:         rows[k] = k + nc*(slot);
1191:       }
1192:       if (removedups) {
1193:         MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1194:       } else {
1195:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1196:       }
1197:     }
1198:   }
1199:   MatSetBlockSize(J,nc);
1200:   MatSeqAIJSetPreallocation(J,0,dnz);
1201:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1202:   MatPreallocateFinalize(dnz,onz);
1203:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1204:   if (!mltog) {
1205:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1206:   }

1208:   /*
1209:     For each node in the grid: we get the neighbors in the local (on processor ordering
1210:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1211:     PETSc ordering.
1212:   */
1213:   if (!da->prealloc_only) {
1214:     for (i=xs; i<xs+nx; i++) {

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

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

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

1225:         cnt = 0;
1226:         for (l=lstart; l<lend+1; l++) {
1227:           for (p=pstart; p<pend+1; p++) {
1228:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1229:               cols[cnt++] = nc*(slot + gnx*l + p);
1230:               for (k=1; k<nc; k++) {
1231:                 cols[cnt] = 1 + cols[cnt-1];cnt++;
1232:               }
1233:             }
1234:           }
1235:         }
1236:         for (k=0; k<nc; k++) rows[k] = k + nc*(slot);
1237:         MatSetValuesLocal(J,nc,rows,cnt,cols,NULL,INSERT_VALUES);
1238:       }
1239:     }
1240:     /* do not copy values to GPU since they are all zero and not yet needed there */
1241:     MatBindToCPU(J,PETSC_TRUE);
1242:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1243:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1244:     MatBindToCPU(J,PETSC_FALSE);
1245:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1246:     if (bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE) {
1247:       MatSetOption(J,MAT_SORTED_FULL,PETSC_FALSE);
1248:     }
1249:   }
1250:   PetscFree2(rows,cols);
1251:   return(0);
1252: }

1254: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
1255: {
1256:   PetscErrorCode         ierr;
1257:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1258:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
1259:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1260:   DM_DA                  *dd = (DM_DA*)da->data;
1261:   PetscInt               ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
1262:   MPI_Comm               comm;
1263:   DMBoundaryType         bx,by;
1264:   ISLocalToGlobalMapping ltog;
1265:   DMDAStencilType        st;
1266:   PetscBool              removedups = PETSC_FALSE;

1269:   /*
1270:          nc - number of components per grid point
1271:          col - number of colors needed in one direction for single component problem

1273:   */
1274:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,NULL,NULL,&nc,&s,&bx,&by,NULL,&st);
1275:   col  = 2*s + 1;
1276:   /*
1277:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1278:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1279:   */
1280:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1281:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1282:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
1283:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
1284:   PetscObjectGetComm((PetscObject)da,&comm);

1286:   PetscMalloc1(col*col*nc,&cols);
1287:   DMGetLocalToGlobalMapping(da,&ltog);

1289:   MatSetBlockSize(J,nc);
1290:   /* determine the matrix preallocation information */
1291:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1292:   for (i=xs; i<xs+nx; i++) {

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

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

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

1303:       for (k=0; k<nc; k++) {
1304:         cnt = 0;
1305:         for (l=lstart; l<lend+1; l++) {
1306:           for (p=pstart; p<pend+1; p++) {
1307:             if (l || p) {
1308:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1309:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1310:               }
1311:             } else {
1312:               if (dfill) {
1313:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1314:               } else {
1315:                 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1316:               }
1317:             }
1318:           }
1319:         }
1320:         row    = k + nc*(slot);
1321:         maxcnt = PetscMax(maxcnt,cnt);
1322:         if (removedups) {
1323:           MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1324:         } else {
1325:           MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1326:         }
1327:       }
1328:     }
1329:   }
1330:   MatSeqAIJSetPreallocation(J,0,dnz);
1331:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1332:   MatPreallocateFinalize(dnz,onz);
1333:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1335:   /*
1336:     For each node in the grid: we get the neighbors in the local (on processor ordering
1337:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1338:     PETSc ordering.
1339:   */
1340:   if (!da->prealloc_only) {
1341:     for (i=xs; i<xs+nx; i++) {

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

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

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

1352:         for (k=0; k<nc; k++) {
1353:           cnt = 0;
1354:           for (l=lstart; l<lend+1; l++) {
1355:             for (p=pstart; p<pend+1; p++) {
1356:               if (l || p) {
1357:                 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1358:                   for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1359:                 }
1360:               } else {
1361:                 if (dfill) {
1362:                   for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1363:                 } else {
1364:                   for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1365:                 }
1366:               }
1367:             }
1368:           }
1369:           row  = k + nc*(slot);
1370:           MatSetValuesLocal(J,1,&row,cnt,cols,NULL,INSERT_VALUES);
1371:         }
1372:       }
1373:     }
1374:     /* do not copy values to GPU since they are all zero and not yet needed there */
1375:     MatBindToCPU(J,PETSC_TRUE);
1376:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1377:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1378:     MatBindToCPU(J,PETSC_FALSE);
1379:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1380:   }
1381:   PetscFree(cols);
1382:   return(0);
1383: }

1385: /* ---------------------------------------------------------------------------------*/

1387: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J,PetscBool isIS)
1388: {
1389:   PetscErrorCode         ierr;
1390:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1391:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1392:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1393:   MPI_Comm               comm;
1394:   DMBoundaryType         bx,by,bz;
1395:   ISLocalToGlobalMapping ltog,mltog;
1396:   DMDAStencilType        st;
1397:   PetscBool              removedups = PETSC_FALSE;

1400:   /*
1401:          nc - number of components per grid point
1402:          col - number of colors needed in one direction for single component problem

1404:   */
1405:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1406:   if (!isIS && bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE && bz == DM_BOUNDARY_NONE) {
1407:     MatSetOption(J,MAT_SORTED_FULL,PETSC_TRUE);
1408:   }
1409:   col  = 2*s + 1;

1411:   /*
1412:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1413:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1414:   */
1415:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1416:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1417:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

1423:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1424:   DMGetLocalToGlobalMapping(da,&ltog);

1426:   MatSetBlockSize(J,nc);
1427:   /* determine the matrix preallocation information */
1428:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1429:   for (i=xs; i<xs+nx; i++) {
1430:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1431:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1432:     for (j=ys; j<ys+ny; j++) {
1433:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1434:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1435:       for (k=zs; k<zs+nz; k++) {
1436:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1437:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1441:         cnt = 0;
1442:         for (l=0; l<nc; l++) {
1443:           for (ii=istart; ii<iend+1; ii++) {
1444:             for (jj=jstart; jj<jend+1; jj++) {
1445:               for (kk=kstart; kk<kend+1; kk++) {
1446:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1447:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1448:                 }
1449:               }
1450:             }
1451:           }
1452:           rows[l] = l + nc*(slot);
1453:         }
1454:         if (removedups) {
1455:           MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1456:         } else {
1457:           MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1458:         }
1459:       }
1460:     }
1461:   }
1462:   MatSetBlockSize(J,nc);
1463:   MatSeqAIJSetPreallocation(J,0,dnz);
1464:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1465:   MatPreallocateFinalize(dnz,onz);
1466:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1467:   if (!mltog) {
1468:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1469:   }

1471:   /*
1472:     For each node in the grid: we get the neighbors in the local (on processor ordering
1473:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1474:     PETSc ordering.
1475:   */
1476:   if (!da->prealloc_only) {
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 (kk=kstart; kk<kend+1; kk++) {
1491:             for (jj=jstart; jj<jend+1; jj++) {
1492:               for (ii=istart; ii<iend+1; ii++) {
1493:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1494:                   cols[cnt++] = nc*(slot + ii + gnx*jj + gnx*gny*kk);
1495:                     for (l=1; l<nc; l++) {
1496:                       cols[cnt] = 1 + cols[cnt-1];cnt++;
1497:                   }
1498:                 }
1499:               }
1500:             }
1501:           }
1502:           rows[0] = nc*(slot); for (l=1; l<nc; l++) rows[l] = 1 + rows[l-1];
1503:           MatSetValuesLocal(J,nc,rows,cnt,cols,NULL,INSERT_VALUES);
1504:         }
1505:       }
1506:     }
1507:     /* do not copy values to GPU since they are all zero and not yet needed there */
1508:     MatBindToCPU(J,PETSC_TRUE);
1509:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1510:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1511:     if (!isIS && bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE && bz == DM_BOUNDARY_NONE) {
1512:       MatSetOption(J,MAT_SORTED_FULL,PETSC_FALSE);
1513:     }
1514:     MatBindToCPU(J,PETSC_FALSE);
1515:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1516:   }
1517:   PetscFree2(rows,cols);
1518:   return(0);
1519: }

1521: /* ---------------------------------------------------------------------------------*/

1523: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1524: {
1525:   PetscErrorCode         ierr;
1526:   DM_DA                  *dd = (DM_DA*)da->data;
1527:   PetscInt               xs,nx,i,j,gxs,gnx,row,k,l;
1528:   PetscInt               m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1529:   PetscInt               *ofill = dd->ofill,*dfill = dd->dfill;
1530:   DMBoundaryType         bx;
1531:   ISLocalToGlobalMapping ltog;
1532:   PetscMPIInt            rank,size;

1535:   MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1536:   MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);

1538:   /*
1539:          nc - number of components per grid point

1541:   */
1542:   DMDAGetInfo(da,&dim,&m,NULL,NULL,NULL,NULL,NULL,&nc,&s,&bx,NULL,NULL,NULL);
1543:   if (s > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Matrix creation for 1d not implemented correctly for stencil width larger than 1");
1544:   DMDAGetCorners(da,&xs,NULL,NULL,&nx,NULL,NULL);
1545:   DMDAGetGhostCorners(da,&gxs,NULL,NULL,&gnx,NULL,NULL);

1547:   MatSetBlockSize(J,nc);
1548:   PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);

1550:   /*
1551:         note should be smaller for first and last process with no periodic
1552:         does not handle dfill
1553:   */
1554:   cnt = 0;
1555:   /* coupling with process to the left */
1556:   for (i=0; i<s; i++) {
1557:     for (j=0; j<nc; j++) {
1558:       ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1559:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1560:       if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1561:         if (size > 1) ocols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1562:         else cols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1563:       }
1564:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1565:       cnt++;
1566:     }
1567:   }
1568:   for (i=s; i<nx-s; i++) {
1569:     for (j=0; j<nc; j++) {
1570:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1571:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1572:       cnt++;
1573:     }
1574:   }
1575:   /* coupling with process to the right */
1576:   for (i=nx-s; i<nx; i++) {
1577:     for (j=0; j<nc; j++) {
1578:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1579:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1580:       if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1581:         if (size > 1) ocols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1582:         else cols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1583:       }
1584:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1585:       cnt++;
1586:     }
1587:   }

1589:   MatSeqAIJSetPreallocation(J,0,cols);
1590:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1591:   PetscFree2(cols,ocols);

1593:   DMGetLocalToGlobalMapping(da,&ltog);
1594:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1596:   /*
1597:     For each node in the grid: we get the neighbors in the local (on processor ordering
1598:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1599:     PETSc ordering.
1600:   */
1601:   if (!da->prealloc_only) {
1602:     PetscMalloc1(maxcnt,&cols);
1603:     row = xs*nc;
1604:     /* coupling with process to the left */
1605:     for (i=xs; i<xs+s; i++) {
1606:       for (j=0; j<nc; j++) {
1607:         cnt = 0;
1608:         if (rank) {
1609:           for (l=0; l<s; l++) {
1610:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1611:           }
1612:         }
1613:         if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1614:           for (l=0; l<s; l++) {
1615:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (m + i - s - l)*nc + ofill[k];
1616:           }
1617:         }
1618:         if (dfill) {
1619:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1620:             cols[cnt++] = i*nc + dfill[k];
1621:           }
1622:         } else {
1623:           for (k=0; k<nc; k++) {
1624:             cols[cnt++] = i*nc + k;
1625:           }
1626:         }
1627:         for (l=0; l<s; l++) {
1628:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1629:         }
1630:         MatSetValues(J,1,&row,cnt,cols,NULL,INSERT_VALUES);
1631:         row++;
1632:       }
1633:     }
1634:     for (i=xs+s; i<xs+nx-s; i++) {
1635:       for (j=0; j<nc; j++) {
1636:         cnt = 0;
1637:         for (l=0; l<s; l++) {
1638:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1639:         }
1640:         if (dfill) {
1641:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1642:             cols[cnt++] = i*nc + dfill[k];
1643:           }
1644:         } else {
1645:           for (k=0; k<nc; k++) {
1646:             cols[cnt++] = i*nc + k;
1647:           }
1648:         }
1649:         for (l=0; l<s; l++) {
1650:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1651:         }
1652:         MatSetValues(J,1,&row,cnt,cols,NULL,INSERT_VALUES);
1653:         row++;
1654:       }
1655:     }
1656:     /* coupling with process to the right */
1657:     for (i=xs+nx-s; i<xs+nx; i++) {
1658:       for (j=0; j<nc; j++) {
1659:         cnt = 0;
1660:         for (l=0; l<s; l++) {
1661:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
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:         if (rank < size-1) {
1673:           for (l=0; l<s; l++) {
1674:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1675:           }
1676:         }
1677:         if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1678:           for (l=0; l<s; l++) {
1679:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s - l - m + 2)*nc + ofill[k];
1680:           }
1681:         }
1682:         MatSetValues(J,1,&row,cnt,cols,NULL,INSERT_VALUES);
1683:         row++;
1684:       }
1685:     }
1686:     PetscFree(cols);
1687:     /* do not copy values to GPU since they are all zero and not yet needed there */
1688:     MatBindToCPU(J,PETSC_TRUE);
1689:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1690:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1691:     MatBindToCPU(J,PETSC_FALSE);
1692:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1693:   }
1694:   return(0);
1695: }

1697: /* ---------------------------------------------------------------------------------*/

1699: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J,PetscBool isIS)
1700: {
1701:   PetscErrorCode         ierr;
1702:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1703:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1704:   PetscInt               istart,iend;
1705:   DMBoundaryType         bx;
1706:   ISLocalToGlobalMapping ltog,mltog;

1709:   /*
1710:          nc - number of components per grid point
1711:          col - number of colors needed in one direction for single component problem

1713:   */
1714:   DMDAGetInfo(da,&dim,&m,NULL,NULL,NULL,NULL,NULL,&nc,&s,&bx,NULL,NULL,NULL);
1715:   if (!isIS && bx == DM_BOUNDARY_NONE) {
1716:     MatSetOption(J,MAT_SORTED_FULL,PETSC_TRUE);
1717:   }
1718:   col  = 2*s + 1;

1720:   DMDAGetCorners(da,&xs,NULL,NULL,&nx,NULL,NULL);
1721:   DMDAGetGhostCorners(da,&gxs,NULL,NULL,&gnx,NULL,NULL);

1723:   MatSetBlockSize(J,nc);
1724:   MatSeqAIJSetPreallocation(J,col*nc,NULL);
1725:   MatMPIAIJSetPreallocation(J,col*nc,NULL,col*nc,NULL);

1727:   DMGetLocalToGlobalMapping(da,&ltog);
1728:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1729:   if (!mltog) {
1730:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1731:   }

1733:   /*
1734:     For each node in the grid: we get the neighbors in the local (on processor ordering
1735:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1736:     PETSc ordering.
1737:   */
1738:   if (!da->prealloc_only) {
1739:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1740:     for (i=xs; i<xs+nx; i++) {
1741:       istart = PetscMax(-s,gxs - i);
1742:       iend   = PetscMin(s,gxs + gnx - i - 1);
1743:       slot   = i - gxs;

1745:       cnt = 0;
1746:       for (i1=istart; i1<iend+1; i1++) {
1747:         cols[cnt++] = nc*(slot + i1);
1748:         for (l=1; l<nc; l++) {
1749:           cols[cnt] = 1 + cols[cnt-1];cnt++;
1750:         }
1751:       }
1752:       rows[0] = nc*(slot); for (l=1; l<nc; l++) rows[l] = 1 + rows[l-1];
1753:       MatSetValuesLocal(J,nc,rows,cnt,cols,NULL,INSERT_VALUES);
1754:     }
1755:     /* do not copy values to GPU since they are all zero and not yet needed there */
1756:     MatBindToCPU(J,PETSC_TRUE);
1757:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1758:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1759:     if (!isIS && bx == DM_BOUNDARY_NONE) {
1760:       MatSetOption(J,MAT_SORTED_FULL,PETSC_FALSE);
1761:     }
1762:     MatBindToCPU(J,PETSC_FALSE);
1763:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1764:     PetscFree2(rows,cols);
1765:   }
1766:   return(0);
1767: }

1769: /* ---------------------------------------------------------------------------------*/

1771: PetscErrorCode DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(DM da,Mat J,PetscBool isIS)
1772: {
1773:   PetscErrorCode         ierr;
1774:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1775:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1776:   PetscInt               istart,iend;
1777:   DMBoundaryType         bx;
1778:   ISLocalToGlobalMapping ltog,mltog;

1781:   /*
1782:          nc - number of components per grid point
1783:          col - number of colors needed in one direction for single component problem
1784:   */
1785:   DMDAGetInfo(da,&dim,&m,NULL,NULL,NULL,NULL,NULL,&nc,&s,&bx,NULL,NULL,NULL);
1786:   col  = 2*s + 1;

1788:   DMDAGetCorners(da,&xs,NULL,NULL,&nx,NULL,NULL);
1789:   DMDAGetGhostCorners(da,&gxs,NULL,NULL,&gnx,NULL,NULL);

1791:   MatSetBlockSize(J,nc);
1792:   MatSeqAIJSetTotalPreallocation(J,nx*nc*col*nc);

1794:   DMGetLocalToGlobalMapping(da,&ltog);
1795:   MatGetLocalToGlobalMapping(J,&mltog,NULL);
1796:   if (!mltog) {
1797:     MatSetLocalToGlobalMapping(J,ltog,ltog);
1798:   }

1800:   /*
1801:     For each node in the grid: we get the neighbors in the local (on processor ordering
1802:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1803:     PETSc ordering.
1804:   */
1805:   if (!da->prealloc_only) {
1806:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1807:     for (i=xs; i<xs+nx; i++) {
1808:       istart = PetscMax(-s,gxs - i);
1809:       iend   = PetscMin(s,gxs + gnx - i - 1);
1810:       slot   = i - gxs;

1812:       cnt = 0;
1813:       for (i1=istart; i1<iend+1; i1++) {
1814:         cols[cnt++] = nc*(slot + i1);
1815:         for (l=1; l<nc; l++) {
1816:           cols[cnt] = 1 + cols[cnt-1];cnt++;
1817:         }
1818:       }
1819:       rows[0] = nc*(slot); for (l=1; l<nc; l++) rows[l] = 1 + rows[l-1];
1820:       MatSetValuesLocal(J,nc,rows,cnt,cols,NULL,INSERT_VALUES);
1821:     }
1822:     /* do not copy values to GPU since they are all zero and not yet needed there */
1823:     MatBindToCPU(J,PETSC_TRUE);
1824:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1825:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1826:     if (!isIS && bx == DM_BOUNDARY_NONE) {
1827:       MatSetOption(J,MAT_SORTED_FULL,PETSC_FALSE);
1828:     }
1829:     MatBindToCPU(J,PETSC_FALSE);
1830:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1831:     PetscFree2(rows,cols);
1832:   }
1833:   MatSetOption(J,MAT_SORTED_FULL,PETSC_FALSE);
1834:   return(0);
1835: }

1837: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1838: {
1839:   PetscErrorCode         ierr;
1840:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1841:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1842:   PetscInt               istart,iend,jstart,jend,ii,jj;
1843:   MPI_Comm               comm;
1844:   PetscScalar            *values;
1845:   DMBoundaryType         bx,by;
1846:   DMDAStencilType        st;
1847:   ISLocalToGlobalMapping ltog;

1850:   /*
1851:      nc - number of components per grid point
1852:      col - number of colors needed in one direction for single component problem
1853:   */
1854:   DMDAGetInfo(da,&dim,&m,&n,NULL,NULL,NULL,NULL,&nc,&s,&bx,&by,NULL,&st);
1855:   col  = 2*s + 1;

1857:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
1858:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
1859:   PetscObjectGetComm((PetscObject)da,&comm);

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

1863:   DMGetLocalToGlobalMapping(da,&ltog);

1865:   /* determine the matrix preallocation information */
1866:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
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);

1875:       /* Find block columns in block row */
1876:       cnt = 0;
1877:       for (ii=istart; ii<iend+1; ii++) {
1878:         for (jj=jstart; jj<jend+1; jj++) {
1879:           if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1880:             cols[cnt++] = slot + ii + gnx*jj;
1881:           }
1882:         }
1883:       }
1884:       MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1885:     }
1886:   }
1887:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1888:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1889:   MatPreallocateFinalize(dnz,onz);

1891:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1893:   /*
1894:     For each node in the grid: we get the neighbors in the local (on processor ordering
1895:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1896:     PETSc ordering.
1897:   */
1898:   if (!da->prealloc_only) {
1899:     PetscCalloc1(col*col*nc*nc,&values);
1900:     for (i=xs; i<xs+nx; i++) {
1901:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1902:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1903:       for (j=ys; j<ys+ny; j++) {
1904:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1905:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1906:         slot = i - gxs + gnx*(j - gys);
1907:         cnt  = 0;
1908:         for (ii=istart; ii<iend+1; ii++) {
1909:           for (jj=jstart; jj<jend+1; jj++) {
1910:             if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1911:               cols[cnt++] = slot + ii + gnx*jj;
1912:             }
1913:           }
1914:         }
1915:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1916:       }
1917:     }
1918:     PetscFree(values);
1919:     /* do not copy values to GPU since they are all zero and not yet needed there */
1920:     MatBindToCPU(J,PETSC_TRUE);
1921:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1922:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1923:     MatBindToCPU(J,PETSC_FALSE);
1924:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1925:   }
1926:   PetscFree(cols);
1927:   return(0);
1928: }

1930: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1931: {
1932:   PetscErrorCode         ierr;
1933:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1934:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1935:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1936:   MPI_Comm               comm;
1937:   PetscScalar            *values;
1938:   DMBoundaryType         bx,by,bz;
1939:   DMDAStencilType        st;
1940:   ISLocalToGlobalMapping ltog;

1943:   /*
1944:          nc - number of components per grid point
1945:          col - number of colors needed in one direction for single component problem

1947:   */
1948:   DMDAGetInfo(da,&dim,&m,&n,&p,NULL,NULL,NULL,&nc,&s,&bx,&by,&bz,&st);
1949:   col  = 2*s + 1;

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

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

1957:   DMGetLocalToGlobalMapping(da,&ltog);

1959:   /* determine the matrix preallocation information */
1960:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1961:   for (i=xs; i<xs+nx; i++) {
1962:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1963:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1964:     for (j=ys; j<ys+ny; j++) {
1965:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1966:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1967:       for (k=zs; k<zs+nz; k++) {
1968:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1969:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1973:         /* Find block columns in block row */
1974:         cnt = 0;
1975:         for (ii=istart; ii<iend+1; ii++) {
1976:           for (jj=jstart; jj<jend+1; jj++) {
1977:             for (kk=kstart; kk<kend+1; kk++) {
1978:               if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1979:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1980:               }
1981:             }
1982:           }
1983:         }
1984:         MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1985:       }
1986:     }
1987:   }
1988:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1989:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1990:   MatPreallocateFinalize(dnz,onz);

1992:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1994:   /*
1995:     For each node in the grid: we get the neighbors in the local (on processor ordering
1996:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1997:     PETSc ordering.
1998:   */
1999:   if (!da->prealloc_only) {
2000:     PetscCalloc1(col*col*col*nc*nc,&values);
2001:     for (i=xs; i<xs+nx; i++) {
2002:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2003:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2004:       for (j=ys; j<ys+ny; j++) {
2005:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2006:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2007:         for (k=zs; k<zs+nz; k++) {
2008:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2009:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2013:           cnt = 0;
2014:           for (ii=istart; ii<iend+1; ii++) {
2015:             for (jj=jstart; jj<jend+1; jj++) {
2016:               for (kk=kstart; kk<kend+1; kk++) {
2017:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2018:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2019:                 }
2020:               }
2021:             }
2022:           }
2023:           MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2024:         }
2025:       }
2026:     }
2027:     PetscFree(values);
2028:     /* do not copy values to GPU since they are all zero and not yet needed there */
2029:     MatBindToCPU(J,PETSC_TRUE);
2030:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2031:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2032:     MatBindToCPU(J,PETSC_FALSE);
2033:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2034:   }
2035:   PetscFree(cols);
2036:   return(0);
2037: }

2039: /*
2040:   This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
2041:   identify in the local ordering with periodic domain.
2042: */
2043: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
2044: {
2046:   PetscInt       i,n;

2049:   ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
2050:   ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
2051:   for (i=0,n=0; i<*cnt; i++) {
2052:     if (col[i] >= *row) col[n++] = col[i];
2053:   }
2054:   *cnt = n;
2055:   return(0);
2056: }

2058: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
2059: {
2060:   PetscErrorCode         ierr;
2061:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2062:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
2063:   PetscInt               istart,iend,jstart,jend,ii,jj;
2064:   MPI_Comm               comm;
2065:   PetscScalar            *values;
2066:   DMBoundaryType         bx,by;
2067:   DMDAStencilType        st;
2068:   ISLocalToGlobalMapping ltog;

2071:   /*
2072:      nc - number of components per grid point
2073:      col - number of colors needed in one direction for single component problem
2074:   */
2075:   DMDAGetInfo(da,&dim,&m,&n,NULL,NULL,NULL,NULL,&nc,&s,&bx,&by,NULL,&st);
2076:   col  = 2*s + 1;

2078:   DMDAGetCorners(da,&xs,&ys,NULL,&nx,&ny,NULL);
2079:   DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gnx,&gny,NULL);
2080:   PetscObjectGetComm((PetscObject)da,&comm);

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

2084:   DMGetLocalToGlobalMapping(da,&ltog);

2086:   /* determine the matrix preallocation information */
2087:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
2088:   for (i=xs; i<xs+nx; i++) {
2089:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2090:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2091:     for (j=ys; j<ys+ny; j++) {
2092:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2093:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2094:       slot   = i - gxs + gnx*(j - gys);

2096:       /* Find block columns in block row */
2097:       cnt = 0;
2098:       for (ii=istart; ii<iend+1; ii++) {
2099:         for (jj=jstart; jj<jend+1; jj++) {
2100:           if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2101:             cols[cnt++] = slot + ii + gnx*jj;
2102:           }
2103:         }
2104:       }
2105:       L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2106:       MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2107:     }
2108:   }
2109:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2110:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2111:   MatPreallocateFinalize(dnz,onz);

2113:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2115:   /*
2116:     For each node in the grid: we get the neighbors in the local (on processor ordering
2117:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2118:     PETSc ordering.
2119:   */
2120:   if (!da->prealloc_only) {
2121:     PetscCalloc1(col*col*nc*nc,&values);
2122:     for (i=xs; i<xs+nx; i++) {
2123:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2124:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2125:       for (j=ys; j<ys+ny; j++) {
2126:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2127:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2128:         slot   = i - gxs + gnx*(j - gys);

2130:         /* Find block columns in block row */
2131:         cnt = 0;
2132:         for (ii=istart; ii<iend+1; ii++) {
2133:           for (jj=jstart; jj<jend+1; jj++) {
2134:             if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2135:               cols[cnt++] = slot + ii + gnx*jj;
2136:             }
2137:           }
2138:         }
2139:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2140:         MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2141:       }
2142:     }
2143:     PetscFree(values);
2144:     /* do not copy values to GPU since they are all zero and not yet needed there */
2145:     MatBindToCPU(J,PETSC_TRUE);
2146:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2147:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2148:     MatBindToCPU(J,PETSC_FALSE);
2149:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2150:   }
2151:   PetscFree(cols);
2152:   return(0);
2153: }

2155: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
2156: {
2157:   PetscErrorCode         ierr;
2158:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2159:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
2160:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
2161:   MPI_Comm               comm;
2162:   PetscScalar            *values;
2163:   DMBoundaryType         bx,by,bz;
2164:   DMDAStencilType        st;
2165:   ISLocalToGlobalMapping ltog;

2168:   /*
2169:      nc - number of components per grid point
2170:      col - number of colors needed in one direction for single component problem
2171:   */
2172:   DMDAGetInfo(da,&dim,&m,&n,&p,NULL,NULL,NULL,&nc,&s,&bx,&by,&bz,&st);
2173:   col  = 2*s + 1;

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

2179:   /* create the matrix */
2180:   PetscMalloc1(col*col*col,&cols);

2182:   DMGetLocalToGlobalMapping(da,&ltog);

2184:   /* determine the matrix preallocation information */
2185:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2186:   for (i=xs; i<xs+nx; i++) {
2187:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2188:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2189:     for (j=ys; j<ys+ny; j++) {
2190:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2191:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2192:       for (k=zs; k<zs+nz; k++) {
2193:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2194:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2198:         /* Find block columns in block row */
2199:         cnt = 0;
2200:         for (ii=istart; ii<iend+1; ii++) {
2201:           for (jj=jstart; jj<jend+1; jj++) {
2202:             for (kk=kstart; kk<kend+1; kk++) {
2203:               if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2204:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2205:               }
2206:             }
2207:           }
2208:         }
2209:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2210:         MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2211:       }
2212:     }
2213:   }
2214:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2215:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2216:   MatPreallocateFinalize(dnz,onz);

2218:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2220:   /*
2221:     For each node in the grid: we get the neighbors in the local (on processor ordering
2222:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2223:     PETSc ordering.
2224:   */
2225:   if (!da->prealloc_only) {
2226:     PetscCalloc1(col*col*col*nc*nc,&values);
2227:     for (i=xs; i<xs+nx; i++) {
2228:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2229:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2230:       for (j=ys; j<ys+ny; j++) {
2231:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2232:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2233:         for (k=zs; k<zs+nz; k++) {
2234:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2235:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2239:           cnt = 0;
2240:           for (ii=istart; ii<iend+1; ii++) {
2241:             for (jj=jstart; jj<jend+1; jj++) {
2242:               for (kk=kstart; kk<kend+1; kk++) {
2243:                 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2244:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2245:                 }
2246:               }
2247:             }
2248:           }
2249:           L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2250:           MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2251:         }
2252:       }
2253:     }
2254:     PetscFree(values);
2255:     /* do not copy values to GPU since they are all zero and not yet needed there */
2256:     MatBindToCPU(J,PETSC_TRUE);
2257:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2258:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2259:     MatBindToCPU(J,PETSC_FALSE);
2260:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2261:   }
2262:   PetscFree(cols);
2263:   return(0);
2264: }

2266: /* ---------------------------------------------------------------------------------*/

2268: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2269: {
2270:   PetscErrorCode         ierr;
2271:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2272:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2273:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2274:   DM_DA                  *dd = (DM_DA*)da->data;
2275:   PetscInt               ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2276:   MPI_Comm               comm;
2277:   PetscScalar            *values;
2278:   DMBoundaryType         bx,by,bz;
2279:   ISLocalToGlobalMapping ltog;
2280:   DMDAStencilType        st;
2281:   PetscBool              removedups = PETSC_FALSE;

2284:   /*
2285:          nc - number of components per grid point
2286:          col - number of colors needed in one direction for single component problem

2288:   */
2289:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2290:   col  = 2*s + 1;
2291:   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\
2292:                  by 2*stencil_width + 1\n");
2293:   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\
2294:                  by 2*stencil_width + 1\n");
2295:   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\
2296:                  by 2*stencil_width + 1\n");

2298:   /*
2299:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2300:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2301:   */
2302:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2303:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2304:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

2310:   PetscMalloc1(col*col*col*nc,&cols);
2311:   DMGetLocalToGlobalMapping(da,&ltog);

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

2316:   MatSetBlockSize(J,nc);
2317:   for (i=xs; i<xs+nx; i++) {
2318:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2319:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2320:     for (j=ys; j<ys+ny; j++) {
2321:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2322:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2323:       for (k=zs; k<zs+nz; k++) {
2324:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2325:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2329:         for (l=0; l<nc; l++) {
2330:           cnt = 0;
2331:           for (ii=istart; ii<iend+1; ii++) {
2332:             for (jj=jstart; jj<jend+1; jj++) {
2333:               for (kk=kstart; kk<kend+1; kk++) {
2334:                 if (ii || jj || kk) {
2335:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2336:                     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);
2337:                   }
2338:                 } else {
2339:                   if (dfill) {
2340:                     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);
2341:                   } else {
2342:                     for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2343:                   }
2344:                 }
2345:               }
2346:             }
2347:           }
2348:           row  = l + nc*(slot);
2349:           maxcnt = PetscMax(maxcnt,cnt);
2350:           if (removedups) {
2351:             MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2352:           } else {
2353:             MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2354:           }
2355:         }
2356:       }
2357:     }
2358:   }
2359:   MatSeqAIJSetPreallocation(J,0,dnz);
2360:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2361:   MatPreallocateFinalize(dnz,onz);
2362:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2364:   /*
2365:     For each node in the grid: we get the neighbors in the local (on processor ordering
2366:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2367:     PETSc ordering.
2368:   */
2369:   if (!da->prealloc_only) {
2370:     PetscCalloc1(maxcnt,&values);
2371:     for (i=xs; i<xs+nx; i++) {
2372:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2373:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2374:       for (j=ys; j<ys+ny; j++) {
2375:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2376:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2377:         for (k=zs; k<zs+nz; k++) {
2378:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2379:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2383:           for (l=0; l<nc; l++) {
2384:             cnt = 0;
2385:             for (ii=istart; ii<iend+1; ii++) {
2386:               for (jj=jstart; jj<jend+1; jj++) {
2387:                 for (kk=kstart; kk<kend+1; kk++) {
2388:                   if (ii || jj || kk) {
2389:                     if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2390:                       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);
2391:                     }
2392:                   } else {
2393:                     if (dfill) {
2394:                       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);
2395:                     } else {
2396:                       for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2397:                     }
2398:                   }
2399:                 }
2400:               }
2401:             }
2402:             row  = l + nc*(slot);
2403:             MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2404:           }
2405:         }
2406:       }
2407:     }
2408:     PetscFree(values);
2409:     /* do not copy values to GPU since they are all zero and not yet needed there */
2410:     MatBindToCPU(J,PETSC_TRUE);
2411:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2412:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2413:     MatBindToCPU(J,PETSC_FALSE);
2414:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2415:   }
2416:   PetscFree(cols);
2417:   return(0);
2418: }