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