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, const MatType mtype,Mat *J)
545: {
547: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3];
548: Mat A;
549: MPI_Comm comm;
550: const 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,ltogb;
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: DAGetISLocalToGlobalMappingBlck(da,<ogb);
682:
683: /* determine the matrix preallocation information */
684: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
685: for (i=xs; i<xs+nx; i++) {
687: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
688: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
690: for (j=ys; j<ys+ny; j++) {
691: slot = i - gxs + gnx*(j - gys);
693: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
694: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
696: cnt = 0;
697: for (k=0; k<nc; k++) {
698: for (l=lstart; l<lend+1; l++) {
699: for (p=pstart; p<pend+1; p++) {
700: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
701: cols[cnt++] = k + nc*(slot + gnx*l + p);
702: }
703: }
704: }
705: rows[k] = k + nc*(slot);
706: }
707: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
708: }
709: }
710: MatSeqAIJSetPreallocation(J,0,dnz);
711: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
712: MatSetBlockSize(J,nc);
713: MatPreallocateFinalize(dnz,onz);
715: MatSetLocalToGlobalMapping(J,ltog);
716: MatSetLocalToGlobalMappingBlock(J,ltogb);
718: /*
719: For each node in the grid: we get the neighbors in the local (on processor ordering
720: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
721: PETSc ordering.
722: */
723: if (!da->prealloc_only) {
724: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
725: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
726: for (i=xs; i<xs+nx; i++) {
727:
728: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
729: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
730:
731: for (j=ys; j<ys+ny; j++) {
732: slot = i - gxs + gnx*(j - gys);
733:
734: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
735: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
737: cnt = 0;
738: for (k=0; k<nc; k++) {
739: for (l=lstart; l<lend+1; l++) {
740: for (p=pstart; p<pend+1; p++) {
741: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
742: cols[cnt++] = k + nc*(slot + gnx*l + p);
743: }
744: }
745: }
746: rows[k] = k + nc*(slot);
747: }
748: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
749: }
750: }
751: PetscFree(values);
752: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
753: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
754: }
755: PetscFree2(rows,cols);
756: return(0);
757: }
761: PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA da,Mat J)
762: {
763: PetscErrorCode ierr;
764: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
765: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p;
766: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
767: PetscInt ifill_col,*ofill = da->ofill, *dfill = da->dfill;
768: MPI_Comm comm;
769: PetscScalar *values;
770: DAPeriodicType wrap;
771: ISLocalToGlobalMapping ltog,ltogb;
772: DAStencilType st;
775: /*
776: nc - number of components per grid point
777: col - number of colors needed in one direction for single component problem
778:
779: */
780: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
781: col = 2*s + 1;
782: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
783: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
784: PetscObjectGetComm((PetscObject)da,&comm);
786: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
787: DAGetISLocalToGlobalMapping(da,<og);
788: DAGetISLocalToGlobalMappingBlck(da,<ogb);
789:
790: /* determine the matrix preallocation information */
791: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
792: for (i=xs; i<xs+nx; i++) {
794: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
795: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
797: for (j=ys; j<ys+ny; j++) {
798: slot = i - gxs + gnx*(j - gys);
800: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
801: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
803: for (k=0; k<nc; k++) {
804: cnt = 0;
805: for (l=lstart; l<lend+1; l++) {
806: for (p=pstart; p<pend+1; p++) {
807: if (l || p) {
808: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
809: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
810: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
811: }
812: } else {
813: if (dfill) {
814: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
815: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
816: } else {
817: for (ifill_col=0; ifill_col<nc; ifill_col++)
818: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
819: }
820: }
821: }
822: }
823: row = k + nc*(slot);
824: MatPreallocateSetLocal(ltog,1,&row,cnt,cols,dnz,onz);
825: }
826: }
827: }
828: MatSeqAIJSetPreallocation(J,0,dnz);
829: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
830: MatPreallocateFinalize(dnz,onz);
831: MatSetLocalToGlobalMapping(J,ltog);
832: MatSetLocalToGlobalMappingBlock(J,ltogb);
834: /*
835: For each node in the grid: we get the neighbors in the local (on processor ordering
836: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
837: PETSc ordering.
838: */
839: if (!da->prealloc_only) {
840: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
841: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
842: for (i=xs; i<xs+nx; i++) {
843:
844: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
845: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
846:
847: for (j=ys; j<ys+ny; j++) {
848: slot = i - gxs + gnx*(j - gys);
849:
850: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
851: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
853: for (k=0; k<nc; k++) {
854: cnt = 0;
855: for (l=lstart; l<lend+1; l++) {
856: for (p=pstart; p<pend+1; p++) {
857: if (l || p) {
858: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
859: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
860: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
861: }
862: } else {
863: if (dfill) {
864: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
865: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
866: } else {
867: for (ifill_col=0; ifill_col<nc; ifill_col++)
868: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
869: }
870: }
871: }
872: }
873: row = k + nc*(slot);
874: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
875: }
876: }
877: }
878: PetscFree(values);
879: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
880: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
881: }
882: PetscFree(cols);
883: return(0);
884: }
886: /* ---------------------------------------------------------------------------------*/
890: PetscErrorCode DAGetMatrix3d_MPIAIJ(DA da,Mat J)
891: {
892: PetscErrorCode ierr;
893: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
894: PetscInt m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p,*dnz = PETSC_NULL,*onz = PETSC_NULL;
895: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
896: MPI_Comm comm;
897: PetscScalar *values;
898: DAPeriodicType wrap;
899: ISLocalToGlobalMapping ltog,ltogb;
900: DAStencilType st;
903: /*
904: nc - number of components per grid point
905: col - number of colors needed in one direction for single component problem
906:
907: */
908: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
909: col = 2*s + 1;
911: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
912: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
913: PetscObjectGetComm((PetscObject)da,&comm);
915: PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);
916: DAGetISLocalToGlobalMapping(da,<og);
917: DAGetISLocalToGlobalMappingBlck(da,<ogb);
919: /* determine the matrix preallocation information */
920: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
921: for (i=xs; i<xs+nx; i++) {
922: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
923: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
924: for (j=ys; j<ys+ny; j++) {
925: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
926: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
927: for (k=zs; k<zs+nz; k++) {
928: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
929: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
930:
931: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
932:
933: cnt = 0;
934: for (l=0; l<nc; l++) {
935: for (ii=istart; ii<iend+1; ii++) {
936: for (jj=jstart; jj<jend+1; jj++) {
937: for (kk=kstart; kk<kend+1; kk++) {
938: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
939: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
940: }
941: }
942: }
943: }
944: rows[l] = l + nc*(slot);
945: }
946: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
947: }
948: }
949: }
950: MatSeqAIJSetPreallocation(J,0,dnz);
951: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
952: MatPreallocateFinalize(dnz,onz);
953: MatSetBlockSize(J,nc);CHKERRQ(ierr)
954: MatSetLocalToGlobalMapping(J,ltog);
955: MatSetLocalToGlobalMappingBlock(J,ltogb);
957: /*
958: For each node in the grid: we get the neighbors in the local (on processor ordering
959: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
960: PETSc ordering.
961: */
962: if (!da->prealloc_only) {
963: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
964: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
965: for (i=xs; i<xs+nx; i++) {
966: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
967: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
968: for (j=ys; j<ys+ny; j++) {
969: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
970: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
971: for (k=zs; k<zs+nz; k++) {
972: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
973: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
974:
975: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
976:
977: cnt = 0;
978: for (l=0; l<nc; l++) {
979: for (ii=istart; ii<iend+1; ii++) {
980: for (jj=jstart; jj<jend+1; jj++) {
981: for (kk=kstart; kk<kend+1; kk++) {
982: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
983: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
984: }
985: }
986: }
987: }
988: rows[l] = l + nc*(slot);
989: }
990: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
991: }
992: }
993: }
994: PetscFree(values);
995: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
996: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
997: }
998: PetscFree2(rows,cols);
999: return(0);
1000: }
1002: /* ---------------------------------------------------------------------------------*/
1006: PetscErrorCode DAGetMatrix1d_MPIAIJ(DA da,Mat J)
1007: {
1008: PetscErrorCode ierr;
1009: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1010: PetscInt m,dim,s,*cols = PETSC_NULL,nc,*rows = PETSC_NULL,col,cnt,l;
1011: PetscInt istart,iend;
1012: PetscScalar *values;
1013: DAPeriodicType wrap;
1014: ISLocalToGlobalMapping ltog,ltogb;
1017: /*
1018: nc - number of components per grid point
1019: col - number of colors needed in one direction for single component problem
1020:
1021: */
1022: DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
1023: col = 2*s + 1;
1025: DAGetCorners(da,&xs,0,0,&nx,0,0);
1026: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1028: MatSeqAIJSetPreallocation(J,col*nc,0);
1029: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1030: MatSetBlockSize(J,nc);
1031: PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);
1032:
1033: DAGetISLocalToGlobalMapping(da,<og);
1034: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1035: MatSetLocalToGlobalMapping(J,ltog);
1036: MatSetLocalToGlobalMappingBlock(J,ltogb);
1037:
1038: /*
1039: For each node in the grid: we get the neighbors in the local (on processor ordering
1040: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1041: PETSc ordering.
1042: */
1043: if (!da->prealloc_only) {
1044: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
1045: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
1046: for (i=xs; i<xs+nx; i++) {
1047: istart = PetscMax(-s,gxs - i);
1048: iend = PetscMin(s,gxs + gnx - i - 1);
1049: slot = i - gxs;
1050:
1051: cnt = 0;
1052: for (l=0; l<nc; l++) {
1053: for (i1=istart; i1<iend+1; i1++) {
1054: cols[cnt++] = l + nc*(slot + i1);
1055: }
1056: rows[l] = l + nc*(slot);
1057: }
1058: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1059: }
1060: PetscFree(values);
1061: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1062: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1063: }
1064: PetscFree2(rows,cols);
1065: return(0);
1066: }
1070: PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA da,Mat J)
1071: {
1072: PetscErrorCode ierr;
1073: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1074: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1075: PetscInt istart,iend,jstart,jend,ii,jj;
1076: MPI_Comm comm;
1077: PetscScalar *values;
1078: DAPeriodicType wrap;
1079: DAStencilType st;
1080: ISLocalToGlobalMapping ltog,ltogb;
1083: /*
1084: nc - number of components per grid point
1085: col - number of colors needed in one direction for single component problem
1086: */
1087: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1088: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1089: col = 2*s + 1;
1091: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1092: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1093: PetscObjectGetComm((PetscObject)da,&comm);
1095: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1097: DAGetISLocalToGlobalMapping(da,<og);
1098: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1100: /* determine the matrix preallocation information */
1101: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1102: for (i=xs; i<xs+nx; i++) {
1103: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1104: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1105: for (j=ys; j<ys+ny; j++) {
1106: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1107: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1108: slot = i - gxs + gnx*(j - gys);
1110: /* Find block columns in block row */
1111: cnt = 0;
1112: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1113: for (ii=istart; ii<iend+1; ii++) {
1114: for (jj=jstart; jj<jend+1; jj++) {
1115: cols[cnt++] = slot + ii + gnx*jj;
1116: }
1117: }
1118: } else { /* Star stencil */
1119: cnt = 0;
1120: for (ii=istart; ii<iend+1; ii++) {
1121: if (ii) { /* jj must be zero */
1122: cols[cnt++] = slot + ii;
1123: } else {
1124: for (jj=jstart; jj<jend+1; jj++) { /* ii must be zero */
1125: cols[cnt++] = slot + gnx*jj;
1126: }
1127: }
1128: }
1129: }
1130: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1131: }
1132: }
1133: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1134: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1135: MatPreallocateFinalize(dnz,onz);
1137: MatSetLocalToGlobalMapping(J,ltog);
1138: MatSetLocalToGlobalMappingBlock(J,ltogb);
1140: /*
1141: For each node in the grid: we get the neighbors in the local (on processor ordering
1142: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1143: PETSc ordering.
1144: */
1145: if (!da->prealloc_only) {
1146: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1147: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1148: for (i=xs; i<xs+nx; i++) {
1149: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1150: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1151: for (j=ys; j<ys+ny; j++) {
1152: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1153: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1154: slot = i - gxs + gnx*(j - gys);
1155: cnt = 0;
1156: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1157: for (ii=istart; ii<iend+1; ii++) {
1158: for (jj=jstart; jj<jend+1; jj++) {
1159: cols[cnt++] = slot + ii + gnx*jj;
1160: }
1161: }
1162: } else { /* Star stencil */
1163: cnt = 0;
1164: for (ii=istart; ii<iend+1; ii++) {
1165: if (ii) { /* jj must be zero */
1166: cols[cnt++] = slot + ii;
1167: } else {
1168: for (jj=jstart; jj<jend+1; jj++) {/* ii must be zero */
1169: cols[cnt++] = slot + gnx*jj;
1170: }
1171: }
1172: }
1173: }
1174: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1175: }
1176: }
1177: PetscFree(values);
1178: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1179: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1180: }
1181: PetscFree(cols);
1182: return(0);
1183: }
1187: PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA da,Mat J)
1188: {
1189: PetscErrorCode ierr;
1190: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1191: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1192: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1193: MPI_Comm comm;
1194: PetscScalar *values;
1195: DAPeriodicType wrap;
1196: DAStencilType st;
1197: ISLocalToGlobalMapping ltog,ltogb;
1200: /*
1201: nc - number of components per grid point
1202: col - number of colors needed in one direction for single component problem
1203:
1204: */
1205: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1206: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1207: col = 2*s + 1;
1209: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1210: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1211: PetscObjectGetComm((PetscObject)da,&comm);
1213: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1215: DAGetISLocalToGlobalMapping(da,<og);
1216: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1218: /* determine the matrix preallocation information */
1219: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1220: for (i=xs; i<xs+nx; i++) {
1221: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1222: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1223: for (j=ys; j<ys+ny; j++) {
1224: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1225: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1226: for (k=zs; k<zs+nz; k++) {
1227: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1228: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1230: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1232: /* Find block columns in block row */
1233: cnt = 0;
1234: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1235: for (ii=istart; ii<iend+1; ii++) {
1236: for (jj=jstart; jj<jend+1; jj++) {
1237: for (kk=kstart; kk<kend+1; kk++) {
1238: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1239: }
1240: }
1241: }
1242: } else { /* Star stencil */
1243: cnt = 0;
1244: for (ii=istart; ii<iend+1; ii++) {
1245: if (ii) {
1246: /* jj and kk must be zero */
1247: /* cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; */
1248: cols[cnt++] = slot + ii;
1249: } else {
1250: for (jj=jstart; jj<jend+1; jj++) {
1251: if (jj) {
1252: /* ii and kk must be zero */
1253: cols[cnt++] = slot + gnx*jj;
1254: } else {
1255: /* ii and jj must be zero */
1256: for (kk=kstart; kk<kend+1; kk++) {
1257: cols[cnt++] = slot + gnx*gny*kk;
1258: }
1259: }
1260: }
1261: }
1262: }
1263: }
1264: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1265: }
1266: }
1267: }
1268: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1269: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1270: MatPreallocateFinalize(dnz,onz);
1272: MatSetLocalToGlobalMapping(J,ltog);
1273: MatSetLocalToGlobalMappingBlock(J,ltogb);
1275: /*
1276: For each node in the grid: we get the neighbors in the local (on processor ordering
1277: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1278: PETSc ordering.
1279: */
1280: if (!da->prealloc_only) {
1281: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1282: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1283: for (i=xs; i<xs+nx; i++) {
1284: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1285: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1286: for (j=ys; j<ys+ny; j++) {
1287: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1288: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1289: for (k=zs; k<zs+nz; k++) {
1290: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1291: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1292:
1293: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1294:
1295: cnt = 0;
1296: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1297: for (ii=istart; ii<iend+1; ii++) {
1298: for (jj=jstart; jj<jend+1; jj++) {
1299: for (kk=kstart; kk<kend+1; kk++) {
1300: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1301: }
1302: }
1303: }
1304: } else { /* Star stencil */
1305: cnt = 0;
1306: for (ii=istart; ii<iend+1; ii++) {
1307: if (ii) {
1308: /* jj and kk must be zero */
1309: cols[cnt++] = slot + ii;
1310: } else {
1311: for (jj=jstart; jj<jend+1; jj++) {
1312: if (jj) {
1313: /* ii and kk must be zero */
1314: cols[cnt++] = slot + gnx*jj;
1315: } else {
1316: /* ii and jj must be zero */
1317: for (kk=kstart; kk<kend+1; kk++) {
1318: cols[cnt++] = slot + gnx*gny*kk;
1319: }
1320: }
1321: }
1322: }
1323: }
1324: }
1325: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1326: }
1327: }
1328: }
1329: PetscFree(values);
1330: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1331: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1332: }
1333: PetscFree(cols);
1334: return(0);
1335: }
1338: PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA da,Mat J)
1339: {
1340: PetscErrorCode ierr;
1341: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1342: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1343: PetscInt iend,jend,ii,jj;
1344: MPI_Comm comm;
1345: PetscScalar *values;
1346: DAPeriodicType wrap;
1347: DAStencilType st;
1348: ISLocalToGlobalMapping ltog,ltogb;
1351: /*
1352: nc - number of components per grid point
1353: col - number of colors needed in one direction for single component problem
1354: */
1355: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1356: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1357: col = 2*s + 1;
1359: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1360: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1361: PetscObjectGetComm((PetscObject)da,&comm);
1363: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1365: DAGetISLocalToGlobalMapping(da,<og);
1366: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1368: /* determine the matrix preallocation information */
1369: MatPreallocateSymmetricInitialize(comm,nx*ny,nx*ny,dnz,onz);
1370: for (i=xs; i<xs+nx; i++) {
1371: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1372: for (j=ys; j<ys+ny; j++) {
1373: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1374: slot = i - gxs + gnx*(j - gys);
1376: /* Find block columns in block row */
1377: cnt = 0;
1378: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1379: for (ii=0; ii<iend+1; ii++) {
1380: for (jj=0; jj<jend+1; jj++) {
1381: cols[cnt++] = slot + ii + gnx*jj;
1382: }
1383: }
1384: } else { /* Star stencil */
1385: cnt = 0;
1386: for (ii=0; ii<iend+1; ii++) {
1387: if (ii) { /* jj must be zero */
1388: cols[cnt++] = slot + ii;
1389: } else {
1390: for (jj=0; jj<jend+1; jj++) { /* ii must be zero */
1391: cols[cnt++] = slot + gnx*jj;
1392: }
1393: }
1394: }
1395: }
1396: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1397: }
1398: }
1399: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1400: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1401: MatPreallocateFinalize(dnz,onz);
1403: MatSetLocalToGlobalMapping(J,ltog);
1404: MatSetLocalToGlobalMappingBlock(J,ltogb);
1406: /*
1407: For each node in the grid: we get the neighbors in the local (on processor ordering
1408: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1409: PETSc ordering.
1410: */
1411: if (!da->prealloc_only) {
1412: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1413: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1414: for (i=xs; i<xs+nx; i++) {
1415: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1416: for (j=ys; j<ys+ny; j++) {
1417: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1418: slot = i - gxs + gnx*(j - gys);
1419: cnt = 0;
1420: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1421: for (ii=0; ii<iend+1; ii++) {
1422: for (jj=0; jj<jend+1; jj++) {
1423: cols[cnt++] = slot + ii + gnx*jj;
1424: }
1425: }
1426: } else { /* Star stencil */
1427: cnt = 0;
1428: for (ii=0; ii<iend+1; ii++) {
1429: if (ii) { /* jj must be zero */
1430: cols[cnt++] = slot + ii;
1431: } else {
1432: for (jj=0; jj<jend+1; jj++) {/* ii must be zero */
1433: cols[cnt++] = slot + gnx*jj;
1434: }
1435: }
1436: }
1437: }
1438: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1439: }
1440: }
1441: PetscFree(values);
1442: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1443: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1444: }
1445: PetscFree(cols);
1446: return(0);
1447: }
1451: PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA da,Mat J)
1452: {
1453: PetscErrorCode ierr;
1454: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1455: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1456: PetscInt iend,jend,kend,zs,nz,gzs,gnz,ii,jj,kk;
1457: MPI_Comm comm;
1458: PetscScalar *values;
1459: DAPeriodicType wrap;
1460: DAStencilType st;
1461: ISLocalToGlobalMapping ltog,ltogb;
1464: /*
1465: nc - number of components per grid point
1466: col - number of colors needed in one direction for single component problem
1467: */
1468: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1469: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1470: col = 2*s + 1;
1472: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1473: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1474: PetscObjectGetComm((PetscObject)da,&comm);
1476: /* create the matrix */
1477: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1479: DAGetISLocalToGlobalMapping(da,<og);
1480: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1482: /* determine the matrix preallocation information */
1483: MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1484: for (i=xs; i<xs+nx; i++) {
1485: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1486: for (j=ys; j<ys+ny; j++) {
1487: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1488: for (k=zs; k<zs+nz; k++) {
1489: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1491: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1493: /* Find block columns in block row */
1494: cnt = 0;
1495: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1496: for (ii=0; ii<iend+1; ii++) {
1497: for (jj=0; jj<jend+1; jj++) {
1498: for (kk=0; kk<kend+1; kk++) {
1499: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1500: }
1501: }
1502: }
1503: } else { /* Star stencil */
1504: cnt = 0;
1505: for (ii=0; ii<iend+1; ii++) {
1506: if (ii) {
1507: /* jj and kk must be zero */
1508: cols[cnt++] = slot + ii;
1509: } else {
1510: for (jj=0; jj<jend+1; jj++) {
1511: if (jj) {
1512: /* ii and kk must be zero */
1513: cols[cnt++] = slot + gnx*jj;
1514: } else {
1515: /* ii and jj must be zero */
1516: for (kk=0; kk<kend+1; kk++) {
1517: cols[cnt++] = slot + gnx*gny*kk;
1518: }
1519: }
1520: }
1521: }
1522: }
1523: }
1524: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1525: }
1526: }
1527: }
1528: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1529: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1530: MatPreallocateFinalize(dnz,onz);
1532: MatSetLocalToGlobalMapping(J,ltog);
1533: MatSetLocalToGlobalMappingBlock(J,ltogb);
1535: /*
1536: For each node in the grid: we get the neighbors in the local (on processor ordering
1537: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1538: PETSc ordering.
1539: */
1540: if (!da->prealloc_only) {
1541: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1542: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1543: for (i=xs; i<xs+nx; i++) {
1544: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1545: for (j=ys; j<ys+ny; j++) {
1546: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1547: for (k=zs; k<zs+nz; k++) {
1548: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1549:
1550: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1551:
1552: cnt = 0;
1553: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1554: for (ii=0; ii<iend+1; ii++) {
1555: for (jj=0; jj<jend+1; jj++) {
1556: for (kk=0; kk<kend+1; kk++) {
1557: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1558: }
1559: }
1560: }
1561: } else { /* Star stencil */
1562: cnt = 0;
1563: for (ii=0; ii<iend+1; ii++) {
1564: if (ii) {
1565: /* jj and kk must be zero */
1566: cols[cnt++] = slot + ii;
1567: } else {
1568: for (jj=0; jj<jend+1; jj++) {
1569: if (jj) {
1570: /* ii and kk must be zero */
1571: cols[cnt++] = slot + gnx*jj;
1572: } else {
1573: /* ii and jj must be zero */
1574: for (kk=0; kk<kend+1; kk++) {
1575: cols[cnt++] = slot + gnx*gny*kk;
1576: }
1577: }
1578: }
1579: }
1580: }
1581: }
1582: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1583: }
1584: }
1585: }
1586: PetscFree(values);
1587: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1588: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1589: }
1590: PetscFree(cols);
1591: return(0);
1592: }
1594: /* ---------------------------------------------------------------------------------*/
1598: PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA da,Mat J)
1599: {
1600: PetscErrorCode ierr;
1601: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1602: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz;
1603: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1604: PetscInt ifill_col,*dfill = da->dfill,*ofill = da->ofill;
1605: MPI_Comm comm;
1606: PetscScalar *values;
1607: DAPeriodicType wrap;
1608: ISLocalToGlobalMapping ltog,ltogb;
1609: DAStencilType st;
1612: /*
1613: nc - number of components per grid point
1614: col - number of colors needed in one direction for single component problem
1615:
1616: */
1617: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1618: col = 2*s + 1;
1619: if (DAXPeriodic(wrap) && (m % col)){
1620: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1621: by 2*stencil_width + 1\n");
1622: }
1623: if (DAYPeriodic(wrap) && (n % col)){
1624: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1625: by 2*stencil_width + 1\n");
1626: }
1627: if (DAZPeriodic(wrap) && (p % col)){
1628: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1629: by 2*stencil_width + 1\n");
1630: }
1632: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1633: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1634: PetscObjectGetComm((PetscObject)da,&comm);
1636: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
1637: DAGetISLocalToGlobalMapping(da,<og);
1638: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1640: /* determine the matrix preallocation information */
1641: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1644: for (i=xs; i<xs+nx; i++) {
1645: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1646: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1647: for (j=ys; j<ys+ny; j++) {
1648: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1649: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1650: for (k=zs; k<zs+nz; k++) {
1651: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1652: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1653:
1654: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1655:
1656: for (l=0; l<nc; l++) {
1657: cnt = 0;
1658: for (ii=istart; ii<iend+1; ii++) {
1659: for (jj=jstart; jj<jend+1; jj++) {
1660: for (kk=kstart; kk<kend+1; kk++) {
1661: if (ii || jj || kk) {
1662: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1663: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1664: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1665: }
1666: } else {
1667: if (dfill) {
1668: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1669: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1670: } else {
1671: for (ifill_col=0; ifill_col<nc; ifill_col++)
1672: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1673: }
1674: }
1675: }
1676: }
1677: }
1678: row = l + nc*(slot);
1679: MatPreallocateSetLocal(ltog,1,&row,cnt,cols,dnz,onz);
1680: }
1681: }
1682: }
1683: }
1684: MatSeqAIJSetPreallocation(J,0,dnz);
1685: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1686: MatPreallocateFinalize(dnz,onz);
1687: MatSetLocalToGlobalMapping(J,ltog);
1688: MatSetLocalToGlobalMappingBlock(J,ltogb);
1690: /*
1691: For each node in the grid: we get the neighbors in the local (on processor ordering
1692: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1693: PETSc ordering.
1694: */
1695: if (!da->prealloc_only) {
1696: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1697: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1698: for (i=xs; i<xs+nx; i++) {
1699: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1700: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1701: for (j=ys; j<ys+ny; j++) {
1702: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1703: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1704: for (k=zs; k<zs+nz; k++) {
1705: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1706: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1707:
1708: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1709:
1710: for (l=0; l<nc; l++) {
1711: cnt = 0;
1712: for (ii=istart; ii<iend+1; ii++) {
1713: for (jj=jstart; jj<jend+1; jj++) {
1714: for (kk=kstart; kk<kend+1; kk++) {
1715: if (ii || jj || kk) {
1716: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1717: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1718: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1719: }
1720: } else {
1721: if (dfill) {
1722: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1723: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1724: } else {
1725: for (ifill_col=0; ifill_col<nc; ifill_col++)
1726: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1727: }
1728: }
1729: }
1730: }
1731: }
1732: row = l + nc*(slot);
1733: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1734: }
1735: }
1736: }
1737: }
1738: PetscFree(values);
1739: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1740: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1741: }
1742: PetscFree(cols);
1743: return(0);
1744: }