Actual source code: fdda.c

petsc-3.9.2 2018-05-20
Report Typos and Errors

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

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

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

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

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

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

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

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

 55:     Logically Collective on DMDA

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


 63:     Level: developer

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

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

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

 80:    Contributed by Glenn Hammond

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

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

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

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


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

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

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

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

145:   PetscObjectGetComm((PetscObject)da,&comm);
146:   MPI_Comm_size(comm,&size);
147:   if (ctype == IS_COLORING_LOCAL) {
148:     if (size == 1) {
149:       ctype = IS_COLORING_GLOBAL;
150:     } else if (dim > 1) {
151:       if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
152:         SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_LOCAL cannot be used for periodic boundary condition having both ends of the domain  on the same process");
153:       }
154:     }
155:   }

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

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

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

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

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

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

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

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

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

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

268: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
269: {
270:   PetscErrorCode   ierr;
271:   PetscInt         xs,ys,nx,ny,i,j,gxs,gys,gnx,gny,m,n,p,dim,s,k,nc,col,zs,gzs,ii,l,nz,gnz,M,N,P;
272:   PetscInt         ncolors;
273:   MPI_Comm         comm;
274:   DMBoundaryType bx,by,bz;
275:   DMDAStencilType  st;
276:   ISColoringValue  *colors;
277:   DM_DA            *dd = (DM_DA*)da->data;

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

284:   */
285:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
286:   col  = 2*s + 1;
287:   if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
288:                                                          by 2*stencil_width + 1\n");
289:   if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
290:                                                          by 2*stencil_width + 1\n");
291:   if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
292:                                                          by 2*stencil_width + 1\n");

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

496:    Logically Collective on Mat

498:    Input Parameters:
499: +  mat - the matrix
500: -  da - the da

502:    Level: intermediate

504: @*/
505: PetscErrorCode MatSetupDM(Mat mat,DM da)
506: {

512:   PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
513:   return(0);
514: }

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

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

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

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

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

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

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

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

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

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

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

608:   MatInitializePackage();
609:   mtype = da->mattype;

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

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

637:       if (bs < 0) {
638:         if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
639:           PetscInt pStart, pEnd, p, dof;

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

678:   /*
679:          nc - number of components per grid point
680:          col - number of colors needed in one direction for single component problem

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

705:      We could switch based on Atype (or mtype), but we do not since the
706:    specialized setting routines depend only the particular preallocation
707:    details of the matrix, not the type itself.
708:   */
709:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
710:   if (!aij) {
711:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
712:   }
713:   if (!aij) {
714:     PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
715:     if (!baij) {
716:       PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
717:     }
718:     if (!baij) {
719:       PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
720:       if (!sbaij) {
721:         PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
722:       }
723:       if (!sbaij) {
724:         PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);
725:         if (!sell) {
726:           PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);
727:         }
728:       }
729:       if (!sell) {
730:         PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
731:       }
732:     }
733:   }
734:   if (aij) {
735:     if (dim == 1) {
736:       if (dd->ofill) {
737:         DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
738:       } else {
739:         DMCreateMatrix_DA_1d_MPIAIJ(da,A);
740:       }
741:     } else if (dim == 2) {
742:       if (dd->ofill) {
743:         DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
744:       } else {
745:         DMCreateMatrix_DA_2d_MPIAIJ(da,A);
746:       }
747:     } else if (dim == 3) {
748:       if (dd->ofill) {
749:         DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
750:       } else {
751:         DMCreateMatrix_DA_3d_MPIAIJ(da,A);
752:       }
753:     }
754:   } else if (baij) {
755:     if (dim == 2) {
756:       DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
757:     } else if (dim == 3) {
758:       DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
759:     } 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);
760:   } else if (sbaij) {
761:     if (dim == 2) {
762:       DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
763:     } else if (dim == 3) {
764:       DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
765:     } 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);
766:   } else if (sell) {
767:      if (dim == 2) {
768:        DMCreateMatrix_DA_2d_MPISELL(da,A);
769:      } else if (dim == 3) {
770:        DMCreateMatrix_DA_3d_MPISELL(da,A);
771:      } 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);
772:   } else if (is) {
773:     DMCreateMatrix_DA_IS(da,A);
774:   } else {
775:     ISLocalToGlobalMapping ltog;

777:     MatSetBlockSize(A,dof);
778:     MatSetUp(A);
779:     DMGetLocalToGlobalMapping(da,&ltog);
780:     MatSetLocalToGlobalMapping(A,ltog,ltog);
781:   }
782:   DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
783:   MatSetStencil(A,dim,dims,starts,dof);
784:   MatSetDM(A,da);
785:   MPI_Comm_size(comm,&size);
786:   if (size > 1) {
787:     /* change viewer to display matrix in natural ordering */
788:     MatSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
789:     MatSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
790:   }
791:   MatSetFromOptions(A);
792:   *J = A;
793:   return(0);
794: }

