Actual source code: fdda.c

  1: /*$Id: fdda.c,v 1.75 2001/08/07 21:31:51 bsmith Exp $*/
  2: 
 3:  #include src/dm/da/daimpl.h
 4:  #include petscmat.h


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

 12: #undef __FUNCT__  
 14: /*@C
 15:     DAGetColoring - Gets the coloring required for computing the Jacobian via
 16:     finite differences on a function defined using a stencil on the DA.

 18:     Collective on DA

 20:     Input Parameter:
 21: +   da - the distributed array
 22: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GHOSTED

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

 27:     Level: advanced

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


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

 36: @*/
 37: int DAGetColoring(DA da,ISColoringType ctype,ISColoring *coloring)
 38: {
 39:   int        ierr,dim;

 42:   /*
 43:                                   m
 44:           ------------------------------------------------------
 45:          |                                                     |
 46:          |                                                     |
 47:          |               ----------------------                |
 48:          |               |                    |                |
 49:       n  |           yn  |                    |                |
 50:          |               |                    |                |
 51:          |               .---------------------                |
 52:          |             (xs,ys)     xn                          |
 53:          |            .                                        |
 54:          |         (gxs,gys)                                   |
 55:          |                                                     |
 56:           -----------------------------------------------------
 57:   */

 59:   /*     
 60:          nc - number of components per grid point 
 61:          col - number of colors needed in one direction for single component problem
 62:   
 63:   */
 64:   DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);

 66:   /*
 67:      We do not provide a getcoloring function in the DA operations because 
 68:    the basic DA does not know about matrices. We think of DA as being more 
 69:    more low-level then matrices.
 70:   */
 71:   if (dim == 1) {
 72:     DAGetColoring1d_MPIAIJ(da,ctype,coloring);
 73:   } else if (dim == 2) {
 74:      DAGetColoring2d_MPIAIJ(da,ctype,coloring);
 75:   } else if (dim == 3) {
 76:      DAGetColoring3d_MPIAIJ(da,ctype,coloring);
 77:   } else {
 78:       SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
 79:   }
 80:   return(0);
 81: }

 83: /* ---------------------------------------------------------------------------------*/

 85: #undef __FUNCT__  
 87: int DAGetColoring2d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
 88: {
 89:   int                    ierr,xs,ys,nx,ny,*colors,i,j,ii,gxs,gys,gnx,gny;
 90:   int                    m,n,M,N,dim,w,s,k,nc,col,size;
 91:   MPI_Comm               comm;
 92:   DAPeriodicType         wrap;
 93:   DAStencilType          st;

 96:   /*     
 97:          nc - number of components per grid point 
 98:          col - number of colors needed in one direction for single component problem
 99:   
100:   */
101:   DAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&w,&s,&wrap,&st);
102:   nc     = w;
103:   col    = 2*s + 1;
104:   if (DAXPeriodic(wrap) && (m % col)){
105:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
106:                  by 2*stencil_width + 1n");
107:   }
108:   if (DAYPeriodic(wrap) && (n % col)){
109:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
110:                  by 2*stencil_width + 1n");
111:   }
112:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
113:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
114:   PetscObjectGetComm((PetscObject)da,&comm);
115:   MPI_Comm_size(comm,&size);

117:   /* create the coloring */
118:   if (st == DA_STENCIL_STAR && s == 1 && !DAXPeriodic(wrap) && !DAYPeriodic(wrap)) {
119:     DAGetColoring2d_5pt_MPIAIJ(da,ctype,coloring);
120:   } else if (ctype == IS_COLORING_LOCAL) {
121:     if (!da->localcoloring) {
122:     PetscMalloc(nc*nx*ny*sizeof(int),&colors);
123:     ii = 0;
124:     for (j=ys; j<ys+ny; j++) {
125:         for (i=xs; i<xs+nx; i++) {
126:           for (k=0; k<nc; k++) {
127:             colors[ii++] = k + nc*((i % col) + col*(j % col));
128:           }
129:         }
130:       }
131:       ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
132:     }
133:     *coloring = da->localcoloring;
134:   } else if (ctype == IS_COLORING_GHOSTED) {
135:     if (!da->ghostedcoloring) {
136:       PetscMalloc(nc*gnx*gny*sizeof(int),&colors);
137:       ii = 0;
138:       for (j=gys; j<gys+gny; j++) {
139:         for (i=gxs; i<gxs+gnx; i++) {
140:           for (k=0; k<nc; k++) {
141:             /* the complicated stuff is to handle periodic boundaries */
142:             colors[ii++] = k + nc*(    (((i < 0) ? m+i:((i >= m) ? i-m:i)) % col) +
143:                                      col*(((j < 0) ? n+j:((j >= n) ? j-n:j)) % col));
144:           }
145:         }
146:       }
147:       ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
148:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
149:     }
150:     *coloring = da->ghostedcoloring;
151:   } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
152:   ISColoringReference(*coloring);
153:   return(0);
154: }

