Actual source code: fdda.c

      1: #define PETSCDM_DLL
      2:
     3:  #include src/dm/da/daimpl.h
     4:  #include petscmat.h

      7: EXTERN PetscErrorCode DAGetColoring1d_MPIAIJ(DA,ISColoringType,ISColoring *);
      8: EXTERN PetscErrorCode DAGetColoring2d_MPIAIJ(DA,ISColoringType,ISColoring *);
      9: EXTERN PetscErrorCode DAGetColoring2d_5pt_MPIAIJ(DA,ISColoringType,ISColoring *);
     10: EXTERN PetscErrorCode DAGetColoring3d_MPIAIJ(DA,ISColoringType,ISColoring *);

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

     20: static PetscErrorCode DASetBlockFills_Private(PetscInt *dfill,PetscInt w,PetscInt **rfill)
     21: {
     23:   PetscInt       i,j,nz,*fill;

     26:   /* count number nonzeros */
     27:   nz = 0;
     28:   for (i=0; i<w; i++) {
     29:     for (j=0; j<w; j++) {
     30:       if (dfill[w*i+j]) nz++;
     31:     }
     32:   }
     33:   PetscMalloc((nz + w + 1)*sizeof(PetscInt),&fill);
     34:   /* construct modified CSR storage of nonzero structure */
     35:   nz = w + 1;
     36:   for (i=0; i<w; i++) {
     37:     fill[i] = nz;
     38:     for (j=0; j<w; j++) {
     39:       if (dfill[w*i+j]) {
     40:         fill[nz] = j;
     41:         nz++;
     42:       }
     43:     }
     44:   }
     45:   fill[w] = nz;
     46:
     47:   *rfill = fill;
     48:   return(0);
     49: }

     53: /*@C
     54:     DASetBlockFills - Sets the fill pattern in each block for a multi-component problem
     55:     of the matrix returned by DAGetMatrix().

     57:     Collective on DA

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

     65:     Level: developer

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

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

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

     82:    Contributed by Glenn Hammond

     84: .seealso DAGetMatrix(), DASetGetMatrix()

     86: @*/
     87: PetscErrorCode PETSCDM_DLLEXPORT DASetBlockFills(DA da,PetscInt *dfill,PetscInt *ofill)
     88: {

     92:   if (dfill) {
     93:     DASetBlockFills_Private(dfill,da->w,&da->dfill);
     94:   }
     95:   DASetBlockFills_Private(ofill,da->w,&da->ofill);
     96:   return(0);
     97: }

    102: /*@C
    103:     DAGetColoring - Gets the coloring required for computing the Jacobian via
    104:     finite differences on a function defined using a stencil on the DA.

    106:     Collective on DA

    108:     Input Parameter:
    109: +   da - the distributed array
    110: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GHOSTED

    112:     Output Parameters:
    113: .   coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)

    115:     Level: advanced

    117:     Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used
    118:    for efficient (parallel or thread based) triangular solves etc is NOT yet
    119:    available.

    122: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring

    124: @*/
    125: PetscErrorCode PETSCDM_DLLEXPORT DAGetColoring(DA da,ISColoringType ctype,ISColoring *coloring)
    126: {
    128:   PetscInt       dim;

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

    148:   /*
    149:          nc - number of components per grid point
    150:          col - number of colors needed in one direction for single component problem
    151:
    152:   */
    153:   DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);

    155:   /*
    156:      We do not provide a getcoloring function in the DA operations because
    157:    the basic DA does not know about matrices. We think of DA as being more
    158:    more low-level then matrices.
    159:   */
    160:   if (dim == 1) {
    161:     DAGetColoring1d_MPIAIJ(da,ctype,coloring);
    162:   } else if (dim == 2) {
    163:      DAGetColoring2d_MPIAIJ(da,ctype,coloring);
    164:   } else if (dim == 3) {
    165:      DAGetColoring3d_MPIAIJ(da,ctype,coloring);
    166:   } else {
    167:       SETERRQ1(PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
    168:   }
    169:   return(0);
    170: }

    172: /* ---------------------------------------------------------------------------------*/

    176: PetscErrorCode DAGetColoring2d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
    177: {
    178:   PetscErrorCode         ierr;
    179:   PetscInt               xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
    180:   PetscMPIInt            size;
    181:   MPI_Comm               comm;
    182:   DAPeriodicType         wrap;
    183:   DAStencilType          st;
    184:   ISColoringValue        *colors;

    187:   /*
    188:          nc - number of components per grid point
    189:          col - number of colors needed in one direction for single component problem
    190:
    191:   */
    192:   DAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&wrap,&st);
    193:   col    = 2*s + 1;
    194:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    195:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    196:   PetscObjectGetComm((PetscObject)da,&comm);
    197:   MPI_Comm_size(comm,&size);

    199:   /* special case as taught to us by Paul Hovland */
    200:   if (st == DA_STENCIL_STAR && s == 1) {
    201:     DAGetColoring2d_5pt_MPIAIJ(da,ctype,coloring);
    202:   } else {

    204:     if (DAXPeriodic(wrap) && (m % col)){
    205:       SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
    206:                  by 2*stencil_width + 1\n");
    207:     }
    208:     if (DAYPeriodic(wrap) && (n % col)){
    209:       SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
    210:                  by 2*stencil_width + 1\n");
    211:     }
    212:     if (ctype == IS_COLORING_LOCAL) {
    213:       if (!da->localcoloring) {
    214:         PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
    215:         ii = 0;
    216:         for (j=ys; j<ys+ny; j++) {
    217:           for (i=xs; i<xs+nx; i++) {
    218:             for (k=0; k<nc; k++) {
    219:               colors[ii++] = k + nc*((i % col) + col*(j % col));
    220:             }
    221:           }
    222:         }
    223:         ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
    224:       }
    225:       *coloring = da->localcoloring;
    226:     } else if (ctype == IS_COLORING_GHOSTED) {
    227:       if (!da->ghostedcoloring) {
    228:         PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
    229:         ii = 0;
    230:         for (j=gys; j<gys+gny; j++) {
    231:           for (i=gxs; i<gxs+gnx; i++) {
    232:             for (k=0; k<nc; k++) {
    233:               /* the complicated stuff is to handle periodic boundaries */
    234:               colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
    235:             }
    236:           }
    237:         }
    238:         ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
    239:         ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
    240:       }
    241:       *coloring = da->ghostedcoloring;
    242:     } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
    243:   }
    244:   ISColoringReference(*coloring);
    245:   return(0);
    246: }

    248: /* ---------------------------------------------------------------------------------*/

    252: PetscErrorCode DAGetColoring3d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
    253: {
    254:   PetscErrorCode  ierr;
    255:   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;
    256:   PetscMPIInt     size;
    257:   MPI_Comm        comm;
    258:   DAPeriodicType  wrap;
    259:   DAStencilType   st;
    260:   ISColoringValue *colors;

    263:   /*
    264:          nc - number of components per grid point
    265:          col - number of colors needed in one direction for single component problem
    266:
    267:   */
    268:   DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&wrap,&st);
    269:   col    = 2*s + 1;
    270:   if (DAXPeriodic(wrap) && (m % col)){
    271:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
    272:                  by 2*stencil_width + 1\n");
    273:   }
    274:   if (DAYPeriodic(wrap) && (n % col)){
    275:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
    276:                  by 2*stencil_width + 1\n");
    277:   }
    278:   if (DAZPeriodic(wrap) && (p % col)){
    279:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
    280:                  by 2*stencil_width + 1\n");
    281:   }

    283:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
    284:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
    285:   PetscObjectGetComm((PetscObject)da,&comm);
    286:   MPI_Comm_size(comm,&size);

    288:   /* create the coloring */
    289:   if (ctype == IS_COLORING_LOCAL) {
    290:     if (!da->localcoloring) {
    291:       PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);
    292:       ii = 0;
    293:       for (k=zs; k<zs+nz; k++) {
    294:         for (j=ys; j<ys+ny; j++) {
    295:           for (i=xs; i<xs+nx; i++) {
    296:             for (l=0; l<nc; l++) {
    297:               colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
    298:             }
    299:           }
    300:         }
    301:       }
    302:       ISColoringCreate(comm,nc*nx*ny*nz,colors,&da->localcoloring);
    303:     }
    304:     *coloring = da->localcoloring;
    305:   } else if (ctype == IS_COLORING_GHOSTED) {
    306:     if (!da->ghostedcoloring) {
    307:       PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);
    308:       ii = 0;
    309:       for (k=gzs; k<gzs+gnz; k++) {
    310:         for (j=gys; j<gys+gny; j++) {
    311:           for (i=gxs; i<gxs+gnx; i++) {
    312:             for (l=0; l<nc; l++) {
    313:               /* the complicated stuff is to handle periodic boundaries */
    314:               colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
    315:             }
    316:           }
    317:         }
    318:       }
    319:       ISColoringCreate(comm,nc*gnx*gny*gnz,colors,&da->ghostedcoloring);
    320:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
    321:     }
    322:     *coloring = da->ghostedcoloring;
    323:   } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
    324:   ISColoringReference(*coloring);
    325:   return(0);
    326: }

    328: /* ---------------------------------------------------------------------------------*/

    332: PetscErrorCode DAGetColoring1d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
    333: {
    334:   PetscErrorCode  ierr;
    335:   PetscInt        xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
    336:   PetscMPIInt     size;
    337:   MPI_Comm        comm;
    338:   DAPeriodicType  wrap;
    339:   ISColoringValue *colors;

    342:   /*
    343:          nc - number of components per grid point
    344:          col - number of colors needed in one direction for single component problem
    345:
    346:   */
    347:   DAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&wrap,0);
    348:   col    = 2*s + 1;

    350:   if (DAXPeriodic(wrap) && (m % col)) {
    351:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisible\n\
    352:                  by 2*stencil_width + 1\n");
    353:   }

    355:   DAGetCorners(da,&xs,0,0,&nx,0,0);
    356:   DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
    357:   PetscObjectGetComm((PetscObject)da,&comm);
    358:   MPI_Comm_size(comm,&size);

    360:   /* create the coloring */
    361:   if (ctype == IS_COLORING_LOCAL) {
    362:     if (!da->localcoloring) {
    363:       PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);
    364:       i1 = 0;
    365:       for (i=xs; i<xs+nx; i++) {
    366:         for (l=0; l<nc; l++) {
    367:           colors[i1++] = l + nc*(i % col);
    368:         }
    369:       }
    370:       ISColoringCreate(comm,nc*nx,colors,&da->localcoloring);
    371:     }
    372:     *coloring = da->localcoloring;
    373:   } else if (ctype == IS_COLORING_GHOSTED) {
    374:     if (!da->ghostedcoloring) {
    375:       PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);
    376:       i1 = 0;
    377:       for (i=gxs; i<gxs+gnx; i++) {
    378:         for (l=0; l<nc; l++) {
    379:           /* the complicated stuff is to handle periodic boundaries */
    380:           colors[i1++] = l + nc*(SetInRange(i,m) % col);
    381:         }
    382:       }
    383:       ISColoringCreate(comm,nc*gnx,colors,&da->ghostedcoloring);
    384:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
    385:     }
    386:     *coloring = da->ghostedcoloring;
    387:   } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
    388:   ISColoringReference(*coloring);
    389:   return(0);
    390: }

    394: PetscErrorCode DAGetColoring2d_5pt_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
    395: {
    396:   PetscErrorCode  ierr;
    397:   PetscInt        xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
    398:   MPI_Comm        comm;
    399:   DAPeriodicType  wrap;
    400:   ISColoringValue *colors;

    403:   /*
    404:          nc - number of components per grid point
    405:          col - number of colors needed in one direction for single component problem
    406:
    407:   */
    408:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,0);
    409:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    410:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    411:   PetscObjectGetComm((PetscObject)da,&comm);

    413:   if (DAXPeriodic(wrap) && (m % 5)){
    414:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
    415:                  by 5\n");
    416:   }
    417:   if (DAYPeriodic(wrap) && (n % 5)){
    418:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
    419:                  by 5\n");
    420:   }

    422:   /* create the coloring */
    423:   if (ctype == IS_COLORING_LOCAL) {
    424:     if (!da->localcoloring) {
    425:       PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
    426:       ii = 0;
    427:       for (j=ys; j<ys+ny; j++) {
    428:         for (i=xs; i<xs+nx; i++) {
    429:           for (k=0; k<nc; k++) {
    430:             colors[ii++] = k + nc*((3*j+i) % 5);
    431:           }
    432:         }
    433:       }
    434:       ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
    435:     }
    436:     *coloring = da->localcoloring;
    437:   } else if (ctype == IS_COLORING_GHOSTED) {
    438:     if (!da->ghostedcoloring) {
    439:       PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
    440:       ii = 0;
    441:       for (j=gys; j<gys+gny; j++) {
    442:         for (i=gxs; i<gxs+gnx; i++) {
    443:           for (k=0; k<nc; k++) {
    444:             colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
    445:           }
    446:         }
    447:       }
    448:       ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
    449:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
    450:     }
    451:     *coloring = da->ghostedcoloring;
    452:   } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
    453:   return(0);
    454: }

    456: /* =========================================================================== */
    457: EXTERN PetscErrorCode DAGetMatrix1d_MPIAIJ(DA,Mat);
    458: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ(DA,Mat);
    459: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA,Mat);
    460: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ(DA,Mat);
    461: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA,Mat);
    462: EXTERN PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA,Mat);
    463: EXTERN PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA,Mat);
    464: EXTERN PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA,Mat);
    465: EXTERN PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA,Mat);

    469: /*@C
    470:     DAGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for
    471:       computing the Jacobian on a function defined using the stencil set in the DA.

    473:     Collective on DA

    475:     Input Parameter:
    476: +   da - the distributed array
    477: -   mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
    478:             or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).

    480:     Output Parameters:
    481: .   J  - matrix with the correct nonzero structure
    482:         (obviously without the correct Jacobian values)

    484:     Level: advanced

    486:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
    487:        do not need to do it yourself.

    489: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills()

    491: @*/
    492: PetscErrorCode PETSCDM_DLLEXPORT DAGetMatrix(DA da, MatType mtype,Mat *J)
    493: {
    495:   PetscInt       dim,dof,nx,ny,nz,dims[3],starts[3];
    496:   Mat            A;
    497:   MPI_Comm       comm;
    498:   MatType        Atype;
    499:   void           (*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL;

    502:   /*
    503:                                   m
    504:           ------------------------------------------------------
    505:          |                                                     |
    506:          |                                                     |
    507:          |               ----------------------                |
    508:          |               |                    |                |
    509:       n  |           ny  |                    |                |
    510:          |               |                    |                |
    511:          |               .---------------------                |
    512:          |             (xs,ys)     nx                          |
    513:          |            .                                        |
    514:          |         (gxs,gys)                                   |
    515:          |                                                     |
    516:           -----------------------------------------------------
    517:   */

    519:   /*
    520:          nc - number of components per grid point
    521:          col - number of colors needed in one direction for single component problem
    522:
    523:   */
    524:   DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0);
    525:   DAGetCorners(da,0,0,0,&nx,&ny,&nz);
    526:   PetscObjectGetComm((PetscObject)da,&comm);
    527:   MatCreate(comm,&A);
    528:   MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
    529:   MatSetType(A,mtype);
    530:   MatSetFromOptions(A);
    531:   MatGetType(A,&Atype);
    532:   /*
    533:      We do not provide a getmatrix function in the DA operations because
    534:    the basic DA does not know about matrices. We think of DA as being more
    535:    more low-level than matrices. This is kind of cheating but, cause sometimes
    536:    we think of DA has higher level than matrices.

    538:      We could switch based on Atype (or mtype), but we do not since the
    539:    specialized setting routines depend only the particular preallocation
    540:    details of the matrix, not the type itself.
    541:   */
    542:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
    543:   if (!aij) {
    544:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
    545:   }
    546:   if (!aij) {
    547:     PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
    548:     if (!baij) {
    549:       PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
    550:     }
    551:     if (!baij){
    552:       PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
    553:       if (!sbaij) {
    554:         PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
    555:       }
    556:       if (!sbaij) SETERRQ2(PETSC_ERR_SUP,"Not implemented for the matrix type: %s!\n" \
    557:                            "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
    558:     }
    559:   }
    560:   if (aij) {
    561:     if (dim == 1) {
    562:       DAGetMatrix1d_MPIAIJ(da,A);
    563:     } else if (dim == 2) {
    564:       if (da->ofill) {
    565:         DAGetMatrix2d_MPIAIJ_Fill(da,A);
    566:       } else {
    567:         DAGetMatrix2d_MPIAIJ(da,A);
    568:       }
    569:     } else if (dim == 3) {
    570:       if (da->ofill) {
    571:         DAGetMatrix3d_MPIAIJ_Fill(da,A);
    572:       } else {
    573:         DAGetMatrix3d_MPIAIJ(da,A);
    574:       }
    575:     }
    576:   } else if (baij) {
    577:     if (dim == 2) {
    578:       DAGetMatrix2d_MPIBAIJ(da,A);
    579:     } else if (dim == 3) {
    580:       DAGetMatrix3d_MPIBAIJ(da,A);
    581:     } else {
    582:       SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
    583:                                "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
    584:     }
    585:   } else if (sbaij) {
    586:     if (dim == 2) {
    587:       DAGetMatrix2d_MPISBAIJ(da,A);
    588:     } else if (dim == 3) {
    589:       DAGetMatrix3d_MPISBAIJ(da,A);
    590:     } else {
    591:       SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
    592:                                "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
    593:     }
    594:   }
    595:   DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
    596:   MatSetStencil(A,dim,dims,starts,dof);
    597:   *J = A;
    598:   return(0);
    599: }

    601: /* ---------------------------------------------------------------------------------*/
    604: PetscErrorCode DAGetMatrix2d_MPIAIJ(DA da,Mat J)
    605: {
    606:   PetscErrorCode         ierr;
    607:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
    608:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
    609:   MPI_Comm               comm;
    610:   PetscScalar            *values;
    611:   DAPeriodicType         wrap;
    612:   ISLocalToGlobalMapping ltog;
    613:   DAStencilType          st;

    616:   /*
    617:          nc - number of components per grid point
    618:          col - number of colors needed in one direction for single component problem
    619:
    620:   */
    621:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
    622:   col = 2*s + 1;
    623:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    624:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    625:   PetscObjectGetComm((PetscObject)da,&comm);

    627:   PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
    628:   PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
    629:   PetscMalloc(nc*sizeof(PetscInt),&rows);
    630:   PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
    631:   DAGetISLocalToGlobalMapping(da,&ltog);
    632:
    633:   /* determine the matrix preallocation information */
    634:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
    635:   for (i=xs; i<xs+nx; i++) {

    637:     pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    638:     pend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));

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

    643:       lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    644:       lend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));

    646:       cnt  = 0;
    647:       for (k=0; k<nc; k++) {
    648:         for (l=lstart; l<lend+1; l++) {
    649:           for (p=pstart; p<pend+1; p++) {
    650:             if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
    651:               cols[cnt++]  = k + nc*(slot + gnx*l + p);
    652:             }
    653:           }
    654:         }
    655:         rows[k] = k + nc*(slot);
    656:       }
    657:       MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
    658:     }
    659:   }
    660:   MatSeqAIJSetPreallocation(J,0,dnz);
    661:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
    662:   MatSetBlockSize(J,nc);
    663:   MatPreallocateFinalize(dnz,onz);

    665:   MatSetLocalToGlobalMapping(J,ltog);

    667:   /*
    668:     For each node in the grid: we get the neighbors in the local (on processor ordering
    669:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    670:     PETSc ordering.
    671:   */
    672:   for (i=xs; i<xs+nx; i++) {
    673:
    674:     pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    675:     pend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    676:
    677:     for (j=ys; j<ys+ny; j++) {
    678:       slot = i - gxs + gnx*(j - gys);
    679:
    680:       lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    681:       lend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));

    683:       cnt  = 0;
    684:       for (k=0; k<nc; k++) {
    685:         for (l=lstart; l<lend+1; l++) {
    686:           for (p=pstart; p<pend+1; p++) {
    687:             if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
    688:               cols[cnt++]  = k + nc*(slot + gnx*l + p);
    689:             }
    690:           }
    691:         }
    692:         rows[k]      = k + nc*(slot);
    693:       }
    694:       MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
    695:     }
    696:   }
    697:   PetscFree(values);
    698:   PetscFree(rows);
    699:   PetscFree(cols);
    700:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    701:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    702:   return(0);
    703: }

    707: PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA da,Mat J)
    708: {
    709:   PetscErrorCode         ierr;
    710:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    711:   PetscInt               m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
    712:   PetscInt               lstart,lend,pstart,pend,*dnz,*onz;
    713:   PetscInt               ifill_col,*ofill = da->ofill, *dfill = da->dfill;
    714:   MPI_Comm               comm;
    715:   PetscScalar            *values;
    716:   DAPeriodicType         wrap;
    717:   ISLocalToGlobalMapping ltog;
    718:   DAStencilType          st;

    721:   /*
    722:          nc - number of components per grid point
    723:          col - number of colors needed in one direction for single component problem
    724:
    725:   */
    726:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
    727:   col = 2*s + 1;
    728:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    729:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    730:   PetscObjectGetComm((PetscObject)da,&comm);

    732:   PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
    733:   PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
    734:   PetscMalloc(nc*sizeof(PetscInt),&rows);
    735:   PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
    736:   DAGetISLocalToGlobalMapping(da,&ltog);
    737:
    738:   /* determine the matrix preallocation information */
    739:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
    740:   for (i=xs; i<xs+nx; i++) {

    742:     pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    743:     pend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));

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

    748:       lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    749:       lend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));

    751:       for (k=0; k<nc; k++) {
    752:         cnt  = 0;
    753:         for (l=lstart; l<lend+1; l++) {
    754:           for (p=pstart; p<pend+1; p++) {
    755:             if (l || p) {
    756:               if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
    757:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
    758:                   cols[cnt++]  = ofill[ifill_col] + nc*(slot + gnx*l + p);
    759:               }
    760:             } else {
    761:               if (dfill) {
    762:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
    763:                   cols[cnt++]  = dfill[ifill_col] + nc*(slot + gnx*l + p);
    764:               } else {
    765:                 for (ifill_col=0; ifill_col<nc; ifill_col++)
    766:                   cols[cnt++]  = ifill_col + nc*(slot + gnx*l + p);
    767:               }
    768:             }
    769:           }
    770:         }
    771:         rows[0] = k + nc*(slot);
    772:         MatPreallocateSetLocal(ltog,1,rows,cnt,cols,dnz,onz);
    773:       }
    774:     }
    775:   }
    776:   MatSeqAIJSetPreallocation(J,0,dnz);
    777:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
    778:   MatPreallocateFinalize(dnz,onz);
    779:   MatSetLocalToGlobalMapping(J,ltog);

    781:   /*
    782:     For each node in the grid: we get the neighbors in the local (on processor ordering
    783:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    784:     PETSc ordering.
    785:   */
    786:   for (i=xs; i<xs+nx; i++) {
    787:
    788:     pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    789:     pend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    790:
    791:     for (j=ys; j<ys+ny; j++) {
    792:       slot = i - gxs + gnx*(j - gys);
    793:
    794:       lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    795:       lend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));

    797:       for (k=0; k<nc; k++) {
    798:         cnt  = 0;
    799:         for (l=lstart; l<lend+1; l++) {
    800:           for (p=pstart; p<pend+1; p++) {
    801:             if (l || p) {
    802:               if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star */
    803:                 for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
    804:                   cols[cnt++]  = ofill[ifill_col] + nc*(slot + gnx*l + p);
    805:               }
    806:             } else {
    807:               if (dfill) {
    808:                 for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
    809:                   cols[cnt++]  = dfill[ifill_col] + nc*(slot + gnx*l + p);
    810:               } else {
    811:                 for (ifill_col=0; ifill_col<nc; ifill_col++)
    812:                   cols[cnt++]  = ifill_col + nc*(slot + gnx*l + p);
    813:               }
    814:             }
    815:           }
    816:         }
    817:         rows[0]      = k + nc*(slot);
    818:         MatSetValuesLocal(J,1,rows,cnt,cols,values,INSERT_VALUES);
    819:       }
    820:     }
    821:   }
    822:   PetscFree(values);
    823:   PetscFree(rows);
    824:   PetscFree(cols);
    825:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    826:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    827:   return(0);
    828: }

    830: /* ---------------------------------------------------------------------------------*/

    834: PetscErrorCode DAGetMatrix3d_MPIAIJ(DA da,Mat J)
    835: {
    836:   PetscErrorCode         ierr;
    837:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    838:   PetscInt               m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
    839:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
    840:   MPI_Comm               comm;
    841:   PetscScalar            *values;
    842:   DAPeriodicType         wrap;
    843:   ISLocalToGlobalMapping ltog;
    844:   DAStencilType          st;

    847:   /*
    848:          nc - number of components per grid point
    849:          col - number of colors needed in one direction for single component problem
    850:
    851:   */
    852:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
    853:   col    = 2*s + 1;

    855:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
    856:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
    857:   PetscObjectGetComm((PetscObject)da,&comm);

    859:   PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
    860:   PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
    861:   PetscMalloc(nc*sizeof(PetscInt),&rows);
    862:   PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
    863:   DAGetISLocalToGlobalMapping(da,&ltog);

    865:   /* determine the matrix preallocation information */
    866:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
    867:   for (i=xs; i<xs+nx; i++) {
    868:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    869:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    870:     for (j=ys; j<ys+ny; j++) {
    871:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    872:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    873:       for (k=zs; k<zs+nz; k++) {
    874:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    875:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    876:
    877:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    878:
    879:         cnt  = 0;
    880:         for (l=0; l<nc; l++) {
    881:           for (ii=istart; ii<iend+1; ii++) {
    882:             for (jj=jstart; jj<jend+1; jj++) {
    883:               for (kk=kstart; kk<kend+1; kk++) {
    884:                 if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
    885:                   cols[cnt++]  = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    886:                 }
    887:               }
    888:             }
    889:           }
    890:           rows[l] = l + nc*(slot);
    891:         }
    892:         MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
    893:       }
    894:     }
    895:   }
    896:   MatSeqAIJSetPreallocation(J,0,dnz);
    897:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
    898:   MatPreallocateFinalize(dnz,onz);
    899:   MatSetBlockSize(J,nc);CHKERRQ(ierr)
    900:   MatSetLocalToGlobalMapping(J,ltog);

    902:   /*
    903:     For each node in the grid: we get the neighbors in the local (on processor ordering
    904:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    905:     PETSc ordering.
    906:   */
    907:   for (i=xs; i<xs+nx; i++) {
    908:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    909:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    910:     for (j=ys; j<ys+ny; j++) {
    911:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    912:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    913:       for (k=zs; k<zs+nz; k++) {
    914:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    915:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    916:
    917:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    918:
    919:         cnt  = 0;
    920:         for (l=0; l<nc; l++) {
    921:           for (ii=istart; ii<iend+1; ii++) {
    922:             for (jj=jstart; jj<jend+1; jj++) {
    923:               for (kk=kstart; kk<kend+1; kk++) {
    924:                 if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
    925:                   cols[cnt++]  = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    926:                 }
    927:               }
    928:             }
    929:           }
    930:           rows[l]      = l + nc*(slot);
    931:         }
    932:         MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
    933:       }
    934:     }
    935:   }
    936:   PetscFree(values);
    937:   PetscFree(rows);
    938:   PetscFree(cols);
    939:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    940:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    941:   return(0);
    942: }

    944: /* ---------------------------------------------------------------------------------*/

    948: PetscErrorCode DAGetMatrix1d_MPIAIJ(DA da,Mat J)
    949: {
    950:   PetscErrorCode         ierr;
    951:   PetscInt               xs,nx,i,i1,slot,gxs,gnx;
    952:   PetscInt               m,dim,s,*cols,nc,*rows,col,cnt,l;
    953:   PetscInt               istart,iend;
    954:   PetscScalar            *values;
    955:   DAPeriodicType         wrap;
    956:   ISLocalToGlobalMapping ltog;

    959:   /*
    960:          nc - number of components per grid point
    961:          col - number of colors needed in one direction for single component problem
    962:
    963:   */
    964:   DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
    965:   col    = 2*s + 1;

    967:   DAGetCorners(da,&xs,0,0,&nx,0,0);
    968:   DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);

    970:   MatSeqAIJSetPreallocation(J,col*nc,0);
    971:   MatMPIAIJSetPreallocation(J,col*nc,0,0,0);
    972:   MatSetBlockSize(J,nc);
    973:   PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
    974:   PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
    975:   PetscMalloc(nc*sizeof(PetscInt),&rows);
    976:   PetscMalloc(col*nc*sizeof(PetscInt),&cols);
    977:
    978:   DAGetISLocalToGlobalMapping(da,&ltog);
    979:   MatSetLocalToGlobalMapping(J,ltog);
    980:
    981:   /*
    982:     For each node in the grid: we get the neighbors in the local (on processor ordering
    983:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    984:     PETSc ordering.
    985:   */
    986:   for (i=xs; i<xs+nx; i++) {
    987:     istart = PetscMax(-s,gxs - i);
    988:     iend   = PetscMin(s,gxs + gnx - i - 1);
    989:     slot   = i - gxs;
    990:
    991:     cnt  = 0;
    992:     for (l=0; l<nc; l++) {
    993:       for (i1=istart; i1<iend+1; i1++) {
    994:         cols[cnt++] = l + nc*(slot + i1);
    995:       }
    996:       rows[l]      = l + nc*(slot);
    997:     }
    998:     MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
    999:   }
    1000:   PetscFree(values);
    1001:   PetscFree(rows);
    1002:   PetscFree(cols);
    1003:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1004:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1005:   return(0);
    1006: }

    1010: PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA da,Mat J)
    1011: {
    1012:   PetscErrorCode         ierr;
    1013:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    1014:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
    1015:   PetscInt               istart,iend,jstart,jend,ii,jj;
    1016:   MPI_Comm               comm;
    1017:   PetscScalar            *values;
    1018:   DAPeriodicType         wrap;
    1019:   DAStencilType          st;
    1020:   ISLocalToGlobalMapping ltog,ltogb;

    1023:   /*
    1024:      nc - number of components per grid point
    1025:      col - number of colors needed in one direction for single component problem
    1026:   */
    1027:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
    1028:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
    1029:   col = 2*s + 1;

    1031:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    1032:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    1033:   PetscObjectGetComm((PetscObject)da,&comm);

    1035:   PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
    1036:   PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
    1037:   PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);

    1039:   DAGetISLocalToGlobalMapping(da,&ltog);
    1040:   DAGetISLocalToGlobalMappingBlck(da,&ltogb);

    1042:   /* determine the matrix preallocation information */
    1043:   MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
    1044:   for (i=xs; i<xs+nx; i++) {
    1045:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1046:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1047:     for (j=ys; j<ys+ny; j++) {
    1048:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1049:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1050:       slot = i - gxs + gnx*(j - gys);

    1052:       /* Find block columns in block row */
    1053:       cnt  = 0;
    1054:       if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1055:         for (ii=istart; ii<iend+1; ii++) {
    1056:           for (jj=jstart; jj<jend+1; jj++) {
    1057:             cols[cnt++]  = slot + ii + gnx*jj;
    1058:           }
    1059:         }
    1060:       } else {  /* Star stencil */
    1061:         cnt  = 0;
    1062:         for (ii=istart; ii<iend+1; ii++) {
    1063:           if (ii) { /* jj must be zero */
    1064:             cols[cnt++]  = slot + ii;
    1065:           } else {
    1066:             for (jj=jstart; jj<jend+1; jj++) { /* ii must be zero */
    1067:               cols[cnt++] = slot + gnx*jj;
    1068:             }
    1069:           }
    1070:         }
    1071:       }
    1072:       MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
    1073:     }
    1074:   }
    1075:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
    1076:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
    1077:   MatPreallocateFinalize(dnz,onz);

    1079:   MatSetLocalToGlobalMapping(J,ltog);
    1080:   MatSetLocalToGlobalMappingBlock(J,ltogb);

    1082:   /*
    1083:     For each node in the grid: we get the neighbors in the local (on processor ordering
    1084:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    1085:     PETSc ordering.
    1086:   */
    1087:   for (i=xs; i<xs+nx; i++) {
    1088:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1089:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1090:     for (j=ys; j<ys+ny; j++) {
    1091:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1092:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1093:       slot = i - gxs + gnx*(j - gys);
    1094:       cnt  = 0;
    1095:       if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1096:         for (ii=istart; ii<iend+1; ii++) {
    1097:           for (jj=jstart; jj<jend+1; jj++) {
    1098:             cols[cnt++]  = slot + ii + gnx*jj;
    1099:           }
    1100:         }
    1101:       } else {  /* Star stencil */
    1102:         cnt  = 0;
    1103:         for (ii=istart; ii<iend+1; ii++) {
    1104:           if (ii) { /* jj must be zero */
    1105:             cols[cnt++]  = slot + ii;
    1106:           } else {
    1107:             for (jj=jstart; jj<jend+1; jj++) {/* ii must be zero */
    1108:               cols[cnt++]  = slot + gnx*jj;
    1109:             }
    1110:           }
    1111:         }
    1112:       }
    1113:       MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
    1114:     }
    1115:   }
    1116:   PetscFree(values);
    1117:   PetscFree(cols);
    1118:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1119:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1120:   return(0);
    1121: }

    1125: PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA da,Mat J)
    1126: {
    1127:   PetscErrorCode         ierr;
    1128:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    1129:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
    1130:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
    1131:   MPI_Comm               comm;
    1132:   PetscScalar            *values;
    1133:   DAPeriodicType         wrap;
    1134:   DAStencilType          st;
    1135:   ISLocalToGlobalMapping ltog,ltogb;

    1138:   /*
    1139:          nc - number of components per grid point
    1140:          col - number of colors needed in one direction for single component problem
    1141:
    1142:   */
    1143:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
    1144:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
    1145:   col    = 2*s + 1;

    1147:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
    1148:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
    1149:   PetscObjectGetComm((PetscObject)da,&comm);

    1151:   PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
    1152:   PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
    1153:   PetscMalloc(col*col*col*sizeof(PetscInt),&cols);

    1155:   DAGetISLocalToGlobalMapping(da,&ltog);
    1156:   DAGetISLocalToGlobalMappingBlck(da,&ltogb);

    1158:   /* determine the matrix preallocation information */
    1159:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
    1160:   for (i=xs; i<xs+nx; i++) {
    1161:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1162:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1163:     for (j=ys; j<ys+ny; j++) {
    1164:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1165:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1166:       for (k=zs; k<zs+nz; k++) {
    1167:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    1168:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));

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

    1172:         /* Find block columns in block row */
    1173:         cnt  = 0;
    1174:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1175:           for (ii=istart; ii<iend+1; ii++) {
    1176:             for (jj=jstart; jj<jend+1; jj++) {
    1177:               for (kk=kstart; kk<kend+1; kk++) {
    1178:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
    1179:               }
    1180:             }
    1181:           }
    1182:         } else {  /* Star stencil */
    1183:           cnt  = 0;
    1184:           for (ii=istart; ii<iend+1; ii++) {
    1185:             if (ii) {
    1186:               /* jj and kk must be zero */
    1187:               /* cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk; */
    1188:               cols[cnt++]  = slot + ii;
    1189:             } else {
    1190:               for (jj=jstart; jj<jend+1; jj++) {
    1191:                 if (jj) {
    1192:                   /* ii and kk must be zero */
    1193:                   cols[cnt++]  = slot + gnx*jj;
    1194:                 } else {
    1195:                   /* ii and jj must be zero */
    1196:                   for (kk=kstart; kk<kend+1; kk++) {
    1197:                     cols[cnt++]  = slot + gnx*gny*kk;
    1198:                   }
    1199:                 }
    1200:               }
    1201:             }
    1202:           }
    1203:         }
    1204:         MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
    1205:       }
    1206:     }
    1207:   }
    1208:   MatSeqBAIJSetPreallocation(J,nc,0,dnz);
    1209:   MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
    1210:   MatPreallocateFinalize(dnz,onz);

    1212:   MatSetLocalToGlobalMapping(J,ltog);
    1213:   MatSetLocalToGlobalMappingBlock(J,ltogb);

    1215:   /*
    1216:     For each node in the grid: we get the neighbors in the local (on processor ordering
    1217:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    1218:     PETSc ordering.
    1219:   */

    1221:   for (i=xs; i<xs+nx; i++) {
    1222:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1223:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1224:     for (j=ys; j<ys+ny; j++) {
    1225:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1226:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1227:       for (k=zs; k<zs+nz; k++) {
    1228:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    1229:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    1230:
    1231:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    1232:
    1233:         cnt  = 0;
    1234:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1235:           for (ii=istart; ii<iend+1; ii++) {
    1236:             for (jj=jstart; jj<jend+1; jj++) {
    1237:               for (kk=kstart; kk<kend+1; kk++) {
    1238:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
    1239:               }
    1240:             }
    1241:           }
    1242:         } else {  /* Star stencil */
    1243:           cnt  = 0;
    1244:           for (ii=istart; ii<iend+1; ii++) {
    1245:             if (ii) {
    1246:               /* jj and kk must be zero */
    1247:               cols[cnt++]  = slot + ii;
    1248:             } else {
    1249:               for (jj=jstart; jj<jend+1; jj++) {
    1250:                 if (jj) {
    1251:                   /* ii and kk must be zero */
    1252:                   cols[cnt++]  = slot + gnx*jj;
    1253:                 } else {
    1254:                   /* ii and jj must be zero */
    1255:                   for (kk=kstart; kk<kend+1; kk++) {
    1256:                     cols[cnt++]  = slot + gnx*gny*kk;
    1257:                   }
    1258:                 }
    1259:               }
    1260:             }
    1261:           }
    1262:         }
    1263:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
    1264:       }
    1265:     }
    1266:   }
    1267:   PetscFree(values);
    1268:   PetscFree(cols);
    1269:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1270:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1271:   return(0);
    1272: }
    1275: PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA da,Mat J)
    1276: {
    1277:   PetscErrorCode         ierr;
    1278:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    1279:   PetscInt               m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
    1280:   PetscInt               iend,jend,ii,jj;
    1281:   MPI_Comm               comm;
    1282:   PetscScalar            *values;
    1283:   DAPeriodicType         wrap;
    1284:   DAStencilType          st;
    1285:   ISLocalToGlobalMapping ltog,ltogb;

    1288:   /*
    1289:      nc - number of components per grid point
    1290:      col - number of colors needed in one direction for single component problem
    1291:   */
    1292:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
    1293:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
    1294:   col = 2*s + 1;

    1296:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
    1297:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
    1298:   PetscObjectGetComm((PetscObject)da,&comm);

    1300:   PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
    1301:   PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
    1302:   PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);

    1304:   DAGetISLocalToGlobalMapping(da,&ltog);
    1305:   DAGetISLocalToGlobalMappingBlck(da,&ltogb);

    1307:   /* determine the matrix preallocation information */
    1308:   MatPreallocateSymmetricInitialize(comm,nx*ny,nx*ny,dnz,onz);
    1309:   for (i=xs; i<xs+nx; i++) {
    1310:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1311:     for (j=ys; j<ys+ny; j++) {
    1312:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1313:       slot = i - gxs + gnx*(j - gys);

    1315:       /* Find block columns in block row */
    1316:       cnt  = 0;
    1317:       if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1318:         for (ii=0; ii<iend+1; ii++) {
    1319:           for (jj=0; jj<jend+1; jj++) {
    1320:             cols[cnt++]  = slot + ii + gnx*jj;
    1321:           }
    1322:         }
    1323:       } else {  /* Star stencil */
    1324:         cnt  = 0;
    1325:         for (ii=0; ii<iend+1; ii++) {
    1326:           if (ii) { /* jj must be zero */
    1327:             cols[cnt++]  = slot + ii;
    1328:           } else {
    1329:             for (jj=0; jj<jend+1; jj++) { /* ii must be zero */
    1330:               cols[cnt++] = slot + gnx*jj;
    1331:             }
    1332:           }
    1333:         }
    1334:       }
    1335:       MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
    1336:     }
    1337:   }
    1338:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
    1339:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
    1340:   MatPreallocateFinalize(dnz,onz);

    1342:   MatSetLocalToGlobalMapping(J,ltog);
    1343:   MatSetLocalToGlobalMappingBlock(J,ltogb);

    1345:   /*
    1346:     For each node in the grid: we get the neighbors in the local (on processor ordering
    1347:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    1348:     PETSc ordering.
    1349:   */
    1350:   for (i=xs; i<xs+nx; i++) {
    1351:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1352:     for (j=ys; j<ys+ny; j++) {
    1353:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1354:       slot = i - gxs + gnx*(j - gys);
    1355:       cnt  = 0;
    1356:       if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1357:         for (ii=0; ii<iend+1; ii++) {
    1358:           for (jj=0; jj<jend+1; jj++) {
    1359:             cols[cnt++]  = slot + ii + gnx*jj;
    1360:           }
    1361:         }
    1362:       } else {  /* Star stencil */
    1363:         cnt  = 0;
    1364:         for (ii=0; ii<iend+1; ii++) {
    1365:           if (ii) { /* jj must be zero */
    1366:             cols[cnt++]  = slot + ii;
    1367:           } else {
    1368:             for (jj=0; jj<jend+1; jj++) {/* ii must be zero */
    1369:               cols[cnt++]  = slot + gnx*jj;
    1370:             }
    1371:           }
    1372:         }
    1373:       }
    1374:       MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
    1375:     }
    1376:   }
    1377:   PetscFree(values);
    1378:   PetscFree(cols);
    1379:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1380:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1381:   return(0);
    1382: }

    1386: PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA da,Mat J)
    1387: {
    1388:   PetscErrorCode         ierr;
    1389:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    1390:   PetscInt               m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
    1391:   PetscInt               iend,jend,kend,zs,nz,gzs,gnz,ii,jj,kk;
    1392:   MPI_Comm               comm;
    1393:   PetscScalar            *values;
    1394:   DAPeriodicType         wrap;
    1395:   DAStencilType          st;
    1396:   ISLocalToGlobalMapping ltog,ltogb;

    1399:   /*
    1400:      nc - number of components per grid point
    1401:      col - number of colors needed in one direction for single component problem
    1402:   */
    1403:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
    1404:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
    1405:   col = 2*s + 1;

    1407:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
    1408:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
    1409:   PetscObjectGetComm((PetscObject)da,&comm);

    1411:   /* create the matrix */
    1412:   PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
    1413:   PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
    1414:   PetscMalloc(col*col*col*sizeof(PetscInt),&cols);

    1416:   DAGetISLocalToGlobalMapping(da,&ltog);
    1417:   DAGetISLocalToGlobalMappingBlck(da,&ltogb);

    1419:   /* determine the matrix preallocation information */
    1420:   MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
    1421:   for (i=xs; i<xs+nx; i++) {
    1422:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1423:     for (j=ys; j<ys+ny; j++) {
    1424:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1425:       for (k=zs; k<zs+nz; k++) {
    1426:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));

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

    1430:         /* Find block columns in block row */
    1431:         cnt  = 0;
    1432:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1433:           for (ii=0; ii<iend+1; ii++) {
    1434:             for (jj=0; jj<jend+1; jj++) {
    1435:               for (kk=0; kk<kend+1; kk++) {
    1436:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
    1437:               }
    1438:             }
    1439:           }
    1440:         } else {  /* Star stencil */
    1441:           cnt  = 0;
    1442:           for (ii=0; ii<iend+1; ii++) {
    1443:             if (ii) {
    1444:               /* jj and kk must be zero */
    1445:               cols[cnt++]  = slot + ii;
    1446:             } else {
    1447:               for (jj=0; jj<jend+1; jj++) {
    1448:                 if (jj) {
    1449:                   /* ii and kk must be zero */
    1450:                   cols[cnt++]  = slot + gnx*jj;
    1451:                 } else {
    1452:                   /* ii and jj must be zero */
    1453:                   for (kk=0; kk<kend+1; kk++) {
    1454:                     cols[cnt++]  = slot + gnx*gny*kk;
    1455:                   }
    1456:                 }
    1457:               }
    1458:             }
    1459:           }
    1460:         }
    1461:         MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
    1462:       }
    1463:     }
    1464:   }
    1465:   MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
    1466:   MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
    1467:   MatPreallocateFinalize(dnz,onz);

    1469:   MatSetLocalToGlobalMapping(J,ltog);
    1470:   MatSetLocalToGlobalMappingBlock(J,ltogb);

    1472:   /*
    1473:     For each node in the grid: we get the neighbors in the local (on processor ordering
    1474:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    1475:     PETSc ordering.
    1476:   */

    1478:   for (i=xs; i<xs+nx; i++) {
    1479:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1480:     for (j=ys; j<ys+ny; j++) {
    1481:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1482:       for (k=zs; k<zs+nz; k++) {
    1483:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    1484:
    1485:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    1486:
    1487:         cnt  = 0;
    1488:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
    1489:           for (ii=0; ii<iend+1; ii++) {
    1490:             for (jj=0; jj<jend+1; jj++) {
    1491:               for (kk=0; kk<kend+1; kk++) {
    1492:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
    1493:               }
    1494:             }
    1495:           }
    1496:         } else {  /* Star stencil */
    1497:           cnt  = 0;
    1498:           for (ii=0; ii<iend+1; ii++) {
    1499:             if (ii) {
    1500:               /* jj and kk must be zero */
    1501:               cols[cnt++]  = slot + ii;
    1502:             } else {
    1503:               for (jj=0; jj<jend+1; jj++) {
    1504:                 if (jj) {
    1505:                   /* ii and kk must be zero */
    1506:                   cols[cnt++]  = slot + gnx*jj;
    1507:                 } else {
    1508:                   /* ii and jj must be zero */
    1509:                   for (kk=0; kk<kend+1; kk++) {
    1510:                     cols[cnt++]  = slot + gnx*gny*kk;
    1511:                   }
    1512:                 }
    1513:               }
    1514:             }
    1515:           }
    1516:         }
    1517:         MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
    1518:       }
    1519:     }
    1520:   }
    1521:   PetscFree(values);
    1522:   PetscFree(cols);
    1523:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1524:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1525:   return(0);
    1526: }

    1528: /* ---------------------------------------------------------------------------------*/

    1532: PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA da,Mat J)
    1533: {
    1534:   PetscErrorCode         ierr;
    1535:   PetscInt               xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
    1536:   PetscInt               m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
    1537:   PetscInt               istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
    1538:   PetscInt               ifill_col,*dfill = da->dfill,*ofill = da->ofill;
    1539:   MPI_Comm               comm;
    1540:   PetscScalar            *values;
    1541:   DAPeriodicType         wrap;
    1542:   ISLocalToGlobalMapping ltog;
    1543:   DAStencilType          st;

    1546:   /*
    1547:          nc - number of components per grid point
    1548:          col - number of colors needed in one direction for single component problem
    1549:
    1550:   */
    1551:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
    1552:   col    = 2*s + 1;
    1553:   if (DAXPeriodic(wrap) && (m % col)){
    1554:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
    1555:                  by 2*stencil_width + 1\n");
    1556:   }
    1557:   if (DAYPeriodic(wrap) && (n % col)){
    1558:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
    1559:                  by 2*stencil_width + 1\n");
    1560:   }
    1561:   if (DAZPeriodic(wrap) && (p % col)){
    1562:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
    1563:                  by 2*stencil_width + 1\n");
    1564:   }

    1566:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
    1567:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
    1568:   PetscObjectGetComm((PetscObject)da,&comm);

    1570:   PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
    1571:   PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
    1572:   PetscMalloc(nc*sizeof(PetscInt),&rows);
    1573:   PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
    1574:   DAGetISLocalToGlobalMapping(da,&ltog);

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

    1580:   for (i=xs; i<xs+nx; i++) {
    1581:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1582:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1583:     for (j=ys; j<ys+ny; j++) {
    1584:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1585:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1586:       for (k=zs; k<zs+nz; k++) {
    1587:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    1588:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    1589:
    1590:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    1591:
    1592:         for (l=0; l<nc; l++) {
    1593:           cnt  = 0;
    1594:           for (ii=istart; ii<iend+1; ii++) {
    1595:             for (jj=jstart; jj<jend+1; jj++) {
    1596:               for (kk=kstart; kk<kend+1; kk++) {
    1597:                 if (ii || jj || kk) {
    1598:                   if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
    1599:                     for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
    1600:                       cols[cnt++]  = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1601:                   }
    1602:                 } else {
    1603:                   if (dfill) {
    1604:                     for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
    1605:                       cols[cnt++]  = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1606:                   } else {
    1607:                     for (ifill_col=0; ifill_col<nc; ifill_col++)
    1608:                       cols[cnt++]  = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1609:                   }
    1610:                 }
    1611:               }
    1612:             }
    1613:           }
    1614:           rows[0] = l + nc*(slot);
    1615:           MatPreallocateSetLocal(ltog,1,rows,cnt,cols,dnz,onz);
    1616:         }
    1617:       }
    1618:     }
    1619:   }
    1620:   MatSeqAIJSetPreallocation(J,0,dnz);
    1621:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
    1622:   MatPreallocateFinalize(dnz,onz);
    1623:   MatSetLocalToGlobalMapping(J,ltog);

    1625:   /*
    1626:     For each node in the grid: we get the neighbors in the local (on processor ordering
    1627:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
    1628:     PETSc ordering.
    1629:   */
    1630:   for (i=xs; i<xs+nx; i++) {
    1631:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
    1632:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
    1633:     for (j=ys; j<ys+ny; j++) {
    1634:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
    1635:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
    1636:       for (k=zs; k<zs+nz; k++) {
    1637:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
    1638:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
    1639:
    1640:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
    1641:
    1642:         for (l=0; l<nc; l++) {
    1643:           cnt  = 0;
    1644:           for (ii=istart; ii<iend+1; ii++) {
    1645:             for (jj=jstart; jj<jend+1; jj++) {
    1646:               for (kk=kstart; kk<kend+1; kk++) {
    1647:                 if (ii || jj || kk) {
    1648:                   if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
    1649:                     for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
    1650:                       cols[cnt++]  = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1651:                   }
    1652:                 } else {
    1653:                   if (dfill) {
    1654:                     for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
    1655:                       cols[cnt++]  = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1656:                   } else {
    1657:                     for (ifill_col=0; ifill_col<nc; ifill_col++)
    1658:                       cols[cnt++]  = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
    1659:                   }
    1660:                 }
    1661:               }
    1662:             }
    1663:           }
    1664:           rows[0] = l + nc*(slot);
    1665:           MatSetValuesLocal(J,1,rows,cnt,cols,values,INSERT_VALUES);
    1666:         }
    1667:       }
    1668:     }
    1669:   }
    1670:   PetscFree(values);
    1671:   PetscFree(rows);
    1672:   PetscFree(cols);
    1673:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    1674:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    1675:   return(0);
    1676: }