796: /* ---------------------------------------------------------------------------------*/
797: PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J)
798: {
799:   DM_DA                  *da = (DM_DA*)dm->data;
800:   Mat                    lJ;
801:   ISLocalToGlobalMapping ltog;
802:   IS                     is_loc_filt, is_glob;
803:   const PetscInt         *e_loc,*idx;
804:   PetscInt               i,nel,nen,dnz,nv,dof,dim,*gidx,nb;
805:   PetscErrorCode         ierr;

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

814:   MatSetBlockSize(J,dof);

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

821:   /* obtain a consistent local ordering for MATIS */
822:   ISSortRemoveDups(is_loc_filt);
823:   ISBlockGetLocalSize(is_loc_filt,&nb);
824:   DMGetLocalToGlobalMapping(dm,&ltog);
825:   ISLocalToGlobalMappingGetSize(ltog,&nv);
826:   PetscMalloc1(PetscMax(nb,nv/dof),&gidx);
827:   ISBlockGetIndices(is_loc_filt,&idx);
828:   ISLocalToGlobalMappingApplyBlock(ltog,nb,idx,gidx);
829:   ISBlockRestoreIndices(is_loc_filt,&idx);
830:   ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nb,gidx,PETSC_USE_POINTER,&is_glob);
831:   ISLocalToGlobalMappingCreateIS(is_glob,&ltog);
832:   ISDestroy(&is_glob);
833:   MatSetLocalToGlobalMapping(J,ltog,ltog);
834:   ISLocalToGlobalMappingDestroy(&ltog);

836:   /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */
837:   MatISGetLocalMat(J,&lJ);
838:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
839:   ISDestroy(&is_loc_filt);
840:   ISCreateStride(PetscObjectComm((PetscObject)lJ),nv/dof,0,1,&is_glob);
841:   ISGetIndices(is_glob,&idx);
842:   ISGlobalToLocalMappingApplyBlock(ltog,IS_GTOLM_MASK,nv/dof,idx,&nb,gidx);
843:   ISRestoreIndices(is_glob,&idx);
844:   ISDestroy(&is_glob);
845:   ISLocalToGlobalMappingDestroy(&ltog);
846:   ISCreateBlock(PETSC_COMM_SELF,dof,nb,gidx,PETSC_USE_POINTER,&is_loc_filt);
847:   ISLocalToGlobalMappingCreateIS(is_loc_filt,&ltog);
848:   ISDestroy(&is_loc_filt);
849:   MatSetLocalToGlobalMapping(lJ,ltog,ltog);
850:   ISLocalToGlobalMappingDestroy(&ltog);
851:   PetscFree(gidx);

853:   /* Preallocation (not exact) */
854:   switch (da->elementtype) {
855:   case DMDA_ELEMENT_P1:
856:   case DMDA_ELEMENT_Q1:
857:     dnz = 1;
858:     for (i=0; i<dim; i++) dnz *= 3;
859:     dnz *= dof;
860:     break;
861:   default:
862:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled element type %d",da->elementtype);
863:     break;
864:   }
865:   MatSeqAIJSetPreallocation(lJ,dnz,NULL);
866:   MatSeqBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);
867:   MatSeqSBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);
868:   MatISRestoreLocalMat(J,&lJ);
869:   return(0);
870: }

872: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J)
873: {
874:   PetscErrorCode         ierr;
875:   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;
876:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
877:   MPI_Comm               comm;
878:   PetscScalar            *values;
879:   DMBoundaryType         bx,by;
880:   ISLocalToGlobalMapping ltog;
881:   DMDAStencilType        st;

884:   /*
885:          nc - number of components per grid point
886:          col - number of colors needed in one direction for single component problem

888:   */
889:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
890:   col  = 2*s + 1;
891:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
892:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
893:   PetscObjectGetComm((PetscObject)da,&comm);

895:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
896:   DMGetLocalToGlobalMapping(da,&ltog);

898:   MatSetBlockSize(J,nc);
899:   /* determine the matrix preallocation information */
900:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
901:   for (i=xs; i<xs+nx; i++) {

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

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

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

912:       cnt = 0;
913:       for (k=0; k<nc; k++) {
914:         for (l=lstart; l<lend+1; l++) {
915:           for (p=pstart; p<pend+1; p++) {
916:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
917:               cols[cnt++] = k + nc*(slot + gnx*l + p);
918:             }
919:           }
920:         }
921:         rows[k] = k + nc*(slot);
922:       }
923:       MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
924:     }
925:   }
926:   MatSetBlockSize(J,nc);
927:   MatSeqSELLSetPreallocation(J,0,dnz);
928:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
929:   MatPreallocateFinalize(dnz,onz);

931:   MatSetLocalToGlobalMapping(J,ltog,ltog);

933:   /*
934:     For each node in the grid: we get the neighbors in the local (on processor ordering
935:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
936:     PETSc ordering.
937:   */
938:   if (!da->prealloc_only) {
939:     PetscCalloc1(col*col*nc*nc,&values);
940:     for (i=xs; i<xs+nx; i++) {

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

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

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

951:         cnt = 0;
952:         for (k=0; k<nc; k++) {
953:           for (l=lstart; l<lend+1; l++) {
954:             for (p=pstart; p<pend+1; p++) {
955:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
956:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
957:               }
958:             }
959:           }
960:           rows[k] = k + nc*(slot);
961:         }
962:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
963:       }
964:     }
965:     PetscFree(values);
966:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
967:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
968:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
969:   }
970:   PetscFree2(rows,cols);
971:   return(0);
972: }