156: /* ---------------------------------------------------------------------------------*/

158: #undef __FUNCT__  
160: int DAGetColoring3d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
161: {
162:   int                    ierr,xs,ys,nx,ny,*colors,i,j,gxs,gys,gnx,gny;
163:   int                    m,n,p,dim,s,k,nc,col,size,zs,gzs,ii,l,nz,gnz,M,N,P;
164:   MPI_Comm               comm;
165:   DAPeriodicType         wrap;
166:   DAStencilType          st;

169:   /*     
170:          nc - number of components per grid point 
171:          col - number of colors needed in one direction for single component problem
172:   
173:   */
174:   DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&wrap,&st);
175:   col    = 2*s + 1;
176:   if (DAXPeriodic(wrap) && (m % col)){
177:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
178:                  by 2*stencil_width + 1n");
179:   }
180:   if (DAYPeriodic(wrap) && (n % col)){
181:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
182:                  by 2*stencil_width + 1n");
183:   }
184:   if (DAZPeriodic(wrap) && (p % col)){
185:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisiblen
186:                  by 2*stencil_width + 1n");
187:   }

189:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
190:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
191:   PetscObjectGetComm((PetscObject)da,&comm);
192:   MPI_Comm_size(comm,&size);

194:   /* create the coloring */
195:   if (ctype == IS_COLORING_LOCAL) {
196:     if (!da->localcoloring) {
197:       PetscMalloc(nc*nx*ny*nz*sizeof(int),&colors);
198:       ii = 0;
199:       for (k=zs; k<zs+nz; k++) {
200:         for (j=ys; j<ys+ny; j++) {
201:           for (i=xs; i<xs+nx; i++) {
202:             for (l=0; l<nc; l++) {
203:               colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
204:             }
205:           }
206:         }
207:       }
208:       ISColoringCreate(comm,nc*nx*ny*nz,colors,&da->localcoloring);
209:     }
210:     *coloring = da->localcoloring;
211:   } else if (ctype == IS_COLORING_GHOSTED) {
212:     if (!da->ghostedcoloring) {
213:       PetscMalloc(nc*gnx*gny*gnz*sizeof(int),&colors);
214:       ii = 0;
215:       for (k=gzs; k<gzs+gnz; k++) {
216:         for (j=gys; j<gys+gny; j++) {
217:           for (i=gxs; i<gxs+gnx; i++) {
218:             for (l=0; l<nc; l++) {
219:               /* the complicated stuff is to handle periodic boundaries */
220:               colors[ii++] = l + nc*(        (((i < 0) ? m+i:((i >= m) ? i-m:i)) % col) +
221:                                            col*(((j < 0) ? n+j:((j >= n) ? j-n:j)) % col) +
222:                                      col*col*(((k < 0) ? p+k:((k >= p) ? k-p:k)) % col));
223:             }
224:           }
225:         }
226:       }
227:       ISColoringCreate(comm,nc*gnx*gny*gnz,colors,&da->ghostedcoloring);
228:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
229:     }
230:     *coloring = da->ghostedcoloring;
231:   } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
232:   ISColoringReference(*coloring);
233:   return(0);
234: }

236: /* ---------------------------------------------------------------------------------*/

238: #undef __FUNCT__  
240: int DAGetColoring1d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
241: {
242:   int                    ierr,xs,nx,*colors,i,i1,gxs,gnx,l;
243:   int                    m,M,dim,s,nc,col,size;
244:   MPI_Comm               comm;
245:   DAPeriodicType         wrap;

248:   /*     
249:          nc - number of components per grid point 
250:          col - number of colors needed in one direction for single component problem
251:   
252:   */
253:   DAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&wrap,0);
254:   col    = 2*s + 1;

256:   if (DAXPeriodic(wrap) && (m % col)) {
257:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisiblen
258:                  by 2*stencil_width + 1n");
259:   }

261:   DAGetCorners(da,&xs,0,0,&nx,0,0);
262:   DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
263:   PetscObjectGetComm((PetscObject)da,&comm);
264:   MPI_Comm_size(comm,&size);

