Actual source code: fdda.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: #include <petsc/private/dmdaimpl.h> /*I      "petscdmda.h"     I*/
  3: #include <petscmat.h>
  4: #include <petsc/private/matimpl.h>

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

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

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

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

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

 50:   *rfill = fill;
 51:   return(0);
 52: }

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

 60:     Logically Collective on DMDA

 62:     Input Parameter:
 63: +   da - the distributed array
 64: .   dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
 65: -   ofill - the fill pattern in the off-diagonal blocks


 68:     Level: developer

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

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

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

 85:    Contributed by Glenn Hammond

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

 89: @*/
 90: PetscErrorCode  DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
 91: {
 92:   DM_DA          *dd = (DM_DA*)da->data;
 94:   PetscInt       i,k,cnt = 1;

 97:   DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
 98:   DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);

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


117: PetscErrorCode  DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
118: {
119:   PetscErrorCode   ierr;
120:   PetscInt         dim,m,n,p,nc;
121:   DMBoundaryType bx,by,bz;
122:   MPI_Comm         comm;
123:   PetscMPIInt      size;
124:   PetscBool        isBAIJ;
125:   DM_DA            *dd = (DM_DA*)da->data;

128:   /*
129:                                   m
130:           ------------------------------------------------------
131:          |                                                     |
132:          |                                                     |
133:          |               ----------------------                |
134:          |               |                    |                |
135:       n  |           yn  |                    |                |
136:          |               |                    |                |
137:          |               .---------------------                |
138:          |             (xs,ys)     xn                          |
139:          |            .                                        |
140:          |         (gxs,gys)                                   |
141:          |                                                     |
142:           -----------------------------------------------------
143:   */

145:   /*
146:          nc - number of components per grid point
147:          col - number of colors needed in one direction for single component problem

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

152:   PetscObjectGetComm((PetscObject)da,&comm);
153:   MPI_Comm_size(comm,&size);
154:   if (ctype == IS_COLORING_GHOSTED) {
155:     if (size == 1) {
156:       ctype = IS_COLORING_GLOBAL;
157:     } else if (dim > 1) {
158:       if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
159:         SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_GHOSTED cannot be used for periodic boundary condition having both ends of the domain  on the same process");
160:       }
161:     }
162:   }

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

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

199: /* ---------------------------------------------------------------------------------*/

203: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
204: {
205:   PetscErrorCode   ierr;
206:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
207:   PetscInt         ncolors;
208:   MPI_Comm         comm;
209:   DMBoundaryType bx,by;
210:   DMDAStencilType  st;
211:   ISColoringValue  *colors;
212:   DM_DA            *dd = (DM_DA*)da->data;

215:   /*
216:          nc - number of components per grid point
217:          col - number of colors needed in one direction for single component problem

219:   */
220:   DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
221:   col  = 2*s + 1;
222:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
223:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
224:   PetscObjectGetComm((PetscObject)da,&comm);

226:   /* special case as taught to us by Paul Hovland */
227:   if (st == DMDA_STENCIL_STAR && s == 1) {
228:     DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
229:   } else {

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

266:         ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
267:       }
268:       *coloring = dd->ghostedcoloring;
269:     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
270:   }
271:   ISColoringReference(*coloring);
272:   return(0);
273: }

275: /* ---------------------------------------------------------------------------------*/