974: PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da,Mat J)
975: {
976:   PetscErrorCode         ierr;
977:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
978:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
979:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
980:   MPI_Comm               comm;
981:   PetscScalar            *values;
982:   DMBoundaryType         bx,by,bz;
983:   ISLocalToGlobalMapping ltog;
984:   DMDAStencilType        st;

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

991:   */
992:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
993:   col  = 2*s + 1;
994:   DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
995:   DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
996:   PetscObjectGetComm((PetscObject)da,&comm);

998:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
999:   DMGetLocalToGlobalMapping(da,&ltog);

1001:   MatSetBlockSize(J,nc);
1002:   /* determine the matrix preallocation information */
1003:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1004:   for (i=xs; i<xs+nx; i++) {
1005:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1006:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1007:     for (j=ys; j<ys+ny; j++) {
1008:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1009:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1010:       for (k=zs; k<zs+nz; k++) {
1011:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1012:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1016:         cnt = 0;
1017:         for (l=0; l<nc; l++) {
1018:           for (ii=istart; ii<iend+1; ii++) {
1019:             for (jj=jstart; jj<jend+1; jj++) {
1020:               for (kk=kstart; kk<kend+1; kk++) {
1021:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1022:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1023:                 }
1024:               }
1025:             }
1026:           }
1027:           rows[l] = l + nc*(slot);
1028:         }
1029:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1030:       }
1031:     }
1032:   }
1033:   MatSetBlockSize(J,nc);
1034:   MatSeqSELLSetPreallocation(J,0,dnz);
1035:   MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1036:   MatPreallocateFinalize(dnz,onz);
1037:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1039:   /*
1040:     For each node in the grid: we get the neighbors in the local (on processor ordering
1041:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1042:     PETSc ordering.
1043:   */
1044:   if (!da->prealloc_only) {
1045:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1046:     for (i=xs; i<xs+nx; i++) {
1047:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1048:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1049:       for (j=ys; j<ys+ny; j++) {
1050:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1051:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1052:         for (k=zs; k<zs+nz; k++) {
1053:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1054:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1058:           cnt = 0;
1059:           for (l=0; l<nc; l++) {
1060:             for (ii=istart; ii<iend+1; ii++) {
1061:               for (jj=jstart; jj<jend+1; jj++) {
1062:                 for (kk=kstart; kk<kend+1; kk++) {
1063:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1064:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1065:                   }
1066:                 }
1067:               }
1068:             }
1069:             rows[l] = l + nc*(slot);
1070:           }
1071:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1072:         }
1073:       }
1074:     }
1075:     PetscFree(values);
1076:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1077:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1078:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1079:   }
1080:   PetscFree2(rows,cols);
1081:   return(0);
1082: }

1084: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
1085: {
1086:   PetscErrorCode         ierr;
1087:   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;
1088:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1089:   MPI_Comm               comm;
1090:   PetscScalar            *values;
1091:   DMBoundaryType         bx,by;
1092:   ISLocalToGlobalMapping ltog;
1093:   DMDAStencilType        st;
1094:   PetscBool              removedups = PETSC_FALSE;

1097:   /*
1098:          nc - number of components per grid point
1099:          col - number of colors needed in one direction for single component problem

1101:   */
1102:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1103:   col  = 2*s + 1;
1104:   /*
1105:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1106:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1107:   */
1108:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1109:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1110:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1111:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1112:   PetscObjectGetComm((PetscObject)da,&comm);

1114:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1115:   DMGetLocalToGlobalMapping(da,&ltog);

1117:   MatSetBlockSize(J,nc);
1118:   /* determine the matrix preallocation information */
1119:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1120:   for (i=xs; i<xs+nx; i++) {

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

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

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

1131:       cnt = 0;
1132:       for (k=0; k<nc; k++) {
1133:         for (l=lstart; l<lend+1; l++) {
1134:           for (p=pstart; p<pend+1; p++) {
1135:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1136:               cols[cnt++] = k + nc*(slot + gnx*l + p);
1137:             }
1138:           }
1139:         }
1140:         rows[k] = k + nc*(slot);
1141:       }
1142:       if (removedups) {
1143:         MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1144:       } else {
1145:         MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1146:       }
1147:     }
1148:   }
1149:   MatSetBlockSize(J,nc);
1150:   MatSeqAIJSetPreallocation(J,0,dnz);
1151:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1152:   MatPreallocateFinalize(dnz,onz);

1154:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1156:   /*
1157:     For each node in the grid: we get the neighbors in the local (on processor ordering
1158:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1159:     PETSc ordering.
1160:   */
1161:   if (!da->prealloc_only) {
1162:     PetscCalloc1(col*col*nc*nc,&values);
1163:     for (i=xs; i<xs+nx; i++) {

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

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

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

1174:         cnt = 0;
1175:         for (k=0; k<nc; k++) {
1176:           for (l=lstart; l<lend+1; l++) {
1177:             for (p=pstart; p<pend+1; p++) {
1178:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
1179:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
1180:               }
1181:             }
1182:           }
1183:           rows[k] = k + nc*(slot);
1184:         }
1185:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1186:       }
1187:     }
1188:     PetscFree(values);
1189:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1190:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1191:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1192:   }
1193:   PetscFree2(rows,cols);
1194:   return(0);
1195: }