266:   /* create the coloring */
267:   if (ctype == IS_COLORING_LOCAL) {
268:     if (!da->localcoloring) {
269:       PetscMalloc(nc*nx*sizeof(int),&colors);
270:       i1 = 0;
271:       for (i=xs; i<xs+nx; i++) {
272:         for (l=0; l<nc; l++) {
273:           colors[i1++] = l + nc*(i % col);
274:         }
275:       }
276:       ISColoringCreate(comm,nc*nx,colors,&da->localcoloring);
277:     }
278:     *coloring = da->localcoloring;
279:   } else if (ctype == IS_COLORING_GHOSTED) {
280:     if (!da->ghostedcoloring) {
281:       PetscMalloc(nc*gnx*sizeof(int),&colors);
282:       i1 = 0;
283:       for (i=gxs; i<gxs+gnx; i++) {
284:         for (l=0; l<nc; l++) {
285:           /* the complicated stuff is to handle periodic boundaries */
286:           colors[i1++] = l + nc*((((i < 0) ? m+i:((i >= m) ? i-m:i)) % col));
287:         }
288:       }
289:       ISColoringCreate(comm,nc*gnx,colors,&da->ghostedcoloring);
290:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
291:     }
292:     *coloring = da->ghostedcoloring;
293:   } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
294:   ISColoringReference(*coloring);
295:   return(0);
296: }

298: #undef __FUNCT__  
300: int DAGetColoring2d_5pt_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
301: {
302:   int      ierr,xs,ys,nx,ny,*colors,i,j,ii,gxs,gys,gnx,gny;
303:   int      m,n,dim,w,s,k,nc;
304:   MPI_Comm comm;

307:   /*     
308:          nc - number of components per grid point 
309:          col - number of colors needed in one direction for single component problem
310:   
311:   */
312:   ierr   = DAGetInfo(da,&dim,&m,&n,0,0,0,0,&w,&s,0,0);
313:   nc     = w;
314:   ierr   = DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
315:   ierr   = DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
316:   ierr   = PetscObjectGetComm((PetscObject)da,&comm);

318:   /* create the coloring */
319:   if (ctype == IS_COLORING_LOCAL) {
320:     if (!da->localcoloring) {
321:       PetscMalloc(nc*nx*ny*sizeof(int),&colors);
322:       ii = 0;
323:       for (j=ys; j<ys+ny; j++) {
324:         for (i=xs; i<xs+nx; i++) {
325:           for (k=0; k<nc; k++) {
326:             colors[ii++] = k + nc*((3*j+i) % 5);
327:           }
328:         }
329:       }
330:       ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
331:     }
332:     *coloring = da->localcoloring;
333:   } else if (ctype == IS_COLORING_GHOSTED) {
334:     if (!da->ghostedcoloring) {
335:       PetscMalloc(nc*gnx*gny*sizeof(int),&colors);
336:       ii = 0;
337:       for (j=gys; j<gys+gny; j++) {
338:         for (i=gxs; i<gxs+gnx; i++) {
339:           for (k=0; k<nc; k++) {
340:             colors[ii++] = k + nc*((3*j+i) % 5);
341:           }
342:         }
343:       }
344:       ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
345:       ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
346:     }
347:     *coloring = da->ghostedcoloring;
348:   } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
349:   return(0);
350: }

352: /* =========================================================================== */
353: EXTERN int DAGetMatrix1d_MPIAIJ(DA,Mat *);
354: EXTERN int DAGetMatrix2d_MPIAIJ(DA,Mat *);
355: EXTERN int DAGetMatrix3d_MPIAIJ(DA,Mat *);
356: EXTERN int DAGetMatrix3d_MPIBAIJ(DA,Mat *);
357: EXTERN int DAGetMatrix3d_MPISBAIJ(DA,Mat *);

359: #undef __FUNCT__  
361: /*@C
362:     DAGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for 
363:       computing the Jacobian on a function defined using the stencil set in the DA.

365:     Collective on DA

367:     Input Parameter:
368: +   da - the distributed array
369: -   mtype - either MATMPIAIJ or MATMPIBAIJ

371:     Output Parameters:
372: .   J  - matrix with the correct nonzero structure
373:         (obviously without the correct Jacobian values)

375:     Level: advanced

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

380: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate()

382: @*/
383: int DAGetMatrix(DA da,MatType mtype,Mat *J)
384: {
385:   int        ierr,dim;
386:   PetscTruth aij,baij,sbaij;

389:   /*
390:                                   m
391:           ------------------------------------------------------
392:          |                                                     |
393:          |                                                     |
394:          |               ----------------------                |
395:          |               |                    |                |
396:       n  |           yn  |                    |                |
397:          |               |                    |                |
398:          |               .---------------------                |
399:          |             (xs,ys)     xn                          |
400:          |            .                                        |
401:          |         (gxs,gys)                                   |
402:          |                                                     |
403:           -----------------------------------------------------
404:   */

406:   /*     
407:          nc - number of components per grid point 
408:          col - number of colors needed in one direction for single component problem
409:   
410:   */
411:   DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);