279: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
280: {
281:   PetscErrorCode   ierr;
282:   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;
283:   PetscInt         ncolors;
284:   MPI_Comm         comm;
285:   DMBoundaryType bx,by,bz;
286:   DMDAStencilType  st;
287:   ISColoringValue  *colors;
288:   DM_DA            *dd = (DM_DA*)da->data;

291:   /*
292:          nc - number of components per grid point
293:          col - number of colors needed in one direction for single component problem

295:   */
296:   DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
297:   col  = 2*s + 1;
298:   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\
299:                                                          by 2*stencil_width + 1\n");
300:   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\
301:                                                          by 2*stencil_width + 1\n");
302:   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\
303:                                                          by 2*stencil_width + 1\n");

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

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

351: /* ---------------------------------------------------------------------------------*/

355: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
356: {
357:   PetscErrorCode   ierr;
358:   PetscInt         xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
359:   PetscInt         ncolors;
360:   MPI_Comm         comm;
361:   DMBoundaryType bx;
362:   ISColoringValue  *colors;
363:   DM_DA            *dd = (DM_DA*)da->data;

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

370:   */
371:   DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
372:   col  = 2*s + 1;

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

377:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
378:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
379:   PetscObjectGetComm((PetscObject)da,&comm);

381:   /* create the coloring */
382:   if (ctype == IS_COLORING_GLOBAL) {
383:     if (!dd->localcoloring) {
384:       PetscMalloc1(nc*nx,&colors);
385:       if (dd->ofillcols) {
386:         PetscInt tc = 0;
387:         for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
388:         i1 = 0;
389:         for (i=xs; i<xs+nx; i++) {
390:           for (l=0; l<nc; l++) {
391:             if (dd->ofillcols[l] && (i % col)) {
392:               colors[i1++] =  nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
393:             } else {
394:               colors[i1++] = l;
395:             }
396:           }
397:         }
398:         ncolors = nc + 2*s*tc;
399:       } else {
400:         i1 = 0;
401:         for (i=xs; i<xs+nx; i++) {
402:           for (l=0; l<nc; l++) {
403:             colors[i1++] = l + nc*(i % col);
404:           }
405:         }
406:         ncolors = nc + nc*(col-1);
407:       }
408:       ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
409:     }
410:     *coloring = dd->localcoloring;
411:   } else if (ctype == IS_COLORING_GHOSTED) {
412:     if (!dd->ghostedcoloring) {
413:       PetscMalloc1(nc*gnx,&colors);
414:       i1   = 0;
415:       for (i=gxs; i<gxs+gnx; i++) {
416:         for (l=0; l<nc; l++) {
417:           /* the complicated stuff is to handle periodic boundaries */
418:           colors[i1++] = l + nc*(SetInRange(i,m) % col);
419:         }
420:       }
421:       ncolors = nc + nc*(col-1);
422:       ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
423:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
424:     }
425:     *coloring = dd->ghostedcoloring;
426:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
427:   ISColoringReference(*coloring);
428:   return(0);
429: }

433: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
434: {
435:   PetscErrorCode   ierr;
436:   PetscInt         xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
437:   PetscInt         ncolors;
438:   MPI_Comm         comm;
439:   DMBoundaryType bx,by;
440:   ISColoringValue  *colors;
441:   DM_DA            *dd = (DM_DA*)da->data;

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

448:   */
449:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
450:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
451:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
452:   PetscObjectGetComm((PetscObject)da,&comm);

454:   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");
455:   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");

457:   /* create the coloring */
458:   if (ctype == IS_COLORING_GLOBAL) {
459:     if (!dd->localcoloring) {
460:       PetscMalloc1(nc*nx*ny,&colors);
461:       ii   = 0;
462:       for (j=ys; j<ys+ny; j++) {
463:         for (i=xs; i<xs+nx; i++) {
464:           for (k=0; k<nc; k++) {
465:             colors[ii++] = k + nc*((3*j+i) % 5);
466:           }
467:         }
468:       }
469:       ncolors = 5*nc;
470:       ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
471:     }
472:     *coloring = dd->localcoloring;
473:   } else if (ctype == IS_COLORING_GHOSTED) {
474:     if (!dd->ghostedcoloring) {
475:       PetscMalloc1(nc*gnx*gny,&colors);
476:       ii = 0;
477:       for (j=gys; j<gys+gny; j++) {
478:         for (i=gxs; i<gxs+gnx; i++) {
479:           for (k=0; k<nc; k++) {
480:             colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
481:           }
482:         }
483:       }
484:       ncolors = 5*nc;
485:       ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
486:       ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
487:     }
488:     *coloring = dd->ghostedcoloring;
489:   } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
490:   return(0);
491: }

493: /* =========================================================================== */
494: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
495: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
496: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
497: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
498: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
499: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
500: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
501: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
502: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
503: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);

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

510:    Logically Collective on Mat

512:    Input Parameters:
513: +  mat - the matrix
514: -  da - the da

516:    Level: intermediate

518: @*/
519: PetscErrorCode MatSetupDM(Mat mat,DM da)
520: {

526:   PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
527:   return(0);
528: }

532: PetscErrorCode  MatView_MPI_DA(Mat A,PetscViewer viewer)
533: {
534:   DM                da;
535:   PetscErrorCode    ierr;
536:   const char        *prefix;
537:   Mat               Anatural;
538:   AO                ao;
539:   PetscInt          rstart,rend,*petsc,i;
540:   IS                is;
541:   MPI_Comm          comm;
542:   PetscViewerFormat format;

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

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

553:   DMDAGetAO(da,&ao);
554:   MatGetOwnershipRange(A,&rstart,&rend);
555:   PetscMalloc1(rend-rstart,&petsc);
556:   for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
557:   AOApplicationToPetsc(ao,rend-rstart,petsc);
558:   ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);

560:   /* call viewer on natural ordering */
561:   MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
562:   ISDestroy(&is);
563:   PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
564:   PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
565:   PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
566:   (Anatural->ops->view)(Anatural,viewer);
567:   MatDestroy(&Anatural);
568:   return(0);
569: }