1197: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
1198: {
1199:   PetscErrorCode         ierr;
1200:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1201:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
1202:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
1203:   DM_DA                  *dd = (DM_DA*)da->data;
1204:   PetscInt               ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
1205:   MPI_Comm               comm;
1206:   PetscScalar            *values;
1207:   DMBoundaryType         bx,by;
1208:   ISLocalToGlobalMapping ltog;
1209:   DMDAStencilType        st;
1210:   PetscBool              removedups = PETSC_FALSE;

1213:   /*
1214:          nc - number of components per grid point
1215:          col - number of colors needed in one direction for single component problem

1217:   */
1218:   DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1219:   col  = 2*s + 1;
1220:   /*
1221:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1222:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1223:   */
1224:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1225:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1226:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1227:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1228:   PetscObjectGetComm((PetscObject)da,&comm);

1230:   PetscMalloc1(col*col*nc,&cols);
1231:   DMGetLocalToGlobalMapping(da,&ltog);

1233:   MatSetBlockSize(J,nc);
1234:   /* determine the matrix preallocation information */
1235:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1236:   for (i=xs; i<xs+nx; i++) {

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

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

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

1247:       for (k=0; k<nc; k++) {
1248:         cnt = 0;
1249:         for (l=lstart; l<lend+1; l++) {
1250:           for (p=pstart; p<pend+1; p++) {
1251:             if (l || p) {
1252:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1253:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1254:               }
1255:             } else {
1256:               if (dfill) {
1257:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1258:               } else {
1259:                 for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1260:               }
1261:             }
1262:           }
1263:         }
1264:         row    = k + nc*(slot);
1265:         maxcnt = PetscMax(maxcnt,cnt);
1266:         if (removedups) {
1267:           MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1268:         } else {
1269:           MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1270:         }
1271:       }
1272:     }
1273:   }
1274:   MatSeqAIJSetPreallocation(J,0,dnz);
1275:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1276:   MatPreallocateFinalize(dnz,onz);
1277:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1279:   /*
1280:     For each node in the grid: we get the neighbors in the local (on processor ordering
1281:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1282:     PETSc ordering.
1283:   */
1284:   if (!da->prealloc_only) {
1285:     PetscCalloc1(maxcnt,&values);
1286:     for (i=xs; i<xs+nx; i++) {

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

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

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

1297:         for (k=0; k<nc; k++) {
1298:           cnt = 0;
1299:           for (l=lstart; l<lend+1; l++) {
1300:             for (p=pstart; p<pend+1; p++) {
1301:               if (l || p) {
1302:                 if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
1303:                   for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1304:                 }
1305:               } else {
1306:                 if (dfill) {
1307:                   for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1308:                 } else {
1309:                   for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1310:                 }
1311:               }
1312:             }
1313:           }
1314:           row  = k + nc*(slot);
1315:           MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1316:         }
1317:       }
1318:     }
1319:     PetscFree(values);
1320:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1321:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1322:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1323:   }
1324:   PetscFree(cols);
1325:   return(0);
1326: }

1328: /* ---------------------------------------------------------------------------------*/

1330: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1331: {
1332:   PetscErrorCode         ierr;
1333:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1334:   PetscInt               m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1335:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1336:   MPI_Comm               comm;
1337:   PetscScalar            *values;
1338:   DMBoundaryType         bx,by,bz;
1339:   ISLocalToGlobalMapping ltog;
1340:   DMDAStencilType        st;
1341:   PetscBool              removedups = PETSC_FALSE;

1344:   /*
1345:          nc - number of components per grid point
1346:          col - number of colors needed in one direction for single component problem

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

1352:   /*
1353:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1354:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1355:   */
1356:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1357:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1358:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

1364:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1365:   DMGetLocalToGlobalMapping(da,&ltog);

1367:   MatSetBlockSize(J,nc);
1368:   /* determine the matrix preallocation information */
1369:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1370:   for (i=xs; i<xs+nx; i++) {
1371:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1372:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1373:     for (j=ys; j<ys+ny; j++) {
1374:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1375:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1376:       for (k=zs; k<zs+nz; k++) {
1377:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1378:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1382:         cnt = 0;
1383:         for (l=0; l<nc; l++) {
1384:           for (ii=istart; ii<iend+1; ii++) {
1385:             for (jj=jstart; jj<jend+1; jj++) {
1386:               for (kk=kstart; kk<kend+1; kk++) {
1387:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1388:                   cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1389:                 }
1390:               }
1391:             }
1392:           }
1393:           rows[l] = l + nc*(slot);
1394:         }
1395:         if (removedups) {
1396:           MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1397:         } else {
1398:           MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1399:         }
1400:       }
1401:     }
1402:   }
1403:   MatSetBlockSize(J,nc);
1404:   MatSeqAIJSetPreallocation(J,0,dnz);
1405:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1406:   MatPreallocateFinalize(dnz,onz);
1407:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1409:   /*
1410:     For each node in the grid: we get the neighbors in the local (on processor ordering
1411:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1412:     PETSc ordering.
1413:   */
1414:   if (!da->prealloc_only) {
1415:     PetscCalloc1(col*col*col*nc*nc*nc,&values);
1416:     for (i=xs; i<xs+nx; i++) {
1417:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1418:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1419:       for (j=ys; j<ys+ny; j++) {
1420:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1421:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1422:         for (k=zs; k<zs+nz; k++) {
1423:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1424:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1428:           cnt = 0;
1429:           for (l=0; l<nc; l++) {
1430:             for (ii=istart; ii<iend+1; ii++) {
1431:               for (jj=jstart; jj<jend+1; jj++) {
1432:                 for (kk=kstart; kk<kend+1; kk++) {
1433:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1434:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1435:                   }
1436:                 }
1437:               }
1438:             }
1439:             rows[l] = l + nc*(slot);
1440:           }
1441:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1442:         }
1443:       }
1444:     }
1445:     PetscFree(values);
1446:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1447:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1448:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1449:   }
1450:   PetscFree2(rows,cols);
1451:   return(0);
1452: }