413:   /*
414:      We do not provide a getcoloring function in the DA operations because 
415:    the basic DA does not know about matrices. We think of DA as being more 
416:    more low-level then matrices.
417:   */
418:   PetscStrcmp(MATMPIAIJ,mtype,&aij);
419:   PetscStrcmp(MATMPIBAIJ,mtype,&baij);
420:   PetscStrcmp(MATMPISBAIJ,mtype,&sbaij);
421:   if (aij) {
422:     if (dim == 1) {
423:       DAGetMatrix1d_MPIAIJ(da,J);
424:     } else if (dim == 2) {
425:        DAGetMatrix2d_MPIAIJ(da,J);
426:     } else if (dim == 3) {
427:        DAGetMatrix3d_MPIAIJ(da,J);
428:     }
429:   } else if(baij) {
430:     if (dim == 3) {
431:        DAGetMatrix3d_MPIBAIJ(da,J);
432:     } else {
433:       SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
434:     }
435:   } else if(sbaij) {
436:     if (dim == 3) {
437:        DAGetMatrix3d_MPISBAIJ(da,J);
438:     } else {
439:       SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
440:     }
441:   } else {
442:     SETERRQ1(1,"Not done for %s matrix type, send us mail petsc-maint@mcs.anl.gov for code",mtype);
443:   }
444:   return(0);
445: }

447: /* ---------------------------------------------------------------------------------*/

449: #undef __FUNCT__  
451: int DAGetMatrix2d_MPIAIJ(DA da,Mat *J)
452: {
453:   int                    ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
454:   int                    m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
455:   int                    lstart,lend,pstart,pend,*dnz,*onz,size;
456:   int                    dims[2],starts[2];
457:   MPI_Comm               comm;
458:   PetscScalar            *values;
459:   DAPeriodicType         wrap;
460:   ISLocalToGlobalMapping ltog;
461:   DAStencilType          st;

464:   /*     
465:          nc - number of components per grid point 
466:          col - number of colors needed in one direction for single component problem
467:   
468:   */
469:   DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
470:   col = 2*s + 1;
471:   if (DAXPeriodic(wrap) && (m % col)){
472:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
473:                  by 2*stencil_width + 1n");
474:   }
475:   if (DAYPeriodic(wrap) && (n % col)){
476:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
477:                  by 2*stencil_width + 1n");
478:   }
479:   DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
480:   DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
481:   PetscObjectGetComm((PetscObject)da,&comm);
482:   MPI_Comm_size(comm,&size);

484:   /* create empty Jacobian matrix */
485:   MatCreate(comm,nc*nx*ny,nc*nx*ny,PETSC_DECIDE,PETSC_DECIDE,J);

487:   PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
488:   PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
489:   PetscMalloc(nc*sizeof(int),&rows);
490:   PetscMalloc(col*col*nc*nc*sizeof(int),&cols);
491:   DAGetISLocalToGlobalMapping(da,&ltog);
492: 
493:   /* determine the matrix preallocation information */
494:   MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
495:   for (i=xs; i<xs+nx; i++) {

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

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

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

506:       cnt  = 0;
507:       for (k=0; k<nc; k++) {
508:         for (l=lstart; l<lend+1; l++) {
509:           for (p=pstart; p<pend+1; p++) {
510:             if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
511:               cols[cnt++]  = k + nc*(slot + gnx*l + p);
512:             }
513:           }
514:         }
515:         rows[k] = k + nc*(slot);
516:       }
517:       MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
518:     }
519:   }
520:   /* set matrix type and preallocation information */
521:   if (size > 1) {
522:     MatSetType(*J,MATMPIAIJ);
523:   } else {
524:     MatSetType(*J,MATSEQAIJ);
525:   }
526:   MatSeqAIJSetPreallocation(*J,0,dnz);
527:   MatSeqBAIJSetPreallocation(*J,nc,0,dnz);
528:   MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
529:   MatMPIBAIJSetPreallocation(*J,nc,0,dnz,0,onz);
530:   MatPreallocateFinalize(dnz,onz);
531:   MatSetLocalToGlobalMapping(*J,ltog);
532:   DAGetGhostCorners(da,&starts[0],&starts[1],PETSC_IGNORE,&dims[0],&dims[1],PETSC_IGNORE);
533:   MatSetStencil(*J,2,dims,starts,nc);