573: PetscErrorCode  MatLoad_MPI_DA(Mat A,PetscViewer viewer)
574: {
575:   DM             da;
577:   Mat            Anatural,Aapp;
578:   AO             ao;
579:   PetscInt       rstart,rend,*app,i;
580:   IS             is;
581:   MPI_Comm       comm;

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

588:   /* Load the matrix in natural ordering */
589:   MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
590:   MatSetType(Anatural,((PetscObject)A)->type_name);
591:   MatSetSizes(Anatural,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
592:   MatLoad(Anatural,viewer);

594:   /* Map natural ordering to application ordering and create IS */
595:   DMDAGetAO(da,&ao);
596:   MatGetOwnershipRange(Anatural,&rstart,&rend);
597:   PetscMalloc1(rend-rstart,&app);
598:   for (i=rstart; i<rend; i++) app[i-rstart] = i;
599:   AOPetscToApplication(ao,rend-rstart,app);
600:   ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);

602:   /* Do permutation and replace header */
603:   MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
604:   MatHeaderReplace(A,&Aapp);
605:   ISDestroy(&is);
606:   MatDestroy(&Anatural);
607:   return(0);
608: }

612: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
613: {
615:   PetscInt       dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
616:   Mat            A;
617:   MPI_Comm       comm;
618:   MatType        Atype;
619:   PetscSection   section, sectionGlobal;
620:   void           (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL;
621:   MatType        mtype;
622:   PetscMPIInt    size;
623:   DM_DA          *dd = (DM_DA*)da->data;

626:   MatInitializePackage();
627:   mtype = da->mattype;

629:   DMGetDefaultSection(da, &section);
630:   if (section) {
631:     PetscInt  bs = -1;
632:     PetscInt  localSize;
633:     PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;

635:     DMGetDefaultGlobalSection(da, &sectionGlobal);
636:     PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
637:     MatCreate(PetscObjectComm((PetscObject)da), J);
638:     MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
639:     MatSetType(*J, mtype);
640:     MatSetFromOptions(*J);
641:     PetscStrcmp(mtype, MATSHELL, &isShell);
642:     PetscStrcmp(mtype, MATBAIJ, &isBlock);
643:     PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
644:     PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
645:     PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
646:     PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
647:     PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
648:     /* Check for symmetric storage */
649:     isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
650:     if (isSymmetric) {
651:       MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
652:     }
653:     if (!isShell) {
654:       PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;

656:       if (bs < 0) {
657:         if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
658:           PetscInt pStart, pEnd, p, dof;

660:           PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
661:           for (p = pStart; p < pEnd; ++p) {
662:             PetscSectionGetDof(sectionGlobal, p, &dof);
663:             if (dof) {
664:               bs = dof;
665:               break;
666:             }
667:           }
668:         } else {
669:           bs = 1;
670:         }
671:         /* Must have same blocksize on all procs (some might have no points) */
672:         bsLocal = bs;
673:         MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
674:       }
675:       PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
676:       /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
677:       PetscFree4(dnz, onz, dnzu, onzu);
678:     }
679:   }
680:   /*
681:                                   m
682:           ------------------------------------------------------
683:          |                                                     |
684:          |                                                     |
685:          |               ----------------------                |
686:          |               |                    |                |
687:       n  |           ny  |                    |                |
688:          |               |                    |                |
689:          |               .---------------------                |
690:          |             (xs,ys)     nx                          |
691:          |            .                                        |
692:          |         (gxs,gys)                                   |
693:          |                                                     |
694:           -----------------------------------------------------
695:   */

697:   /*
698:          nc - number of components per grid point
699:          col - number of colors needed in one direction for single component problem

701:   */
702:   M   = dd->M;
703:   N   = dd->N;
704:   P   = dd->P;
705:   dim = da->dim;
706:   dof = dd->w;
707:   /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
708:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
709:   PetscObjectGetComm((PetscObject)da,&comm);
710:   MatCreate(comm,&A);
711:   MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
712:   MatSetType(A,mtype);
713:   MatSetDM(A,da);
714:   MatSetFromOptions(A);
715:   MatGetType(A,&Atype);
716:   /*
717:      We do not provide a getmatrix function in the DMDA operations because
718:    the basic DMDA does not know about matrices. We think of DMDA as being more
719:    more low-level than matrices. This is kind of cheating but, cause sometimes
720:    we think of DMDA has higher level than matrices.

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

793: /* ---------------------------------------------------------------------------------*/
796: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
797: {
798:   PetscErrorCode         ierr;
799:   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;
800:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
801:   MPI_Comm               comm;
802:   PetscScalar            *values;
803:   DMBoundaryType         bx,by;
804:   ISLocalToGlobalMapping ltog;
805:   DMDAStencilType        st;

808:   /*
809:          nc - number of components per grid point
810:          col - number of colors needed in one direction for single component problem

812:   */
813:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
814:   col  = 2*s + 1;
815:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
816:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
817:   PetscObjectGetComm((PetscObject)da,&comm);

819:   PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
820:   DMGetLocalToGlobalMapping(da,&ltog);

822:   MatSetBlockSize(J,nc);
823:   /* determine the matrix preallocation information */
824:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
825:   for (i=xs; i<xs+nx; i++) {

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

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

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

836:       cnt = 0;
837:       for (k=0; k<nc; k++) {
838:         for (l=lstart; l<lend+1; l++) {
839:           for (p=pstart; p<pend+1; p++) {
840:             if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
841:               cols[cnt++] = k + nc*(slot + gnx*l + p);
842:             }
843:           }
844:         }
845:         rows[k] = k + nc*(slot);
846:       }
847:       MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
848:     }
849:   }
850:   MatSetBlockSize(J,nc);
851:   MatSeqAIJSetPreallocation(J,0,dnz);
852:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
853:   MatPreallocateFinalize(dnz,onz);