1454: /* ---------------------------------------------------------------------------------*/

1456: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1457: {
1458:   PetscErrorCode         ierr;
1459:   DM_DA                  *dd = (DM_DA*)da->data;
1460:   PetscInt               xs,nx,i,j,gxs,gnx,row,k,l;
1461:   PetscInt               m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1462:   PetscInt               *ofill = dd->ofill,*dfill = dd->dfill;
1463:   PetscScalar            *values;
1464:   DMBoundaryType         bx;
1465:   ISLocalToGlobalMapping ltog;
1466:   PetscMPIInt            rank,size;

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

1473:   /*
1474:          nc - number of components per grid point

1476:   */
1477:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1478:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1479:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1481:   MatSetBlockSize(J,nc);
1482:   PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);

1484:   /*
1485:         note should be smaller for first and last process with no periodic
1486:         does not handle dfill
1487:   */
1488:   cnt = 0;
1489:   /* coupling with process to the left */
1490:   for (i=0; i<s; i++) {
1491:     for (j=0; j<nc; j++) {
1492:       ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1493:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1494:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1495:       cnt++;
1496:     }
1497:   }
1498:   for (i=s; i<nx-s; i++) {
1499:     for (j=0; j<nc; j++) {
1500:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1501:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1502:       cnt++;
1503:     }
1504:   }
1505:   /* coupling with process to the right */
1506:   for (i=nx-s; i<nx; i++) {
1507:     for (j=0; j<nc; j++) {
1508:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1509:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1510:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1511:       cnt++;
1512:     }
1513:   }

1515:   MatSeqAIJSetPreallocation(J,0,cols);
1516:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1517:   PetscFree2(cols,ocols);

1519:   DMGetLocalToGlobalMapping(da,&ltog);
1520:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1522:   /*
1523:     For each node in the grid: we get the neighbors in the local (on processor ordering
1524:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1525:     PETSc ordering.
1526:   */
1527:   if (!da->prealloc_only) {
1528:     PetscCalloc2(maxcnt,&values,maxcnt,&cols);

1530:     row = xs*nc;
1531:     /* coupling with process to the left */
1532:     for (i=xs; i<xs+s; i++) {
1533:       for (j=0; j<nc; j++) {
1534:         cnt = 0;
1535:         if (rank) {
1536:           for (l=0; l<s; l++) {
1537:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1538:           }
1539:         }
1540:         if (dfill) {
1541:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1542:             cols[cnt++] = i*nc + dfill[k];
1543:           }
1544:         } else {
1545:           for (k=0; k<nc; k++) {
1546:             cols[cnt++] = i*nc + k;
1547:           }
1548:         }
1549:         for (l=0; l<s; l++) {
1550:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1551:         }
1552:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1553:         row++;
1554:       }
1555:     }
1556:     for (i=xs+s; i<xs+nx-s; i++) {
1557:       for (j=0; j<nc; j++) {
1558:         cnt = 0;
1559:         for (l=0; l<s; l++) {
1560:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1561:         }
1562:         if (dfill) {
1563:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1564:             cols[cnt++] = i*nc + dfill[k];
1565:           }
1566:         } else {
1567:           for (k=0; k<nc; k++) {
1568:             cols[cnt++] = i*nc + k;
1569:           }
1570:         }
1571:         for (l=0; l<s; l++) {
1572:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1573:         }
1574:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1575:         row++;
1576:       }
1577:     }
1578:     /* coupling with process to the right */
1579:     for (i=xs+nx-s; i<xs+nx; i++) {
1580:       for (j=0; j<nc; j++) {
1581:         cnt = 0;
1582:         for (l=0; l<s; l++) {
1583:           for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1584:         }
1585:         if (dfill) {
1586:           for (k=dfill[j]; k<dfill[j+1]; k++) {
1587:             cols[cnt++] = i*nc + dfill[k];
1588:           }
1589:         } else {
1590:           for (k=0; k<nc; k++) {
1591:             cols[cnt++] = i*nc + k;
1592:           }
1593:         }
1594:         if (rank < size-1) {
1595:           for (l=0; l<s; l++) {
1596:             for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1597:           }
1598:         }
1599:         MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1600:         row++;
1601:       }
1602:     }
1603:     PetscFree2(values,cols);
1604:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1605:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1606:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1607:   }
1608:   return(0);
1609: }