535:   /*
536:     For each node in the grid: we get the neighbors in the local (on processor ordering
537:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
538:     PETSc ordering.
539:   */
540:   for (i=xs; i<xs+nx; i++) {
541: 
542:     pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
543:     pend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
544: 
545:     for (j=ys; j<ys+ny; j++) {
546:       slot = i - gxs + gnx*(j - gys);
547: 
548:       lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
549:       lend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));

551:       cnt  = 0;
552:       for (k=0; k<nc; k++) {
553:         for (l=lstart; l<lend+1; l++) {
554:           for (p=pstart; p<pend+1; p++) {
555:             if ((st == DA_STENCIL_BOX) || (!l || !p)) {  /* entries on star have either l = 0 or p = 0 */
556:               cols[cnt++]  = k + nc*(slot + gnx*l + p);
557:             }
558:           }
559:         }
560:         rows[k]      = k + nc*(slot);
561:       }
562:       MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
563:     }
564:   }
565:   PetscFree(values);
566:   PetscFree(rows);
567:   PetscFree(cols);
568:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
569:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
570:   return(0);
571: }

573: /* ---------------------------------------------------------------------------------*/

575: #undef __FUNCT__  
577: int DAGetMatrix3d_MPIAIJ(DA da,Mat *J)
578: {
579:   int                    ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
580:   int                    m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
581:   int                    istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,size;
582:   int                    dims[3],starts[3];
583:   MPI_Comm               comm;
584:   PetscScalar            *values;
585:   DAPeriodicType         wrap;
586:   ISLocalToGlobalMapping ltog;
587:   DAStencilType          st;

590:   /*     
591:          nc - number of components per grid point 
592:          col - number of colors needed in one direction for single component problem
593:   
594:   */
595:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
596:   col    = 2*s + 1;
597:   if (DAXPeriodic(wrap) && (m % col)){
598:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
599:                  by 2*stencil_width + 1n");
600:   }
601:   if (DAYPeriodic(wrap) && (n % col)){
602:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
603:                  by 2*stencil_width + 1n");
604:   }
605:   if (DAZPeriodic(wrap) && (p % col)){
606:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisiblen
607:                  by 2*stencil_width + 1n");
608:   }

610:   DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
611:   DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
612:   PetscObjectGetComm((PetscObject)da,&comm);
613:   MPI_Comm_size(comm,&size);


616:   /* create the matrix */
617:   /* create empty Jacobian matrix */
618:   MatCreate(comm,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,J);
619:   PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
620:   PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
621:   PetscMalloc(nc*sizeof(int),&rows);
622:   PetscMalloc(col*col*col*nc*sizeof(int),&cols);
623:   DAGetISLocalToGlobalMapping(da,&ltog);

625:   /* determine the matrix preallocation information */
626:   MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
627:   for (i=xs; i<xs+nx; i++) {
628:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
629:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
630:     for (j=ys; j<ys+ny; j++) {
631:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
632:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
633:       for (k=zs; k<zs+nz; k++) {
634:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
635:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
636: 
637:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
638: 
639:         cnt  = 0;
640:         for (l=0; l<nc; l++) {
641:           for (ii=istart; ii<iend+1; ii++) {
642:             for (jj=jstart; jj<jend+1; jj++) {
643:               for (kk=kstart; kk<kend+1; kk++) {
644:                 if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
645:                   cols[cnt++]  = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
646:                 }
647:               }
648:             }
649:           }
650:           rows[l] = l + nc*(slot);
651:         }
652:         MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
653:       }
654:     }
655:   }
656:   /* set matrix type and preallocation */
657:   if (size > 1) {
658:     MatSetType(*J,MATMPIAIJ);
659:   } else {
660:     MatSetType(*J,MATSEQAIJ);
661:   }
662:   MatSeqAIJSetPreallocation(*J,0,dnz);
663:   MatSeqBAIJSetPreallocation(*J,nc,0,dnz);
664:   MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
665:   MatMPIBAIJSetPreallocation(*J,nc,0,dnz,0,onz);
666:   MatPreallocateFinalize(dnz,onz);
667:   MatSetLocalToGlobalMapping(*J,ltog);
668:   DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
669:   MatSetStencil(*J,3,dims,starts,nc);