855:   MatSetLocalToGlobalMapping(J,ltog,ltog);

857:   /*
858:     For each node in the grid: we get the neighbors in the local (on processor ordering
859:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
860:     PETSc ordering.
861:   */
862:   if (!da->prealloc_only) {
863:     PetscCalloc1(col*col*nc*nc,&values);
864:     for (i=xs; i<xs+nx; i++) {

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

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

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

875:         cnt = 0;
876:         for (k=0; k<nc; k++) {
877:           for (l=lstart; l<lend+1; l++) {
878:             for (p=pstart; p<pend+1; p++) {
879:               if ((st == DMDA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
880:                 cols[cnt++] = k + nc*(slot + gnx*l + p);
881:               }
882:             }
883:           }
884:           rows[k] = k + nc*(slot);
885:         }
886:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
887:       }
888:     }
889:     PetscFree(values);
890:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
891:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
892:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
893:   }
894:   PetscFree2(rows,cols);
895:   return(0);
896: }

900: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
901: {
902:   PetscErrorCode         ierr;
903:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
904:   PetscInt               m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p;
905:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
906:   DM_DA                  *dd = (DM_DA*)da->data;
907:   PetscInt               ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
908:   MPI_Comm               comm;
909:   PetscScalar            *values;
910:   DMBoundaryType         bx,by;
911:   ISLocalToGlobalMapping ltog;
912:   DMDAStencilType        st;

915:   /*
916:          nc - number of components per grid point
917:          col - number of colors needed in one direction for single component problem

919:   */
920:   DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
921:   col  = 2*s + 1;
922:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
923:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
924:   PetscObjectGetComm((PetscObject)da,&comm);

926:   PetscMalloc1(col*col*nc,&cols);
927:   DMGetLocalToGlobalMapping(da,&ltog);

929:   MatSetBlockSize(J,nc);
930:   /* determine the matrix preallocation information */
931:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
932:   for (i=xs; i<xs+nx; i++) {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1049:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1050:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

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

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

1109:           cnt = 0;
1110:           for (l=0; l<nc; l++) {
1111:             for (ii=istart; ii<iend+1; ii++) {
1112:               for (jj=jstart; jj<jend+1; jj++) {
1113:                 for (kk=kstart; kk<kend+1; kk++) {
1114:                   if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1115:                     cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1116:                   }
1117:                 }
1118:               }
1119:             }
1120:             rows[l] = l + nc*(slot);
1121:           }
1122:           MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1123:         }
1124:       }
1125:     }
1126:     PetscFree(values);
1127:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1128:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1129:     MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1130:   }
1131:   PetscFree2(rows,cols);
1132:   return(0);
1133: }

