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: /*@
54: DASetMatPreallocateOnly - When DAGetMatrix() is called the matrix will be properly
55: preallocated but the nonzero structure and zero values will not be set.
57: Collective on DA
59: Input Parameter:
60: + da - the distributed array
61: - only - PETSC_TRUE if only want preallocation
64: Level: developer
66: .seealso DAGetMatrix(), DASetGetMatrix(), DASetBlockSize(), DASetBlockFills()
68: @*/
69: PetscErrorCode PETSCDM_DLLEXPORT DASetMatPreallocateOnly(DA da,PetscTruth only)
70: {
72: da->prealloc_only = only;
73: return(0);
74: }
78: /*@
79: DASetBlockFills - Sets the fill pattern in each block for a multi-component problem
80: of the matrix returned by DAGetMatrix().
82: Collective on DA
84: Input Parameter:
85: + da - the distributed array
86: . dfill - the fill pattern in the diagonal block (may be PETSC_NULL, means use dense block)
87: - ofill - the fill pattern in the off-diagonal blocks
90: Level: developer
92: Notes: This only makes sense when you are doing multicomponent problems but using the
93: MPIAIJ matrix format
95: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
96: representing coupling and 0 entries for missing coupling. For example
97: $ dfill[3][3] = {1, 0, 0,
98: $ 1, 1, 0,
99: $ 0, 1, 1}
100: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
101: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
102: diagonal block.
104: DASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
105: can be represented in the dfill, ofill format
107: Contributed by Glenn Hammond
109: .seealso DAGetMatrix(), DASetGetMatrix(), DASetMatPreallocateOnly()
111: @*/
112: PetscErrorCode PETSCDM_DLLEXPORT DASetBlockFills(DA da,PetscInt *dfill,PetscInt *ofill)
113: {
117: if (dfill) {
118: DASetBlockFills_Private(dfill,da->w,&da->dfill);
119: }
120: DASetBlockFills_Private(ofill,da->w,&da->ofill);
121: return(0);
122: }
127: /*@
128: DAGetColoring - Gets the coloring required for computing the Jacobian via
129: finite differences on a function defined using a stencil on the DA.
131: Collective on DA
133: Input Parameter:
134: + da - the distributed array
135: - ctype - IS_COLORING_LOCAL or IS_COLORING_GHOSTED
137: Output Parameters:
138: . coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)
140: Level: advanced
142: Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used
143: for efficient (parallel or thread based) triangular solves etc is NOT yet
144: available.
147: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring
149: @*/
150: PetscErrorCode PETSCDM_DLLEXPORT DAGetColoring(DA da,ISColoringType ctype,ISColoring *coloring)
151: {
153: PetscInt dim;
156: /*
157: m
158: ------------------------------------------------------
159: | |
160: | |
161: | ---------------------- |
162: | | | |
163: n | yn | | |
164: | | | |
165: | .--------------------- |
166: | (xs,ys) xn |
167: | . |
168: | (gxs,gys) |
169: | |
170: -----------------------------------------------------
171: */
173: /*
174: nc - number of components per grid point
175: col - number of colors needed in one direction for single component problem
176:
177: */
178: DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
180: /*
181: We do not provide a getcoloring function in the DA operations because
182: the basic DA does not know about matrices. We think of DA as being more
183: more low-level then matrices.
184: */
185: if (dim == 1) {
186: DAGetColoring1d_MPIAIJ(da,ctype,coloring);
187: } else if (dim == 2) {
188: DAGetColoring2d_MPIAIJ(da,ctype,coloring);
189: } else if (dim == 3) {
190: DAGetColoring3d_MPIAIJ(da,ctype,coloring);
191: } else {
192: SETERRQ1(PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
193: }
194: return(0);
195: }
197: /* ---------------------------------------------------------------------------------*/
201: PetscErrorCode DAGetColoring2d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
202: {
203: PetscErrorCode ierr;
204: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
205: PetscInt ncolors;
206: MPI_Comm comm;
207: DAPeriodicType wrap;
208: DAStencilType st;
209: ISColoringValue *colors;
212: /*
213: nc - number of components per grid point
214: col - number of colors needed in one direction for single component problem
215:
216: */
217: DAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&wrap,&st);
218: col = 2*s + 1;
219: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
220: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
221: PetscObjectGetComm((PetscObject)da,&comm);
223: /* special case as taught to us by Paul Hovland */
224: if (st == DA_STENCIL_STAR && s == 1) {
225: DAGetColoring2d_5pt_MPIAIJ(da,ctype,coloring);
226: } else {
228: if (DAXPeriodic(wrap) && (m % col)){
229: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
230: by 2*stencil_width + 1\n");
231: }
232: if (DAYPeriodic(wrap) && (n % col)){
233: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
234: by 2*stencil_width + 1\n");
235: }
236: if (ctype == IS_COLORING_LOCAL) {
237: if (!da->localcoloring) {
238: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
239: ii = 0;
240: for (j=ys; j<ys+ny; j++) {
241: for (i=xs; i<xs+nx; i++) {
242: for (k=0; k<nc; k++) {
243: colors[ii++] = k + nc*((i % col) + col*(j % col));
244: }
245: }
246: }
247: ncolors = nc + nc*(col-1 + col*(col-1));
248: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&da->localcoloring);
249: }
250: *coloring = da->localcoloring;
251: } else if (ctype == IS_COLORING_GHOSTED) {
252: if (!da->ghostedcoloring) {
253: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
254: ii = 0;
255: for (j=gys; j<gys+gny; j++) {
256: for (i=gxs; i<gxs+gnx; i++) {
257: for (k=0; k<nc; k++) {
258: /* the complicated stuff is to handle periodic boundaries */
259: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
260: }
261: }
262: }
263: ncolors = nc + nc*(col - 1 + col*(col-1));
264: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&da->ghostedcoloring);
265: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
266: }
267: *coloring = da->ghostedcoloring;
268: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
269: }
270: ISColoringReference(*coloring);
271: return(0);
272: }
274: /* ---------------------------------------------------------------------------------*/
278: PetscErrorCode DAGetColoring3d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
279: {
280: PetscErrorCode ierr;
281: 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;
282: PetscInt ncolors;
283: MPI_Comm comm;
284: DAPeriodicType wrap;
285: DAStencilType st;
286: ISColoringValue *colors;
289: /*
290: nc - number of components per grid point
291: col - number of colors needed in one direction for single component problem
292:
293: */
294: DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&wrap,&st);
295: col = 2*s + 1;
296: if (DAXPeriodic(wrap) && (m % col)){
297: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
298: by 2*stencil_width + 1\n");
299: }
300: if (DAYPeriodic(wrap) && (n % col)){
301: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
302: by 2*stencil_width + 1\n");
303: }
304: if (DAZPeriodic(wrap) && (p % col)){
305: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
306: by 2*stencil_width + 1\n");
307: }
309: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
310: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
311: PetscObjectGetComm((PetscObject)da,&comm);
313: /* create the coloring */
314: if (ctype == IS_COLORING_LOCAL) {
315: if (!da->localcoloring) {
316: PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);
317: ii = 0;
318: for (k=zs; k<zs+nz; k++) {
319: for (j=ys; j<ys+ny; j++) {
320: for (i=xs; i<xs+nx; i++) {
321: for (l=0; l<nc; l++) {
322: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
323: }
324: }
325: }
326: }
327: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
328: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,&da->localcoloring);
329: }
330: *coloring = da->localcoloring;
331: } else if (ctype == IS_COLORING_GHOSTED) {
332: if (!da->ghostedcoloring) {
333: PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);
334: ii = 0;
335: for (k=gzs; k<gzs+gnz; k++) {
336: for (j=gys; j<gys+gny; j++) {
337: for (i=gxs; i<gxs+gnx; i++) {
338: for (l=0; l<nc; l++) {
339: /* the complicated stuff is to handle periodic boundaries */
340: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
341: }
342: }
343: }
344: }
345: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
346: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,&da->ghostedcoloring);
347: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
348: }
349: *coloring = da->ghostedcoloring;
350: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
351: ISColoringReference(*coloring);
352: return(0);
353: }
355: /* ---------------------------------------------------------------------------------*/
359: PetscErrorCode DAGetColoring1d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
360: {
361: PetscErrorCode ierr;
362: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
363: PetscInt ncolors;
364: MPI_Comm comm;
365: DAPeriodicType wrap;
366: ISColoringValue *colors;
369: /*
370: nc - number of components per grid point
371: col - number of colors needed in one direction for single component problem
372:
373: */
374: DAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&wrap,0);
375: col = 2*s + 1;
377: if (DAXPeriodic(wrap) && (m % col)) {
378: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisible\n\
379: by 2*stencil_width + 1\n");
380: }
382: DAGetCorners(da,&xs,0,0,&nx,0,0);
383: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
384: PetscObjectGetComm((PetscObject)da,&comm);
386: /* create the coloring */
387: if (ctype == IS_COLORING_LOCAL) {
388: if (!da->localcoloring) {
389: PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);
390: i1 = 0;
391: for (i=xs; i<xs+nx; i++) {
392: for (l=0; l<nc; l++) {
393: colors[i1++] = l + nc*(i % col);
394: }
395: }
396: ncolors = nc + nc*(col-1);
397: ISColoringCreate(comm,ncolors,nc*nx,colors,&da->localcoloring);
398: }
399: *coloring = da->localcoloring;
400: } else if (ctype == IS_COLORING_GHOSTED) {
401: if (!da->ghostedcoloring) {
402: PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);
403: i1 = 0;
404: for (i=gxs; i<gxs+gnx; i++) {
405: for (l=0; l<nc; l++) {
406: /* the complicated stuff is to handle periodic boundaries */
407: colors[i1++] = l + nc*(SetInRange(i,m) % col);
408: }
409: }
410: ncolors = nc + nc*(col-1);
411: ISColoringCreate(comm,ncolors,nc*gnx,colors,&da->ghostedcoloring);
412: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
413: }
414: *coloring = da->ghostedcoloring;
415: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
416: ISColoringReference(*coloring);
417: return(0);
418: }
422: PetscErrorCode DAGetColoring2d_5pt_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
423: {
424: PetscErrorCode ierr;
425: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
426: PetscInt ncolors;
427: MPI_Comm comm;
428: DAPeriodicType wrap;
429: ISColoringValue *colors;
432: /*
433: nc - number of components per grid point
434: col - number of colors needed in one direction for single component problem
435:
436: */
437: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,0);
438: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
439: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
440: PetscObjectGetComm((PetscObject)da,&comm);
442: if (DAXPeriodic(wrap) && (m % 5)){
443: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
444: by 5\n");
445: }
446: if (DAYPeriodic(wrap) && (n % 5)){
447: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
448: by 5\n");
449: }
451: /* create the coloring */
452: if (ctype == IS_COLORING_LOCAL) {
453: if (!da->localcoloring) {
454: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
455: ii = 0;
456: for (j=ys; j<ys+ny; j++) {
457: for (i=xs; i<xs+nx; i++) {
458: for (k=0; k<nc; k++) {
459: colors[ii++] = k + nc*((3*j+i) % 5);
460: }
461: }
462: }
463: ncolors = 5*nc;
464: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&da->localcoloring);
465: }
466: *coloring = da->localcoloring;
467: } else if (ctype == IS_COLORING_GHOSTED) {
468: if (!da->ghostedcoloring) {
469: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
470: ii = 0;
471: for (j=gys; j<gys+gny; j++) {
472: for (i=gxs; i<gxs+gnx; i++) {
473: for (k=0; k<nc; k++) {
474: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
475: }
476: }
477: }
478: ncolors = 5*nc;
479: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&da->ghostedcoloring);
480: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
481: }
482: *coloring = da->ghostedcoloring;
483: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
484: return(0);
485: }
487: /* =========================================================================== */
488: EXTERN PetscErrorCode DAGetMatrix1d_MPIAIJ(DA,Mat);
489: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ(DA,Mat);
490: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA,Mat);
491: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ(DA,Mat);
492: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA,Mat);
493: EXTERN PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA,Mat);
494: EXTERN PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA,Mat);
495: EXTERN PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA,Mat);
496: EXTERN PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA,Mat);
500: /*@C
501: DAGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for
502: computing the Jacobian on a function defined using the stencil set in the DA.
504: Collective on DA
506: Input Parameter:
507: + da - the distributed array
508: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
509: or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
511: Output Parameters:
512: . J - matrix with the correct nonzero structure
513: (obviously without the correct Jacobian values)
515: Level: advanced
517: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
518: do not need to do it yourself.
520: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
521: the nonzero pattern call DASetMatPreallocateOnly()
523: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills(), DASetMatPreallocateOnly()
525: @*/
526: PetscErrorCode PETSCDM_DLLEXPORT DAGetMatrix(DA da, MatType mtype,Mat *J)
527: {
529: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3];
530: Mat A;
531: MPI_Comm comm;
532: MatType Atype;
533: void (*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL;
536: /*
537: m
538: ------------------------------------------------------
539: | |
540: | |
541: | ---------------------- |
542: | | | |
543: n | ny | | |
544: | | | |
545: | .--------------------- |
546: | (xs,ys) nx |
547: | . |
548: | (gxs,gys) |
549: | |
550: -----------------------------------------------------
551: */
553: /*
554: nc - number of components per grid point
555: col - number of colors needed in one direction for single component problem
556:
557: */
558: DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0);
559: DAGetCorners(da,0,0,0,&nx,&ny,&nz);
560: PetscObjectGetComm((PetscObject)da,&comm);
561: MatCreate(comm,&A);
562: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
563: MatSetType(A,mtype);
564: MatSetFromOptions(A);
565: MatGetType(A,&Atype);
566: /*
567: We do not provide a getmatrix function in the DA operations because
568: the basic DA does not know about matrices. We think of DA as being more
569: more low-level than matrices. This is kind of cheating but, cause sometimes
570: we think of DA has higher level than matrices.
572: We could switch based on Atype (or mtype), but we do not since the
573: specialized setting routines depend only the particular preallocation
574: details of the matrix, not the type itself.
575: */
576: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
577: if (!aij) {
578: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
579: }
580: if (!aij) {
581: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
582: if (!baij) {
583: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
584: }
585: if (!baij){
586: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
587: if (!sbaij) {
588: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
589: }
590: if (!sbaij) SETERRQ2(PETSC_ERR_SUP,"Not implemented for the matrix type: %s!\n" \
591: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
592: }
593: }
594: if (aij) {
595: if (dim == 1) {
596: DAGetMatrix1d_MPIAIJ(da,A);
597: } else if (dim == 2) {
598: if (da->ofill) {
599: DAGetMatrix2d_MPIAIJ_Fill(da,A);
600: } else {
601: DAGetMatrix2d_MPIAIJ(da,A);
602: }
603: } else if (dim == 3) {
604: if (da->ofill) {
605: DAGetMatrix3d_MPIAIJ_Fill(da,A);
606: } else {
607: DAGetMatrix3d_MPIAIJ(da,A);
608: }
609: }
610: } else if (baij) {
611: if (dim == 2) {
612: DAGetMatrix2d_MPIBAIJ(da,A);
613: } else if (dim == 3) {
614: DAGetMatrix3d_MPIBAIJ(da,A);
615: } else {
616: SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
617: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
618: }
619: } else if (sbaij) {
620: if (dim == 2) {
621: DAGetMatrix2d_MPISBAIJ(da,A);
622: } else if (dim == 3) {
623: DAGetMatrix3d_MPISBAIJ(da,A);
624: } else {
625: SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
626: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
627: }
628: }
629: DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
630: MatSetStencil(A,dim,dims,starts,dof);
631: *J = A;
632: return(0);
633: }
635: /* ---------------------------------------------------------------------------------*/
638: PetscErrorCode DAGetMatrix2d_MPIAIJ(DA da,Mat J)
639: {
640: PetscErrorCode ierr;
641: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
642: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
643: MPI_Comm comm;
644: PetscScalar *values;
645: DAPeriodicType wrap;
646: ISLocalToGlobalMapping ltog;
647: DAStencilType st;
650: /*
651: nc - number of components per grid point
652: col - number of colors needed in one direction for single component problem
653:
654: */
655: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
656: col = 2*s + 1;
657: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
658: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
659: PetscObjectGetComm((PetscObject)da,&comm);
661: PetscMalloc2(nc,PetscInt,&rows,col*col*nc*nc,PetscInt,&cols);
662: DAGetISLocalToGlobalMapping(da,<og);
663:
664: /* determine the matrix preallocation information */
665: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
666: for (i=xs; i<xs+nx; i++) {
668: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
669: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
671: for (j=ys; j<ys+ny; j++) {
672: slot = i - gxs + gnx*(j - gys);
674: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
675: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
677: cnt = 0;
678: for (k=0; k<nc; k++) {
679: for (l=lstart; l<lend+1; l++) {
680: for (p=pstart; p<pend+1; p++) {
681: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
682: cols[cnt++] = k + nc*(slot + gnx*l + p);
683: }
684: }
685: }
686: rows[k] = k + nc*(slot);
687: }
688: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
689: }
690: }
691: MatSeqAIJSetPreallocation(J,0,dnz);
692: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
693: MatSetBlockSize(J,nc);
694: MatPreallocateFinalize(dnz,onz);
696: MatSetLocalToGlobalMapping(J,ltog);
698: /*
699: For each node in the grid: we get the neighbors in the local (on processor ordering
700: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
701: PETSc ordering.
702: */
703: if (!da->prealloc_only) {
704: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
705: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
706: for (i=xs; i<xs+nx; i++) {
707:
708: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
709: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
710:
711: for (j=ys; j<ys+ny; j++) {
712: slot = i - gxs + gnx*(j - gys);
713:
714: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
715: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
717: cnt = 0;
718: for (k=0; k<nc; k++) {
719: for (l=lstart; l<lend+1; l++) {
720: for (p=pstart; p<pend+1; p++) {
721: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
722: cols[cnt++] = k + nc*(slot + gnx*l + p);
723: }
724: }
725: }
726: rows[k] = k + nc*(slot);
727: }
728: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
729: }
730: }
731: PetscFree(values);
732: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
733: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
734: }
735: PetscFree2(rows,cols);
736: return(0);
737: }
741: PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA da,Mat J)
742: {
743: PetscErrorCode ierr;
744: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
745: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p;
746: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
747: PetscInt ifill_col,*ofill = da->ofill, *dfill = da->dfill;
748: MPI_Comm comm;
749: PetscScalar *values;
750: DAPeriodicType wrap;
751: ISLocalToGlobalMapping ltog;
752: DAStencilType st;
755: /*
756: nc - number of components per grid point
757: col - number of colors needed in one direction for single component problem
758:
759: */
760: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
761: col = 2*s + 1;
762: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
763: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
764: PetscObjectGetComm((PetscObject)da,&comm);
766: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
767: DAGetISLocalToGlobalMapping(da,<og);
768:
769: /* determine the matrix preallocation information */
770: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
771: for (i=xs; i<xs+nx; i++) {
773: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
774: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
776: for (j=ys; j<ys+ny; j++) {
777: slot = i - gxs + gnx*(j - gys);
779: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
780: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
782: for (k=0; k<nc; k++) {
783: cnt = 0;
784: for (l=lstart; l<lend+1; l++) {
785: for (p=pstart; p<pend+1; p++) {
786: if (l || p) {
787: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
788: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
789: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
790: }
791: } else {
792: if (dfill) {
793: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
794: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
795: } else {
796: for (ifill_col=0; ifill_col<nc; ifill_col++)
797: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
798: }
799: }
800: }
801: }
802: row = k + nc*(slot);
803: MatPreallocateSetLocal(ltog,1,&row,cnt,cols,dnz,onz);
804: }
805: }
806: }
807: MatSeqAIJSetPreallocation(J,0,dnz);
808: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
809: MatPreallocateFinalize(dnz,onz);
810: MatSetLocalToGlobalMapping(J,ltog);
812: /*
813: For each node in the grid: we get the neighbors in the local (on processor ordering
814: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
815: PETSc ordering.
816: */
817: if (!da->prealloc_only) {
818: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
819: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
820: for (i=xs; i<xs+nx; i++) {
821:
822: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
823: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
824:
825: for (j=ys; j<ys+ny; j++) {
826: slot = i - gxs + gnx*(j - gys);
827:
828: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
829: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
831: for (k=0; k<nc; k++) {
832: cnt = 0;
833: for (l=lstart; l<lend+1; l++) {
834: for (p=pstart; p<pend+1; p++) {
835: if (l || p) {
836: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
837: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
838: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
839: }
840: } else {
841: if (dfill) {
842: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
843: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
844: } else {
845: for (ifill_col=0; ifill_col<nc; ifill_col++)
846: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
847: }
848: }
849: }
850: }
851: row = k + nc*(slot);
852: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
853: }
854: }
855: }
856: PetscFree(values);
857: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
858: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
859: }
860: PetscFree(cols);
861: return(0);
862: }
864: /* ---------------------------------------------------------------------------------*/
868: PetscErrorCode DAGetMatrix3d_MPIAIJ(DA da,Mat J)
869: {
870: PetscErrorCode ierr;
871: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
872: PetscInt m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
873: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
874: MPI_Comm comm;
875: PetscScalar *values;
876: DAPeriodicType wrap;
877: ISLocalToGlobalMapping ltog;
878: DAStencilType st;
881: /*
882: nc - number of components per grid point
883: col - number of colors needed in one direction for single component problem
884:
885: */
886: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
887: col = 2*s + 1;
889: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
890: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
891: PetscObjectGetComm((PetscObject)da,&comm);
893: PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);
894: DAGetISLocalToGlobalMapping(da,<og);
896: /* determine the matrix preallocation information */
897: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
898: for (i=xs; i<xs+nx; i++) {
899: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
900: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
901: for (j=ys; j<ys+ny; j++) {
902: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
903: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
904: for (k=zs; k<zs+nz; k++) {
905: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
906: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
907:
908: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
909:
910: cnt = 0;
911: for (l=0; l<nc; l++) {
912: for (ii=istart; ii<iend+1; ii++) {
913: for (jj=jstart; jj<jend+1; jj++) {
914: for (kk=kstart; kk<kend+1; kk++) {
915: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
916: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
917: }
918: }
919: }
920: }
921: rows[l] = l + nc*(slot);
922: }
923: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
924: }
925: }
926: }
927: MatSeqAIJSetPreallocation(J,0,dnz);
928: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
929: MatPreallocateFinalize(dnz,onz);
930: MatSetBlockSize(J,nc);CHKERRQ(ierr)
931: MatSetLocalToGlobalMapping(J,ltog);
933: /*
934: For each node in the grid: we get the neighbors in the local (on processor ordering
935: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
936: PETSc ordering.
937: */
938: if (!da->prealloc_only) {
939: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
940: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
941: for (i=xs; i<xs+nx; i++) {
942: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
943: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
944: for (j=ys; j<ys+ny; j++) {
945: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
946: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
947: for (k=zs; k<zs+nz; k++) {
948: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
949: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
950:
951: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
952:
953: cnt = 0;
954: for (l=0; l<nc; l++) {
955: for (ii=istart; ii<iend+1; ii++) {
956: for (jj=jstart; jj<jend+1; jj++) {
957: for (kk=kstart; kk<kend+1; kk++) {
958: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
959: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
960: }
961: }
962: }
963: }
964: rows[l] = l + nc*(slot);
965: }
966: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
967: }
968: }
969: }
970: PetscFree(values);
971: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
972: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
973: }
974: PetscFree2(rows,cols);
975: return(0);
976: }
978: /* ---------------------------------------------------------------------------------*/
982: PetscErrorCode DAGetMatrix1d_MPIAIJ(DA da,Mat J)
983: {
984: PetscErrorCode ierr;
985: PetscInt xs,nx,i,i1,slot,gxs,gnx;
986: PetscInt m,dim,s,*cols,nc,*rows,col,cnt,l;
987: PetscInt istart,iend;
988: PetscScalar *values;
989: DAPeriodicType wrap;
990: ISLocalToGlobalMapping ltog;
993: /*
994: nc - number of components per grid point
995: col - number of colors needed in one direction for single component problem
996:
997: */
998: DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
999: col = 2*s + 1;
1001: DAGetCorners(da,&xs,0,0,&nx,0,0);
1002: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1004: MatSeqAIJSetPreallocation(J,col*nc,0);
1005: MatMPIAIJSetPreallocation(J,col*nc,0,0,0);
1006: MatSetBlockSize(J,nc);
1007: PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);
1008:
1009: DAGetISLocalToGlobalMapping(da,<og);
1010: MatSetLocalToGlobalMapping(J,ltog);
1011:
1012: /*
1013: For each node in the grid: we get the neighbors in the local (on processor ordering
1014: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1015: PETSc ordering.
1016: */
1017: if (!da->prealloc_only) {
1018: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
1019: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
1020: for (i=xs; i<xs+nx; i++) {
1021: istart = PetscMax(-s,gxs - i);
1022: iend = PetscMin(s,gxs + gnx - i - 1);
1023: slot = i - gxs;
1024:
1025: cnt = 0;
1026: for (l=0; l<nc; l++) {
1027: for (i1=istart; i1<iend+1; i1++) {
1028: cols[cnt++] = l + nc*(slot + i1);
1029: }
1030: rows[l] = l + nc*(slot);
1031: }
1032: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1033: }
1034: PetscFree(values);
1035: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1036: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1037: }
1038: PetscFree2(rows,cols);
1039: return(0);
1040: }
1044: PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA da,Mat J)
1045: {
1046: PetscErrorCode ierr;
1047: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1048: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1049: PetscInt istart,iend,jstart,jend,ii,jj;
1050: MPI_Comm comm;
1051: PetscScalar *values;
1052: DAPeriodicType wrap;
1053: DAStencilType st;
1054: ISLocalToGlobalMapping ltog,ltogb;
1057: /*
1058: nc - number of components per grid point
1059: col - number of colors needed in one direction for single component problem
1060: */
1061: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1062: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1063: col = 2*s + 1;
1065: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1066: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1067: PetscObjectGetComm((PetscObject)da,&comm);
1069: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1071: DAGetISLocalToGlobalMapping(da,<og);
1072: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1074: /* determine the matrix preallocation information */
1075: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1076: for (i=xs; i<xs+nx; i++) {
1077: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1078: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1079: for (j=ys; j<ys+ny; j++) {
1080: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1081: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1082: slot = i - gxs + gnx*(j - gys);
1084: /* Find block columns in block row */
1085: cnt = 0;
1086: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1087: for (ii=istart; ii<iend+1; ii++) {
1088: for (jj=jstart; jj<jend+1; jj++) {
1089: cols[cnt++] = slot + ii + gnx*jj;
1090: }
1091: }
1092: } else { /* Star stencil */
1093: cnt = 0;
1094: for (ii=istart; ii<iend+1; ii++) {
1095: if (ii) { /* jj must be zero */
1096: cols[cnt++] = slot + ii;
1097: } else {
1098: for (jj=jstart; jj<jend+1; jj++) { /* ii must be zero */
1099: cols[cnt++] = slot + gnx*jj;
1100: }
1101: }
1102: }
1103: }
1104: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1105: }
1106: }
1107: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1108: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1109: MatPreallocateFinalize(dnz,onz);
1111: MatSetLocalToGlobalMapping(J,ltog);
1112: MatSetLocalToGlobalMappingBlock(J,ltogb);
1114: /*
1115: For each node in the grid: we get the neighbors in the local (on processor ordering
1116: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1117: PETSc ordering.
1118: */
1119: if (!da->prealloc_only) {
1120: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1121: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1122: for (i=xs; i<xs+nx; i++) {
1123: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1124: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1125: for (j=ys; j<ys+ny; j++) {
1126: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1127: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1128: slot = i - gxs + gnx*(j - gys);
1129: cnt = 0;
1130: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1131: for (ii=istart; ii<iend+1; ii++) {
1132: for (jj=jstart; jj<jend+1; jj++) {
1133: cols[cnt++] = slot + ii + gnx*jj;
1134: }
1135: }
1136: } else { /* Star stencil */
1137: cnt = 0;
1138: for (ii=istart; ii<iend+1; ii++) {
1139: if (ii) { /* jj must be zero */
1140: cols[cnt++] = slot + ii;
1141: } else {
1142: for (jj=jstart; jj<jend+1; jj++) {/* ii must be zero */
1143: cols[cnt++] = slot + gnx*jj;
1144: }
1145: }
1146: }
1147: }
1148: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1149: }
1150: }
1151: PetscFree(values);
1152: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1153: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1154: }
1155: PetscFree(cols);
1156: return(0);
1157: }
1161: PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA da,Mat J)
1162: {
1163: PetscErrorCode ierr;
1164: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1165: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1166: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1167: MPI_Comm comm;
1168: PetscScalar *values;
1169: DAPeriodicType wrap;
1170: DAStencilType st;
1171: ISLocalToGlobalMapping ltog,ltogb;
1174: /*
1175: nc - number of components per grid point
1176: col - number of colors needed in one direction for single component problem
1177:
1178: */
1179: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1180: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1181: col = 2*s + 1;
1183: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1184: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1185: PetscObjectGetComm((PetscObject)da,&comm);
1187: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1189: DAGetISLocalToGlobalMapping(da,<og);
1190: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1192: /* determine the matrix preallocation information */
1193: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1194: for (i=xs; i<xs+nx; i++) {
1195: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1196: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1197: for (j=ys; j<ys+ny; j++) {
1198: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1199: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1200: for (k=zs; k<zs+nz; k++) {
1201: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1202: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1204: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1206: /* Find block columns in block row */
1207: cnt = 0;
1208: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1209: for (ii=istart; ii<iend+1; ii++) {
1210: for (jj=jstart; jj<jend+1; jj++) {
1211: for (kk=kstart; kk<kend+1; kk++) {
1212: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1213: }
1214: }
1215: }
1216: } else { /* Star stencil */
1217: cnt = 0;
1218: for (ii=istart; ii<iend+1; ii++) {
1219: if (ii) {
1220: /* jj and kk must be zero */
1221: /* cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; */
1222: cols[cnt++] = slot + ii;
1223: } else {
1224: for (jj=jstart; jj<jend+1; jj++) {
1225: if (jj) {
1226: /* ii and kk must be zero */
1227: cols[cnt++] = slot + gnx*jj;
1228: } else {
1229: /* ii and jj must be zero */
1230: for (kk=kstart; kk<kend+1; kk++) {
1231: cols[cnt++] = slot + gnx*gny*kk;
1232: }
1233: }
1234: }
1235: }
1236: }
1237: }
1238: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1239: }
1240: }
1241: }
1242: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1243: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1244: MatPreallocateFinalize(dnz,onz);
1246: MatSetLocalToGlobalMapping(J,ltog);
1247: MatSetLocalToGlobalMappingBlock(J,ltogb);
1249: /*
1250: For each node in the grid: we get the neighbors in the local (on processor ordering
1251: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1252: PETSc ordering.
1253: */
1254: if (!da->prealloc_only) {
1255: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1256: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1257: for (i=xs; i<xs+nx; i++) {
1258: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1259: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1260: for (j=ys; j<ys+ny; j++) {
1261: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1262: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1263: for (k=zs; k<zs+nz; k++) {
1264: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1265: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1266:
1267: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1268:
1269: cnt = 0;
1270: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1271: for (ii=istart; ii<iend+1; ii++) {
1272: for (jj=jstart; jj<jend+1; jj++) {
1273: for (kk=kstart; kk<kend+1; kk++) {
1274: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1275: }
1276: }
1277: }
1278: } else { /* Star stencil */
1279: cnt = 0;
1280: for (ii=istart; ii<iend+1; ii++) {
1281: if (ii) {
1282: /* jj and kk must be zero */
1283: cols[cnt++] = slot + ii;
1284: } else {
1285: for (jj=jstart; jj<jend+1; jj++) {
1286: if (jj) {
1287: /* ii and kk must be zero */
1288: cols[cnt++] = slot + gnx*jj;
1289: } else {
1290: /* ii and jj must be zero */
1291: for (kk=kstart; kk<kend+1; kk++) {
1292: cols[cnt++] = slot + gnx*gny*kk;
1293: }
1294: }
1295: }
1296: }
1297: }
1298: }
1299: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1300: }
1301: }
1302: }
1303: PetscFree(values);
1304: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1305: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1306: }
1307: PetscFree(cols);
1308: return(0);
1309: }
1312: PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA da,Mat J)
1313: {
1314: PetscErrorCode ierr;
1315: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1316: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1317: PetscInt iend,jend,ii,jj;
1318: MPI_Comm comm;
1319: PetscScalar *values;
1320: DAPeriodicType wrap;
1321: DAStencilType st;
1322: ISLocalToGlobalMapping ltog,ltogb;
1325: /*
1326: nc - number of components per grid point
1327: col - number of colors needed in one direction for single component problem
1328: */
1329: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1330: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1331: col = 2*s + 1;
1333: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1334: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1335: PetscObjectGetComm((PetscObject)da,&comm);
1337: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1339: DAGetISLocalToGlobalMapping(da,<og);
1340: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1342: /* determine the matrix preallocation information */
1343: MatPreallocateSymmetricInitialize(comm,nx*ny,nx*ny,dnz,onz);
1344: for (i=xs; i<xs+nx; i++) {
1345: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1346: for (j=ys; j<ys+ny; j++) {
1347: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1348: slot = i - gxs + gnx*(j - gys);
1350: /* Find block columns in block row */
1351: cnt = 0;
1352: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1353: for (ii=0; ii<iend+1; ii++) {
1354: for (jj=0; jj<jend+1; jj++) {
1355: cols[cnt++] = slot + ii + gnx*jj;
1356: }
1357: }
1358: } else { /* Star stencil */
1359: cnt = 0;
1360: for (ii=0; ii<iend+1; ii++) {
1361: if (ii) { /* jj must be zero */
1362: cols[cnt++] = slot + ii;
1363: } else {
1364: for (jj=0; jj<jend+1; jj++) { /* ii must be zero */
1365: cols[cnt++] = slot + gnx*jj;
1366: }
1367: }
1368: }
1369: }
1370: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1371: }
1372: }
1373: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1374: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1375: MatPreallocateFinalize(dnz,onz);
1377: MatSetLocalToGlobalMapping(J,ltog);
1378: MatSetLocalToGlobalMappingBlock(J,ltogb);
1380: /*
1381: For each node in the grid: we get the neighbors in the local (on processor ordering
1382: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1383: PETSc ordering.
1384: */
1385: if (!da->prealloc_only) {
1386: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1387: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1388: for (i=xs; i<xs+nx; i++) {
1389: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1390: for (j=ys; j<ys+ny; j++) {
1391: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1392: slot = i - gxs + gnx*(j - gys);
1393: cnt = 0;
1394: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1395: for (ii=0; ii<iend+1; ii++) {
1396: for (jj=0; jj<jend+1; jj++) {
1397: cols[cnt++] = slot + ii + gnx*jj;
1398: }
1399: }
1400: } else { /* Star stencil */
1401: cnt = 0;
1402: for (ii=0; ii<iend+1; ii++) {
1403: if (ii) { /* jj must be zero */
1404: cols[cnt++] = slot + ii;
1405: } else {
1406: for (jj=0; jj<jend+1; jj++) {/* ii must be zero */
1407: cols[cnt++] = slot + gnx*jj;
1408: }
1409: }
1410: }
1411: }
1412: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1413: }
1414: }
1415: PetscFree(values);
1416: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1417: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1418: }
1419: PetscFree(cols);
1420: return(0);
1421: }
1425: PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA da,Mat J)
1426: {
1427: PetscErrorCode ierr;
1428: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1429: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1430: PetscInt iend,jend,kend,zs,nz,gzs,gnz,ii,jj,kk;
1431: MPI_Comm comm;
1432: PetscScalar *values;
1433: DAPeriodicType wrap;
1434: DAStencilType st;
1435: ISLocalToGlobalMapping ltog,ltogb;
1438: /*
1439: nc - number of components per grid point
1440: col - number of colors needed in one direction for single component problem
1441: */
1442: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1443: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1444: col = 2*s + 1;
1446: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1447: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1448: PetscObjectGetComm((PetscObject)da,&comm);
1450: /* create the matrix */
1451: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1453: DAGetISLocalToGlobalMapping(da,<og);
1454: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1456: /* determine the matrix preallocation information */
1457: MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1458: for (i=xs; i<xs+nx; i++) {
1459: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1460: for (j=ys; j<ys+ny; j++) {
1461: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1462: for (k=zs; k<zs+nz; k++) {
1463: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1465: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1467: /* Find block columns in block row */
1468: cnt = 0;
1469: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1470: for (ii=0; ii<iend+1; ii++) {
1471: for (jj=0; jj<jend+1; jj++) {
1472: for (kk=0; kk<kend+1; kk++) {
1473: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1474: }
1475: }
1476: }
1477: } else { /* Star stencil */
1478: cnt = 0;
1479: for (ii=0; ii<iend+1; ii++) {
1480: if (ii) {
1481: /* jj and kk must be zero */
1482: cols[cnt++] = slot + ii;
1483: } else {
1484: for (jj=0; jj<jend+1; jj++) {
1485: if (jj) {
1486: /* ii and kk must be zero */
1487: cols[cnt++] = slot + gnx*jj;
1488: } else {
1489: /* ii and jj must be zero */
1490: for (kk=0; kk<kend+1; kk++) {
1491: cols[cnt++] = slot + gnx*gny*kk;
1492: }
1493: }
1494: }
1495: }
1496: }
1497: }
1498: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1499: }
1500: }
1501: }
1502: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1503: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1504: MatPreallocateFinalize(dnz,onz);
1506: MatSetLocalToGlobalMapping(J,ltog);
1507: MatSetLocalToGlobalMappingBlock(J,ltogb);
1509: /*
1510: For each node in the grid: we get the neighbors in the local (on processor ordering
1511: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1512: PETSc ordering.
1513: */
1514: if (!da->prealloc_only) {
1515: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1516: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1517: for (i=xs; i<xs+nx; i++) {
1518: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1519: for (j=ys; j<ys+ny; j++) {
1520: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1521: for (k=zs; k<zs+nz; k++) {
1522: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1523:
1524: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1525:
1526: cnt = 0;
1527: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1528: for (ii=0; ii<iend+1; ii++) {
1529: for (jj=0; jj<jend+1; jj++) {
1530: for (kk=0; kk<kend+1; kk++) {
1531: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1532: }
1533: }
1534: }
1535: } else { /* Star stencil */
1536: cnt = 0;
1537: for (ii=0; ii<iend+1; ii++) {
1538: if (ii) {
1539: /* jj and kk must be zero */
1540: cols[cnt++] = slot + ii;
1541: } else {
1542: for (jj=0; jj<jend+1; jj++) {
1543: if (jj) {
1544: /* ii and kk must be zero */
1545: cols[cnt++] = slot + gnx*jj;
1546: } else {
1547: /* ii and jj must be zero */
1548: for (kk=0; kk<kend+1; kk++) {
1549: cols[cnt++] = slot + gnx*gny*kk;
1550: }
1551: }
1552: }
1553: }
1554: }
1555: }
1556: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1557: }
1558: }
1559: }
1560: PetscFree(values);
1561: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1562: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1563: }
1564: PetscFree(cols);
1565: return(0);
1566: }
1568: /* ---------------------------------------------------------------------------------*/
1572: PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA da,Mat J)
1573: {
1574: PetscErrorCode ierr;
1575: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1576: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz;
1577: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1578: PetscInt ifill_col,*dfill = da->dfill,*ofill = da->ofill;
1579: MPI_Comm comm;
1580: PetscScalar *values;
1581: DAPeriodicType wrap;
1582: ISLocalToGlobalMapping ltog;
1583: DAStencilType st;
1586: /*
1587: nc - number of components per grid point
1588: col - number of colors needed in one direction for single component problem
1589:
1590: */
1591: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1592: col = 2*s + 1;
1593: if (DAXPeriodic(wrap) && (m % col)){
1594: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1595: by 2*stencil_width + 1\n");
1596: }
1597: if (DAYPeriodic(wrap) && (n % col)){
1598: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1599: by 2*stencil_width + 1\n");
1600: }
1601: if (DAZPeriodic(wrap) && (p % col)){
1602: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1603: by 2*stencil_width + 1\n");
1604: }
1606: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1607: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1608: PetscObjectGetComm((PetscObject)da,&comm);
1610: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
1611: DAGetISLocalToGlobalMapping(da,<og);
1613: /* determine the matrix preallocation information */
1614: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1617: for (i=xs; i<xs+nx; i++) {
1618: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1619: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1620: for (j=ys; j<ys+ny; j++) {
1621: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1622: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1623: for (k=zs; k<zs+nz; k++) {
1624: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1625: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1626:
1627: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1628:
1629: for (l=0; l<nc; l++) {
1630: cnt = 0;
1631: for (ii=istart; ii<iend+1; ii++) {
1632: for (jj=jstart; jj<jend+1; jj++) {
1633: for (kk=kstart; kk<kend+1; kk++) {
1634: if (ii || jj || kk) {
1635: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1636: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1637: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1638: }
1639: } else {
1640: if (dfill) {
1641: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1642: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1643: } else {
1644: for (ifill_col=0; ifill_col<nc; ifill_col++)
1645: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1646: }
1647: }
1648: }
1649: }
1650: }
1651: row = l + nc*(slot);
1652: MatPreallocateSetLocal(ltog,1,&row,cnt,cols,dnz,onz);
1653: }
1654: }
1655: }
1656: }
1657: MatSeqAIJSetPreallocation(J,0,dnz);
1658: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1659: MatPreallocateFinalize(dnz,onz);
1660: MatSetLocalToGlobalMapping(J,ltog);
1662: /*
1663: For each node in the grid: we get the neighbors in the local (on processor ordering
1664: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1665: PETSc ordering.
1666: */
1667: if (!da->prealloc_only) {
1668: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1669: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1670: for (i=xs; i<xs+nx; i++) {
1671: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1672: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1673: for (j=ys; j<ys+ny; j++) {
1674: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1675: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1676: for (k=zs; k<zs+nz; k++) {
1677: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1678: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1679:
1680: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1681:
1682: for (l=0; l<nc; l++) {
1683: cnt = 0;
1684: for (ii=istart; ii<iend+1; ii++) {
1685: for (jj=jstart; jj<jend+1; jj++) {
1686: for (kk=kstart; kk<kend+1; kk++) {
1687: if (ii || jj || kk) {
1688: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1689: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1690: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1691: }
1692: } else {
1693: if (dfill) {
1694: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1695: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1696: } else {
1697: for (ifill_col=0; ifill_col<nc; ifill_col++)
1698: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1699: }
1700: }
1701: }
1702: }
1703: }
1704: row = l + nc*(slot);
1705: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1706: }
1707: }
1708: }
1709: }
1710: PetscFree(values);
1711: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1712: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1713: }
1714: PetscFree(cols);
1715: return(0);
1716: }