671:   /*
672:     For each node in the grid: we get the neighbors in the local (on processor ordering
673:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
674:     PETSc ordering.
675:   */
676:   for (i=xs; i<xs+nx; i++) {
677:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
678:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
679:     for (j=ys; j<ys+ny; j++) {
680:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
681:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
682:       for (k=zs; k<zs+nz; k++) {
683:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
684:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
685: 
686:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
687: 
688:         cnt  = 0;
689:         for (l=0; l<nc; l++) {
690:           for (ii=istart; ii<iend+1; ii++) {
691:             for (jj=jstart; jj<jend+1; jj++) {
692:               for (kk=kstart; kk<kend+1; kk++) {
693:                 if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
694:                   cols[cnt++]  = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
695:                 }
696:               }
697:             }
698:           }
699:           rows[l]      = l + nc*(slot);
700:         }
701:         MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
702:       }
703:     }
704:   }
705:   PetscFree(values);
706:   PetscFree(rows);
707:   PetscFree(cols);
708:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
709:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
710:   return(0);
711: }

713: /* ---------------------------------------------------------------------------------*/

715: #undef __FUNCT__  
717: int DAGetMatrix1d_MPIAIJ(DA da,Mat *J)
718: {
719:   int                    ierr,xs,nx,i,i1,slot,gxs,gnx;
720:   int                    m,dim,s,*cols,nc,*rows,col,cnt,l;
721:   int                    istart,iend,size;
722:   int                    dims[1],starts[1];
723:   MPI_Comm               comm;
724:   PetscScalar            *values;
725:   DAPeriodicType         wrap;
726:   ISLocalToGlobalMapping ltog;

729:   /*     
730:          nc - number of components per grid point 
731:          col - number of colors needed in one direction for single component problem
732:   
733:   */
734:   DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
735:   col    = 2*s + 1;

737:   if (DAXPeriodic(wrap) && (m % col)) {
738:     SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisiblen
739:                  by 2*stencil_width + 1n");
740:   }


743:   DAGetCorners(da,&xs,0,0,&nx,0,0);
744:   DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
745:   PetscObjectGetComm((PetscObject)da,&comm);
746:   MPI_Comm_size(comm,&size);

748:   /* create empty Jacobian matrix */
749: 
750:   ierr    = MatCreate(comm,nc*nx,nc*nx,PETSC_DECIDE,PETSC_DECIDE,J);
751:   if (size > 1) {
752:     MatSetType(*J,MATMPIAIJ);
753:   } else {
754:     MatSetType(*J,MATSEQAIJ);
755:   }
756:   MatSeqAIJSetPreallocation(*J,col*nc,0);
757:   MatSeqBAIJSetPreallocation(*J,nc,col,0);
758:   MatMPIAIJSetPreallocation(*J,col*nc,0,0,0);
759:   MatMPIBAIJSetPreallocation(*J,nc,col,0,0,0);
760:   DAGetGhostCorners(da,&starts[0],PETSC_IGNORE,PETSC_IGNORE,&dims[0],PETSC_IGNORE,PETSC_IGNORE);
761:   MatSetStencil(*J,1,dims,starts,nc);
762: 
763:   PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
764:   PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
765:   PetscMalloc(nc*sizeof(int),&rows);
766:   PetscMalloc(col*nc*sizeof(int),&cols);
767: 
768:   DAGetISLocalToGlobalMapping(da,&ltog);
769:   MatSetLocalToGlobalMapping(*J,ltog);
770: 
771:   /*
772:     For each node in the grid: we get the neighbors in the local (on processor ordering
773:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
774:     PETSc ordering.
775:   */
776:   for (i=xs; i<xs+nx; i++) {
777:     istart = PetscMax(-s,gxs - i);
778:     iend   = PetscMin(s,gxs + gnx - i - 1);
779:     slot   = i - gxs;
780: 
781:     cnt  = 0;
782:     for (l=0; l<nc; l++) {
783:       for (i1=istart; i1<iend+1; i1++) {
784:         cols[cnt++] = l + nc*(slot + i1);
785:       }
786:       rows[l]      = l + nc*(slot);
787:     }
788:     MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
789:   }
790:   PetscFree(values);
791:   PetscFree(rows);
792:   PetscFree(cols);
793:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
794:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
795:   return(0);
796: }