1611: /* ---------------------------------------------------------------------------------*/

1613: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1614: {
1615:   PetscErrorCode         ierr;
1616:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1617:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1618:   PetscInt               istart,iend;
1619:   PetscScalar            *values;
1620:   DMBoundaryType         bx;
1621:   ISLocalToGlobalMapping ltog;

1624:   /*
1625:          nc - number of components per grid point
1626:          col - number of colors needed in one direction for single component problem

1628:   */
1629:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1630:   col  = 2*s + 1;

1632:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1633:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1635:   MatSetBlockSize(J,nc);
1636:   MatSeqAIJSetPreallocation(J,col*nc,0);
1637:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1639:   DMGetLocalToGlobalMapping(da,&ltog);
1640:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1642:   /*
1643:     For each node in the grid: we get the neighbors in the local (on processor ordering
1644:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1645:     PETSc ordering.
1646:   */
1647:   if (!da->prealloc_only) {
1648:     PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1649:     PetscCalloc1(col*nc*nc,&values);
1650:     for (i=xs; i<xs+nx; i++) {
1651:       istart = PetscMax(-s,gxs - i);
1652:       iend   = PetscMin(s,gxs + gnx - i - 1);
1653:       slot   = i - gxs;

1655:       cnt = 0;
1656:       for (l=0; l<nc; l++) {
1657:         for (i1=istart; i1<iend+1; i1++) {
1658:           cols[cnt++] = l + nc*(slot + i1);
1659:         }
1660:         rows[l] = l + nc*(slot);
1661:       }
1662:       MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1663:     }
1664:     PetscFree(values);
1665:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1666:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1667:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1668:     PetscFree2(rows,cols);
1669:   }
1670:   return(0);
1671: }

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

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

1693:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1694:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1695:   PetscObjectGetComm((PetscObject)da,&comm);

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

1699:   DMGetLocalToGlobalMapping(da,&ltog);

1701:   /* determine the matrix preallocation information */
1702:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1703:   for (i=xs; i<xs+nx; i++) {
1704:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1705:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1706:     for (j=ys; j<ys+ny; j++) {
1707:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1708:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1709:       slot   = i - gxs + gnx*(j - gys);

1711:       /* Find block columns in block row */
1712:       cnt = 0;
1713:       for (ii=istart; ii<iend+1; ii++) {
1714:         for (jj=jstart; jj<jend+1; jj++) {
1715:           if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1716:             cols[cnt++] = slot + ii + gnx*jj;
1717:           }
1718:         }
1719:       }
1720:       MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1721:     }
1722:   }
1723:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1724:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1725:   MatPreallocateFinalize(dnz,onz);

1727:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1729:   /*
1730:     For each node in the grid: we get the neighbors in the local (on processor ordering
1731:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1732:     PETSc ordering.
1733:   */
1734:   if (!da->prealloc_only) {
1735:     PetscCalloc1(col*col*nc*nc,&values);
1736:     for (i=xs; i<xs+nx; i++) {
1737:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1738:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1739:       for (j=ys; j<ys+ny; j++) {
1740:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1741:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1742:         slot = i - gxs + gnx*(j - gys);
1743:         cnt  = 0;
1744:         for (ii=istart; ii<iend+1; ii++) {
1745:           for (jj=jstart; jj<jend+1; jj++) {
1746:             if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1747:               cols[cnt++] = slot + ii + gnx*jj;
1748:             }
1749:           }
1750:         }
1751:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1752:       }
1753:     }
1754:     PetscFree(values);
1755:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1756:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1757:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1758:   }
1759:   PetscFree(cols);
1760:   return(0);
1761: }

1763: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1764: {
1765:   PetscErrorCode         ierr;
1766:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1767:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1768:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1769:   MPI_Comm               comm;
1770:   PetscScalar            *values;
1771:   DMBoundaryType         bx,by,bz;
1772:   DMDAStencilType        st;
1773:   ISLocalToGlobalMapping ltog;

1776:   /*
1777:          nc - number of components per grid point
1778:          col - number of colors needed in one direction for single component problem

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

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

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

1790:   DMGetLocalToGlobalMapping(da,&ltog);

1792:   /* determine the matrix preallocation information */
1793:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1794:   for (i=xs; i<xs+nx; i++) {
1795:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1796:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1797:     for (j=ys; j<ys+ny; j++) {
1798:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1799:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1800:       for (k=zs; k<zs+nz; k++) {
1801:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1802:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1806:         /* Find block columns in block row */
1807:         cnt = 0;
1808:         for (ii=istart; ii<iend+1; ii++) {
1809:           for (jj=jstart; jj<jend+1; jj++) {
1810:             for (kk=kstart; kk<kend+1; kk++) {
1811:               if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1812:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1813:               }
1814:             }
1815:           }
1816:         }
1817:         MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1818:       }
1819:     }
1820:   }
1821:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1822:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1823:   MatPreallocateFinalize(dnz,onz);