1135: /* ---------------------------------------------------------------------------------*/

1139: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1140: {
1141:   PetscErrorCode         ierr;
1142:   DM_DA                  *dd = (DM_DA*)da->data;
1143:   PetscInt               xs,nx,i,j,gxs,gnx,row,k,l;
1144:   PetscInt               m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1145:   PetscInt               *ofill = dd->ofill,*dfill = dd->dfill;
1146:   PetscScalar            *values;
1147:   DMBoundaryType         bx;
1148:   ISLocalToGlobalMapping ltog;
1149:   PetscMPIInt            rank,size;

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

1156:   /*
1157:          nc - number of components per grid point

1159:   */
1160:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1161:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1162:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1164:   MatSetBlockSize(J,nc);
1165:   PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);

1167:   if (nx < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Need at least two grid points per process");
1168:   /*
1169:         note should be smaller for first and last process with no periodic
1170:         does not handle dfill
1171:   */
1172:   cnt = 0;
1173:   /* coupling with process to the left */
1174:   for (i=0; i<s; i++) {
1175:     for (j=0; j<nc; j++) {
1176:       ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1177:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1178:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1179:       cnt++;
1180:     }
1181:   }
1182:   for (i=s; i<nx-s; i++) {
1183:     for (j=0; j<nc; j++) {
1184:       cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1185:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1186:       cnt++;
1187:     }
1188:   }
1189:   /* coupling with process to the right */
1190:   for (i=nx-s; i<nx; i++) {
1191:     for (j=0; j<nc; j++) {
1192:       ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1193:       cols[cnt]  = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1194:       maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1195:       cnt++;
1196:     }
1197:   }

1199:   MatSeqAIJSetPreallocation(J,0,cols);
1200:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1201:   PetscFree2(cols,ocols);

1203:   DMGetLocalToGlobalMapping(da,&ltog);
1204:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1206:   /*
1207:     For each node in the grid: we get the neighbors in the local (on processor ordering
1208:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1209:     PETSc ordering.
1210:   */
1211:   if (!da->prealloc_only) {
1212:     PetscCalloc2(maxcnt,&values,maxcnt,&cols);

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

1295: /* ---------------------------------------------------------------------------------*/

1299: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1300: {
1301:   PetscErrorCode         ierr;
1302:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
1303:   PetscInt               m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1304:   PetscInt               istart,iend;
1305:   PetscScalar            *values;
1306:   DMBoundaryType         bx;
1307:   ISLocalToGlobalMapping ltog;

1310:   /*
1311:          nc - number of components per grid point
1312:          col - number of colors needed in one direction for single component problem

1314:   */
1315:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1316:   col  = 2*s + 1;

1318:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1319:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1321:   MatSetBlockSize(J,nc);
1322:   MatSeqAIJSetPreallocation(J,col*nc,0);
1323:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1325:   DMGetLocalToGlobalMapping(da,&ltog);
1326:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

1381:   DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1382:   DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1383:   PetscObjectGetComm((PetscObject)da,&comm);

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

1387:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

1415:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

1480:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

1515:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

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

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

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

1608:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

1637:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

1702:   /* create the matrix */
1703:   PetscMalloc1(col*col*col,&cols);

1705:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

1741:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

1786: /* ---------------------------------------------------------------------------------*/

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

1805:   /*
1806:          nc - number of components per grid point
1807:          col - number of colors needed in one direction for single component problem

1809:   */
1810:   DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1811:   col  = 2*s + 1;
1812:   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\
1813:                  by 2*stencil_width + 1\n");
1814:   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\
1815:                  by 2*stencil_width + 1\n");
1816:   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\
1817:                  by 2*stencil_width + 1\n");

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

1823:   PetscMalloc1(col*col*col*nc,&cols);
1824:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

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

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

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

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