798: #undef __FUNCT__  
800: int DAGetMatrix3d_MPIBAIJ(DA da,Mat *J)
801: {
802:   int                    ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
803:   int                    m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
804:   int                    istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
805:   MPI_Comm               comm;
806:   PetscScalar            *values;
807:   DAPeriodicType         wrap;
808:   DAStencilType          st;
809:   ISLocalToGlobalMapping ltog;

812:   /*     
813:          nc - number of components per grid point 
814:          col - number of colors needed in one direction for single component problem
815:   
816:   */
817:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
818:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
819:   col    = 2*s + 1;

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

825:   /* create the matrix */
826:   ierr  = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
827:   ierr  = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
828:   ierr  = PetscMalloc(col*col*col*sizeof(int),&cols);

830:   DAGetISLocalToGlobalMappingBlck(da,&ltog);

832:   /* determine the matrix preallocation information */
833:   MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
834:   for (i=xs; i<xs+nx; i++) {
835:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
836:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
837:     for (j=ys; j<ys+ny; j++) {
838:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
839:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
840:       for (k=zs; k<zs+nz; k++) {
841:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
842:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));

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

846:         /* Find block columns in block row */
847:         cnt  = 0;
848:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
849:           for (ii=istart; ii<iend+1; ii++) {
850:             for (jj=jstart; jj<jend+1; jj++) {
851:               for (kk=kstart; kk<kend+1; kk++) {
852:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
853:               }
854:             }
855:           }
856:         } else {  /* Star stencil */
857:           cnt  = 0;
858:           for (ii=istart; ii<iend+1; ii++) {
859:             if (ii) {
860:               /* jj and kk must be zero */
861:               /* cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk; */
862:               cols[cnt++]  = slot + ii;
863:             } else {
864:               for (jj=jstart; jj<jend+1; jj++) {
865:                 if (jj) {
866:                   /* ii and kk must be zero */
867:                   cols[cnt++]  = slot + gnx*jj;
868:                 } else {
869:                   /* ii and jj must be zero */
870:                   for (kk=kstart; kk<kend+1; kk++) {
871:                     cols[cnt++]  = slot + gnx*gny*kk;
872:                   }
873:                 }
874:               }
875:             }
876:           }
877:         }
878:         MatPreallocateSetLocal(ltog,1,&slot,cnt,cols,dnz,onz);
879:       }
880:     }
881:   }

883:   /* create empty Jacobian matrix */
884:   MatCreateMPIBAIJ(comm,nc,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,0,dnz,0,onz,J);

886:   MatPreallocateFinalize(dnz,onz);
887:   MatSetLocalToGlobalMappingBlock(*J,ltog);

889:   /*
890:     For each node in the grid: we get the neighbors in the local (on processor ordering
891:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
892:     PETSc ordering.
893:   */

895:   for (i=xs; i<xs+nx; i++) {
896:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
897:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
898:     for (j=ys; j<ys+ny; j++) {
899:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
900:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
901:       for (k=zs; k<zs+nz; k++) {
902:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
903:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
904: 
905:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
906: 
907:         cnt  = 0;
908:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
909:           for (ii=istart; ii<iend+1; ii++) {
910:             for (jj=jstart; jj<jend+1; jj++) {
911:               for (kk=kstart; kk<kend+1; kk++) {
912:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
913:               }
914:             }
915:           }
916:         } else {  /* Star stencil */
917:           cnt  = 0;
918:           for (ii=istart; ii<iend+1; ii++) {
919:             if (ii) {
920:               /* jj and kk must be zero */
921:               cols[cnt++]  = slot + ii;
922:             } else {
923:               for (jj=jstart; jj<jend+1; jj++) {
924:                 if (jj) {
925:                   /* ii and kk must be zero */
926:                   cols[cnt++]  = slot + gnx*jj;
927:                 } else {
928:                   /* ii and jj must be zero */
929:                   for (kk=kstart; kk<kend+1; kk++) {
930:                     cols[cnt++]  = slot + gnx*gny*kk;
931:                   }
932:                 }
933:               }
934:             }
935:           }
936:         }
937:         MatSetValuesBlockedLocal(*J,1,&slot,cnt,cols,values,INSERT_VALUES);
938:       }
939:     }
940:   }
941:   PetscFree(values);
942:   PetscFree(cols);
943:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
944:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
945:   return(0);
946: }