1825:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1827:   /*
1828:     For each node in the grid: we get the neighbors in the local (on processor ordering
1829:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1830:     PETSc ordering.
1831:   */
1832:   if (!da->prealloc_only) {
1833:     PetscCalloc1(col*col*col*nc*nc,&values);
1834:     for (i=xs; i<xs+nx; i++) {
1835:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1836:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1837:       for (j=ys; j<ys+ny; j++) {
1838:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1839:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1840:         for (k=zs; k<zs+nz; k++) {
1841:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1842:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

1846:           cnt = 0;
1847:           for (ii=istart; ii<iend+1; ii++) {
1848:             for (jj=jstart; jj<jend+1; jj++) {
1849:               for (kk=kstart; kk<kend+1; kk++) {
1850:                 if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1851:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1852:                 }
1853:               }
1854:             }
1855:           }
1856:           MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1857:         }
1858:       }
1859:     }
1860:     PetscFree(values);
1861:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1862:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1863:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1864:   }
1865:   PetscFree(cols);
1866:   return(0);
1867: }

1869: /*
1870:   This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1871:   identify in the local ordering with periodic domain.
1872: */
1873: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1874: {
1876:   PetscInt       i,n;

1879:   ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1880:   ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1881:   for (i=0,n=0; i<*cnt; i++) {
1882:     if (col[i] >= *row) col[n++] = col[i];
1883:   }
1884:   *cnt = n;
1885:   return(0);
1886: }

1888: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1889: {
1890:   PetscErrorCode         ierr;
1891:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1892:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1893:   PetscInt               istart,iend,jstart,jend,ii,jj;
1894:   MPI_Comm               comm;
1895:   PetscScalar            *values;
1896:   DMBoundaryType         bx,by;
1897:   DMDAStencilType        st;
1898:   ISLocalToGlobalMapping ltog;

1901:   /*
1902:      nc - number of components per grid point
1903:      col - number of colors needed in one direction for single component problem
1904:   */
1905:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1906:   col  = 2*s + 1;

1908:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1909:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1910:   PetscObjectGetComm((PetscObject)da,&comm);

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

1914:   DMGetLocalToGlobalMapping(da,&ltog);

1916:   /* determine the matrix preallocation information */
1917:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1918:   for (i=xs; i<xs+nx; i++) {
1919:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1920:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
1921:     for (j=ys; j<ys+ny; j++) {
1922:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1923:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
1924:       slot   = i - gxs + gnx*(j - gys);

1926:       /* Find block columns in block row */
1927:       cnt = 0;
1928:       for (ii=istart; ii<iend+1; ii++) {
1929:         for (jj=jstart; jj<jend+1; jj++) {
1930:           if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1931:             cols[cnt++] = slot + ii + gnx*jj;
1932:           }
1933:         }
1934:       }
1935:       L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1936:       MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1937:     }
1938:   }
1939:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1940:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1941:   MatPreallocateFinalize(dnz,onz);

1943:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

1960:         /* Find block columns in block row */
1961:         cnt = 0;
1962:         for (ii=istart; ii<iend+1; ii++) {
1963:           for (jj=jstart; jj<jend+1; jj++) {
1964:             if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1965:               cols[cnt++] = slot + ii + gnx*jj;
1966:             }
1967:           }
1968:         }
1969:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1970:         MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1971:       }
1972:     }
1973:     PetscFree(values);
1974:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1975:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1976:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1977:   }
1978:   PetscFree(cols);
1979:   return(0);
1980: }

1982: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1983: {
1984:   PetscErrorCode         ierr;
1985:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1986:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1987:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1988:   MPI_Comm               comm;
1989:   PetscScalar            *values;
1990:   DMBoundaryType         bx,by,bz;
1991:   DMDAStencilType        st;
1992:   ISLocalToGlobalMapping ltog;

1995:   /*
1996:      nc - number of components per grid point
1997:      col - number of colors needed in one direction for single component problem
1998:   */
1999:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2000:   col  = 2*s + 1;

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

2006:   /* create the matrix */
2007:   PetscMalloc1(col*col*col,&cols);

2009:   DMGetLocalToGlobalMapping(da,&ltog);

2011:   /* determine the matrix preallocation information */
2012:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2013:   for (i=xs; i<xs+nx; i++) {
2014:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2015:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2016:     for (j=ys; j<ys+ny; j++) {
2017:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2018:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2019:       for (k=zs; k<zs+nz; k++) {
2020:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2021:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2025:         /* Find block columns in block row */
2026:         cnt = 0;
2027:         for (ii=istart; ii<iend+1; ii++) {
2028:           for (jj=jstart; jj<jend+1; jj++) {
2029:             for (kk=kstart; kk<kend+1; kk++) {
2030:               if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2031:                 cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2032:               }
2033:             }
2034:           }
2035:         }
2036:         L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2037:         MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2038:       }
2039:     }
2040:   }
2041:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2042:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2043:   MatPreallocateFinalize(dnz,onz);

