Actual source code: fdda.c

petsc-3.5.4 2015-05-23
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,&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,&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,&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,&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,&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,&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,&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,&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:   MatView(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:         MPI_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 = dd->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:   /* determine the matrix preallocation information */
823:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
824:   for (i=xs; i<xs+nx; i++) {

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

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

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

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

854:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1017: /* ---------------------------------------------------------------------------------*/

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

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

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

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

1046:   PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1047:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

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

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

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

1131: /* ---------------------------------------------------------------------------------*/

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

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

1152:   /*
1153:          nc - number of components per grid point
1154:          col - number of colors needed in one direction for single component problem

1156:   */
1157:   DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1158:   col  = 2*s + 1;

1160:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1161:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

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

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

1195:   MatSeqAIJSetPreallocation(J,0,cols);
1196:   MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1197:   PetscFree2(cols,ocols);

1199:   DMGetLocalToGlobalMapping(da,&ltog);
1200:   MatSetLocalToGlobalMapping(J,ltog,ltog);

1202:   /*
1203:     For each node in the grid: we get the neighbors in the local (on processor ordering
1204:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1205:     PETSc ordering.
1206:   */
1207:   if (!da->prealloc_only) {
1208:     PetscMalloc1(col*nc*nc,&cols);
1209:     PetscCalloc1(col*nc*nc,&values);

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

1293: /* ---------------------------------------------------------------------------------*/

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

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

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

1316:   DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1317:   DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

1319:   MatSetBlockSize(J,nc);
1320:   MatSeqAIJSetPreallocation(J,col*nc,0);
1321:   MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);

1323:   DMGetLocalToGlobalMapping(da,&ltog);
1324:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

1385:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

1413:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

1478:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

1513:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

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

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

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

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

1606:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

1635:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

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

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

1700:   /* create the matrix */
1701:   PetscMalloc1(col*col*col,&cols);

1703:   DMGetLocalToGlobalMapping(da,&ltog);

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

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

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

1739:   MatSetLocalToGlobalMapping(J,ltog,ltog);

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

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

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

1784: /* ---------------------------------------------------------------------------------*/

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

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

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

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

1821:   PetscMalloc1(col*col*col*nc,&cols);
1822:   DMGetLocalToGlobalMapping(da,&ltog);

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


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

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

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

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

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

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