948: /* BAD! Almost identical to the BAIJ one */
949: #undef __FUNCT__  
951: int DAGetMatrix3d_MPISBAIJ(DA da,Mat *J)
952: {
953:   int                    ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
954:   int                    m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
955:   int                    istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
956:   MPI_Comm               comm;
957:   PetscScalar            *values;
958:   DAPeriodicType         wrap;
959:   DAStencilType          st;
960:   ISLocalToGlobalMapping ltog;

963:   /*     
964:          nc - number of components per grid point 
965:          col - number of colors needed in one direction for single component problem
966:   
967:   */
968:   DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
969:   if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
970:   col    = 2*s + 1;

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

976:   /* create the matrix */
977:   ierr  = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
978:   ierr  = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
979:   ierr  = PetscMalloc(col*col*col*sizeof(int),&cols);

981:   DAGetISLocalToGlobalMappingBlck(da,&ltog);

983:   /* determine the matrix preallocation information */
984:   MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
985:   for (i=xs; i<xs+nx; i++) {
986:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
987:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
988:     for (j=ys; j<ys+ny; j++) {
989:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
990:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
991:       for (k=zs; k<zs+nz; k++) {
992:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
993:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));

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

997:         /* Find block columns in block row */
998:         cnt  = 0;
999:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
1000:           for (ii=istart; ii<iend+1; ii++) {
1001:             for (jj=jstart; jj<jend+1; jj++) {
1002:               for (kk=kstart; kk<kend+1; kk++) {
1003:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
1004:               }
1005:             }
1006:           }
1007:         } else {  /* Star stencil */
1008:           cnt  = 0;
1009:           for (ii=istart; ii<iend+1; ii++) {
1010:             if (ii) {
1011:               /* jj and kk must be zero */
1012:               /* cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk; */
1013:               cols[cnt++]  = slot + ii;
1014:             } else {
1015:               for (jj=jstart; jj<jend+1; jj++) {
1016:                 if (jj) {
1017:                   /* ii and kk must be zero */
1018:                   cols[cnt++]  = slot + gnx*jj;
1019:                 } else {
1020:                   /* ii and jj must be zero */
1021:                   for (kk=kstart; kk<kend+1; kk++) {
1022:                     cols[cnt++]  = slot + gnx*gny*kk;
1023:                   }
1024:                 }
1025:               }
1026:             }
1027:           }
1028:         }
1029:         MatPreallocateSymmetricSetLocal(ltog,1,&slot,cnt,cols,dnz,onz);
1030:       }
1031:     }
1032:   }

1034:   /* create empty Jacobian matrix */
1035:   MatCreateMPISBAIJ(comm,nc,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,0,dnz,0,onz,J);

1037:   MatPreallocateFinalize(dnz,onz);
1038:   MatSetLocalToGlobalMappingBlock(*J,ltog);

1040:   /*
1041:     For each node in the grid: we get the neighbors in the local (on processor ordering
1042:     that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1043:     PETSc ordering.
1044:   */

1046:   for (i=xs; i<xs+nx; i++) {
1047:     istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1048:     iend   = DAXPeriodic(wrap) ?  s : (PetscMin(s,m-i-1));
1049:     for (j=ys; j<ys+ny; j++) {
1050:       jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1051:       jend   = DAYPeriodic(wrap) ?  s : (PetscMin(s,n-j-1));
1052:       for (k=zs; k<zs+nz; k++) {
1053:         kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1054:         kend   = DAZPeriodic(wrap) ?  s : (PetscMin(s,p-k-1));
1055: 
1056:         slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1057: 
1058:         cnt  = 0;
1059:         if (st == DA_STENCIL_BOX) {   /* if using BOX stencil */
1060:           for (ii=istart; ii<iend+1; ii++) {
1061:             for (jj=jstart; jj<jend+1; jj++) {
1062:               for (kk=kstart; kk<kend+1; kk++) {
1063:                 cols[cnt++]  = slot + ii + gnx*jj + gnx*gny*kk;
1064:               }
1065:             }
1066:           }
1067:         } else {  /* Star stencil */
1068:           cnt  = 0;
1069:           for (ii=istart; ii<iend+1; ii++) {
1070:             if (ii) {
1071:               /* jj and kk must be zero */
1072:               cols[cnt++]  = slot + ii;
1073:             } else {
1074:               for (jj=jstart; jj<jend+1; jj++) {
1075:                 if (jj) {
1076:                   /* ii and kk must be zero */
1077:                   cols[cnt++]  = slot + gnx*jj;
1078:                 } else {
1079:                   /* ii and jj must be zero */
1080:                   for (kk=kstart; kk<kend+1; kk++) {
1081:                     cols[cnt++]  = slot + gnx*gny*kk;
1082:                   }
1083:                 }
1084:               }
1085:             }
1086:           }
1087:         }
1088:         MatSetValuesBlockedLocal(*J,1,&slot,cnt,cols,values,INSERT_VALUES);
1089:       }
1090:     }
1091:   }
1092:   PetscFree(values);
1093:   PetscFree(cols);
1094:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1095:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1096:   return(0);
1097: }