2045:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2047:   /*
2048:     For each node in the grid: we get the neighbors in the local (on processor ordering
2049:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2050:     PETSc ordering.
2051:   */
2052:   if (!da->prealloc_only) {
2053:     PetscCalloc1(col*col*col*nc*nc,&values);
2054:     for (i=xs; i<xs+nx; i++) {
2055:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2056:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2057:       for (j=ys; j<ys+ny; j++) {
2058:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2059:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2060:         for (k=zs; k<zs+nz; k++) {
2061:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2062:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2066:           cnt = 0;
2067:           for (ii=istart; ii<iend+1; ii++) {
2068:             for (jj=jstart; jj<jend+1; jj++) {
2069:               for (kk=kstart; kk<kend+1; kk++) {
2070:                 if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2071:                   cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2072:                 }
2073:               }
2074:             }
2075:           }
2076:           L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2077:           MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2078:         }
2079:       }
2080:     }
2081:     PetscFree(values);
2082:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2083:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2084:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2085:   }
2086:   PetscFree(cols);
2087:   return(0);
2088: }

2090: /* ---------------------------------------------------------------------------------*/

2092: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2093: {
2094:   PetscErrorCode         ierr;
2095:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2096:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2097:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2098:   DM_DA                  *dd = (DM_DA*)da->data;
2099:   PetscInt               ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2100:   MPI_Comm               comm;
2101:   PetscScalar            *values;
2102:   DMBoundaryType         bx,by,bz;
2103:   ISLocalToGlobalMapping ltog;
2104:   DMDAStencilType        st;
2105:   PetscBool              removedups = PETSC_FALSE;

2108:   /*
2109:          nc - number of components per grid point
2110:          col - number of colors needed in one direction for single component problem

2112:   */
2113:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2114:   col  = 2*s + 1;
2115:   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\
2116:                  by 2*stencil_width + 1\n");
2117:   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\
2118:                  by 2*stencil_width + 1\n");
2119:   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\
2120:                  by 2*stencil_width + 1\n");

2122:   /*
2123:        With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2124:        because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2125:   */
2126:   if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2127:   if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2128:   if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;

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

2134:   PetscMalloc1(col*col*col*nc,&cols);
2135:   DMGetLocalToGlobalMapping(da,&ltog);

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

2140:   MatSetBlockSize(J,nc);
2141:   for (i=xs; i<xs+nx; i++) {
2142:     istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2143:     iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2144:     for (j=ys; j<ys+ny; j++) {
2145:       jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2146:       jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2147:       for (k=zs; k<zs+nz; k++) {
2148:         kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2149:         kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2153:         for (l=0; l<nc; l++) {
2154:           cnt = 0;
2155:           for (ii=istart; ii<iend+1; ii++) {
2156:             for (jj=jstart; jj<jend+1; jj++) {
2157:               for (kk=kstart; kk<kend+1; kk++) {
2158:                 if (ii || jj || kk) {
2159:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2160:                     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);
2161:                   }
2162:                 } else {
2163:                   if (dfill) {
2164:                     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);
2165:                   } else {
2166:                     for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2167:                   }
2168:                 }
2169:               }
2170:             }
2171:           }
2172:           row  = l + nc*(slot);
2173:           maxcnt = PetscMax(maxcnt,cnt);
2174:           if (removedups) {
2175:             MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2176:           } else {
2177:             MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2178:           }
2179:         }
2180:       }
2181:     }
2182:   }
2183:   MatSeqAIJSetPreallocation(J,0,dnz);
2184:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2185:   MatPreallocateFinalize(dnz,onz);
2186:   MatSetLocalToGlobalMapping(J,ltog,ltog);

2188:   /*
2189:     For each node in the grid: we get the neighbors in the local (on processor ordering
2190:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2191:     PETSc ordering.
2192:   */
2193:   if (!da->prealloc_only) {
2194:     PetscCalloc1(maxcnt,&values);
2195:     for (i=xs; i<xs+nx; i++) {
2196:       istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2197:       iend   = (bx == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,m-i-1));
2198:       for (j=ys; j<ys+ny; j++) {
2199:         jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2200:         jend   = (by == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,n-j-1));
2201:         for (k=zs; k<zs+nz; k++) {
2202:           kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2203:           kend   = (bz == DM_BOUNDARY_PERIODIC) ?  s : (PetscMin(s,p-k-1));

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

2207:           for (l=0; l<nc; l++) {
2208:             cnt = 0;
2209:             for (ii=istart; ii<iend+1; ii++) {
2210:               for (jj=jstart; jj<jend+1; jj++) {
2211:                 for (kk=kstart; kk<kend+1; kk++) {
2212:                   if (ii || jj || kk) {
2213:                     if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2214:                       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);
2215:                     }
2216:                   } else {
2217:                     if (dfill) {
2218:                       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);
2219:                     } else {
2220:                       for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2221:                     }
2222:                   }
2223:                 }
2224:               }
2225:             }
2226:             row  = l + nc*(slot);
2227:             MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2228:           }
2229:         }
2230:       }
2231:     }
2232:     PetscFree(values);
2233:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2234:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2235:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2236:   }
2237:   PetscFree(cols);
2238:   return(0);
2239: }