Actual source code: fdda.c
petsc-3.11.4 2019-09-28
2: #include <petsc/private/dmdaimpl.h>
3: #include <petscmat.h>
5: extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring*);
6: extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring*);
7: extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring*);
8: extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring*);
10: /*
11: For ghost i that may be negative or greater than the upper bound this
12: maps it into the 0:m-1 range using periodicity
13: */
14: #define SetInRange(i,m) ((i < 0) ? m+i : ((i >= m) ? i-m : i))
16: static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill)
17: {
19: PetscInt i,j,nz,*fill;
22: if (!dfill) return(0);
24: /* count number nonzeros */
25: nz = 0;
26: for (i=0; i<w; i++) {
27: for (j=0; j<w; j++) {
28: if (dfill[w*i+j]) nz++;
29: }
30: }
31: PetscMalloc1(nz + w + 1,&fill);
32: /* construct modified CSR storage of nonzero structure */
33: /* fill[0 -- w] marks starts of each row of column indices (and end of last row)
34: so fill[1] - fill[0] gives number of nonzeros in first row etc */
35: nz = w + 1;
36: for (i=0; i<w; i++) {
37: fill[i] = nz;
38: for (j=0; j<w; j++) {
39: if (dfill[w*i+j]) {
40: fill[nz] = j;
41: nz++;
42: }
43: }
44: }
45: fill[w] = nz;
47: *rfill = fill;
48: return(0);
49: }
52: static PetscErrorCode DMDASetBlockFillsSparse_Private(const PetscInt *dfillsparse,PetscInt w,PetscInt **rfill)
53: {
55: PetscInt nz;
58: if (!dfillsparse) return(0);
60: /* Determine number of non-zeros */
61: nz = (dfillsparse[w] - w - 1);
63: /* Allocate space for our copy of the given sparse matrix representation. */
64: PetscMalloc1(nz + w + 1,rfill);
65: PetscMemcpy(*rfill,dfillsparse,(nz+w+1)*sizeof(PetscInt));
66: return(0);
67: }
70: static PetscErrorCode DMDASetBlockFills_Private2(DM_DA *dd)
71: {
73: PetscInt i,k,cnt = 1;
77: /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of
78: columns to the left with any nonzeros in them plus 1 */
79: PetscCalloc1(dd->w,&dd->ofillcols);
80: for (i=0; i<dd->w; i++) {
81: for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
82: }
83: for (i=0; i<dd->w; i++) {
84: if (dd->ofillcols[i]) {
85: dd->ofillcols[i] = cnt++;
86: }
87: }
88: return(0);
89: }
93: /*@
94: DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
95: of the matrix returned by DMCreateMatrix().
97: Logically Collective on DMDA
99: Input Parameter:
100: + da - the distributed array
101: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
102: - ofill - the fill pattern in the off-diagonal blocks
105: Level: developer
107: Notes:
108: This only makes sense when you are doing multicomponent problems but using the
109: MPIAIJ matrix format
111: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
112: representing coupling and 0 entries for missing coupling. For example
113: $ dfill[9] = {1, 0, 0,
114: $ 1, 1, 0,
115: $ 0, 1, 1}
116: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
117: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
118: diagonal block).
120: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
121: can be represented in the dfill, ofill format
123: Contributed by Glenn Hammond
125: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
127: @*/
128: PetscErrorCode DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
129: {
130: DM_DA *dd = (DM_DA*)da->data;
134: /* save the given dfill and ofill information */
135: DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
136: DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);
138: /* count nonzeros in ofill columns */
139: DMDASetBlockFills_Private2(dd);
141: return(0);
142: }
145: /*@
146: DMDASetBlockFillsSparse - Sets the fill pattern in each block for a multi-component problem
147: of the matrix returned by DMCreateMatrix(), using sparse representations
148: of fill patterns.
150: Logically Collective on DMDA
152: Input Parameter:
153: + da - the distributed array
154: . dfill - the sparse fill pattern in the diagonal block (may be NULL, means use dense block)
155: - ofill - the sparse fill pattern in the off-diagonal blocks
158: Level: developer
160: Notes: This only makes sense when you are doing multicomponent problems but using the
161: MPIAIJ matrix format
163: The format for dfill and ofill is a sparse representation of a
164: dof-by-dof matrix with 1 entries representing coupling and 0 entries
165: for missing coupling. The sparse representation is a 1 dimensional
166: array of length nz + dof + 1, where nz is the number of non-zeros in
167: the matrix. The first dof entries in the array give the
168: starting array indices of each row's items in the rest of the array,
169: the dof+1st item contains the value nz + dof + 1 (i.e. the entire length of the array)
170: and the remaining nz items give the column indices of each of
171: the 1s within the logical 2D matrix. Each row's items within
172: the array are the column indices of the 1s within that row
173: of the 2D matrix. PETSc developers may recognize that this is the
174: same format as that computed by the DMDASetBlockFills_Private()
175: function from a dense 2D matrix representation.
177: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
178: can be represented in the dfill, ofill format
180: Contributed by Philip C. Roth
182: .seealso DMDASetBlockFills(), DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
184: @*/
185: PetscErrorCode DMDASetBlockFillsSparse(DM da,const PetscInt *dfillsparse,const PetscInt *ofillsparse)
186: {
187: DM_DA *dd = (DM_DA*)da->data;
191: /* save the given dfill and ofill information */
192: DMDASetBlockFillsSparse_Private(dfillsparse,dd->w,&dd->dfill);
193: DMDASetBlockFillsSparse_Private(ofillsparse,dd->w,&dd->ofill);
195: /* count nonzeros in ofill columns */
196: DMDASetBlockFills_Private2(dd);
198: return(0);
199: }
202: PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
203: {
204: PetscErrorCode ierr;
205: PetscInt dim,m,n,p,nc;
206: DMBoundaryType bx,by,bz;
207: MPI_Comm comm;
208: PetscMPIInt size;
209: PetscBool isBAIJ;
210: DM_DA *dd = (DM_DA*)da->data;
213: /*
214: m
215: ------------------------------------------------------
216: | |
217: | |
218: | ---------------------- |
219: | | | |
220: n | yn | | |
221: | | | |
222: | .--------------------- |
223: | (xs,ys) xn |
224: | . |
225: | (gxs,gys) |
226: | |
227: -----------------------------------------------------
228: */
230: /*
231: nc - number of components per grid point
232: col - number of colors needed in one direction for single component problem
234: */
235: DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);
237: PetscObjectGetComm((PetscObject)da,&comm);
238: MPI_Comm_size(comm,&size);
239: if (ctype == IS_COLORING_LOCAL) {
240: if (size == 1) {
241: ctype = IS_COLORING_GLOBAL;
242: } else if (dim > 1) {
243: if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
244: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_LOCAL cannot be used for periodic boundary condition having both ends of the domain on the same process");
245: }
246: }
247: }
249: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
250: matrices is for the blocks, not the individual matrix elements */
251: PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);
252: if (!isBAIJ) {PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);}
253: if (!isBAIJ) {PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);}
254: if (isBAIJ) {
255: dd->w = 1;
256: dd->xs = dd->xs/nc;
257: dd->xe = dd->xe/nc;
258: dd->Xs = dd->Xs/nc;
259: dd->Xe = dd->Xe/nc;
260: }
262: /*
263: We do not provide a getcoloring function in the DMDA operations because
264: the basic DMDA does not know about matrices. We think of DMDA as being more
265: more low-level then matrices.
266: */
267: if (dim == 1) {
268: DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
269: } else if (dim == 2) {
270: DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
271: } else if (dim == 3) {
272: DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
273: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
274: if (isBAIJ) {
275: dd->w = nc;
276: dd->xs = dd->xs*nc;
277: dd->xe = dd->xe*nc;
278: dd->Xs = dd->Xs*nc;
279: dd->Xe = dd->Xe*nc;
280: }
281: return(0);
282: }
284: /* ---------------------------------------------------------------------------------*/
286: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
287: {
288: PetscErrorCode ierr;
289: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
290: PetscInt ncolors;
291: MPI_Comm comm;
292: DMBoundaryType bx,by;
293: DMDAStencilType st;
294: ISColoringValue *colors;
295: DM_DA *dd = (DM_DA*)da->data;
298: /*
299: nc - number of components per grid point
300: col - number of colors needed in one direction for single component problem
302: */
303: DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
304: col = 2*s + 1;
305: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
306: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
307: PetscObjectGetComm((PetscObject)da,&comm);
309: /* special case as taught to us by Paul Hovland */
310: if (st == DMDA_STENCIL_STAR && s == 1) {
311: DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
312: } else {
314: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\
315: by 2*stencil_width + 1 (%d)\n", m, col);
316: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\
317: by 2*stencil_width + 1 (%d)\n", n, col);
318: if (ctype == IS_COLORING_GLOBAL) {
319: if (!dd->localcoloring) {
320: PetscMalloc1(nc*nx*ny,&colors);
321: ii = 0;
322: for (j=ys; j<ys+ny; j++) {
323: for (i=xs; i<xs+nx; i++) {
324: for (k=0; k<nc; k++) {
325: colors[ii++] = k + nc*((i % col) + col*(j % col));
326: }
327: }
328: }
329: ncolors = nc + nc*(col-1 + col*(col-1));
330: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
331: }
332: *coloring = dd->localcoloring;
333: } else if (ctype == IS_COLORING_LOCAL) {
334: if (!dd->ghostedcoloring) {
335: PetscMalloc1(nc*gnx*gny,&colors);
336: ii = 0;
337: for (j=gys; j<gys+gny; j++) {
338: for (i=gxs; i<gxs+gnx; i++) {
339: for (k=0; k<nc; k++) {
340: /* the complicated stuff is to handle periodic boundaries */
341: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
342: }
343: }
344: }
345: ncolors = nc + nc*(col - 1 + col*(col-1));
346: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
347: /* PetscIntView(ncolors,(PetscInt*)colors,0); */
349: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
350: }
351: *coloring = dd->ghostedcoloring;
352: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
353: }
354: ISColoringReference(*coloring);
355: return(0);
356: }
358: /* ---------------------------------------------------------------------------------*/
360: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
361: {
362: PetscErrorCode ierr;
363: 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;
364: PetscInt ncolors;
365: MPI_Comm comm;
366: DMBoundaryType bx,by,bz;
367: DMDAStencilType st;
368: ISColoringValue *colors;
369: DM_DA *dd = (DM_DA*)da->data;
372: /*
373: nc - number of components per grid point
374: col - number of colors needed in one direction for single component problem
376: */
377: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
378: col = 2*s + 1;
379: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
380: by 2*stencil_width + 1\n");
381: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
382: by 2*stencil_width + 1\n");
383: if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
384: by 2*stencil_width + 1\n");
386: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
387: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
388: PetscObjectGetComm((PetscObject)da,&comm);
390: /* create the coloring */
391: if (ctype == IS_COLORING_GLOBAL) {
392: if (!dd->localcoloring) {
393: PetscMalloc1(nc*nx*ny*nz,&colors);
394: ii = 0;
395: for (k=zs; k<zs+nz; k++) {
396: for (j=ys; j<ys+ny; j++) {
397: for (i=xs; i<xs+nx; i++) {
398: for (l=0; l<nc; l++) {
399: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
400: }
401: }
402: }
403: }
404: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
405: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
406: }
407: *coloring = dd->localcoloring;
408: } else if (ctype == IS_COLORING_LOCAL) {
409: if (!dd->ghostedcoloring) {
410: PetscMalloc1(nc*gnx*gny*gnz,&colors);
411: ii = 0;
412: for (k=gzs; k<gzs+gnz; k++) {
413: for (j=gys; j<gys+gny; j++) {
414: for (i=gxs; i<gxs+gnx; i++) {
415: for (l=0; l<nc; l++) {
416: /* the complicated stuff is to handle periodic boundaries */
417: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
418: }
419: }
420: }
421: }
422: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
423: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
424: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
425: }
426: *coloring = dd->ghostedcoloring;
427: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
428: ISColoringReference(*coloring);
429: return(0);
430: }
432: /* ---------------------------------------------------------------------------------*/
434: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
435: {
436: PetscErrorCode ierr;
437: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
438: PetscInt ncolors;
439: MPI_Comm comm;
440: DMBoundaryType bx;
441: ISColoringValue *colors;
442: DM_DA *dd = (DM_DA*)da->data;
445: /*
446: nc - number of components per grid point
447: col - number of colors needed in one direction for single component problem
449: */
450: DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
451: col = 2*s + 1;
453: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\
454: by 2*stencil_width + 1 %d\n",(int)m,(int)col);
456: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
457: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
458: PetscObjectGetComm((PetscObject)da,&comm);
460: /* create the coloring */
461: if (ctype == IS_COLORING_GLOBAL) {
462: if (!dd->localcoloring) {
463: PetscMalloc1(nc*nx,&colors);
464: if (dd->ofillcols) {
465: PetscInt tc = 0;
466: for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
467: i1 = 0;
468: for (i=xs; i<xs+nx; i++) {
469: for (l=0; l<nc; l++) {
470: if (dd->ofillcols[l] && (i % col)) {
471: colors[i1++] = nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
472: } else {
473: colors[i1++] = l;
474: }
475: }
476: }
477: ncolors = nc + 2*s*tc;
478: } else {
479: i1 = 0;
480: for (i=xs; i<xs+nx; i++) {
481: for (l=0; l<nc; l++) {
482: colors[i1++] = l + nc*(i % col);
483: }
484: }
485: ncolors = nc + nc*(col-1);
486: }
487: ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
488: }
489: *coloring = dd->localcoloring;
490: } else if (ctype == IS_COLORING_LOCAL) {
491: if (!dd->ghostedcoloring) {
492: PetscMalloc1(nc*gnx,&colors);
493: i1 = 0;
494: for (i=gxs; i<gxs+gnx; i++) {
495: for (l=0; l<nc; l++) {
496: /* the complicated stuff is to handle periodic boundaries */
497: colors[i1++] = l + nc*(SetInRange(i,m) % col);
498: }
499: }
500: ncolors = nc + nc*(col-1);
501: ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
502: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
503: }
504: *coloring = dd->ghostedcoloring;
505: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
506: ISColoringReference(*coloring);
507: return(0);
508: }
510: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
511: {
512: PetscErrorCode ierr;
513: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
514: PetscInt ncolors;
515: MPI_Comm comm;
516: DMBoundaryType bx,by;
517: ISColoringValue *colors;
518: DM_DA *dd = (DM_DA*)da->data;
521: /*
522: nc - number of components per grid point
523: col - number of colors needed in one direction for single component problem
525: */
526: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
527: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
528: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
529: PetscObjectGetComm((PetscObject)da,&comm);
531: if (bx == DM_BOUNDARY_PERIODIC && (m % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible by 5\n");
532: if (by == DM_BOUNDARY_PERIODIC && (n % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible by 5\n");
534: /* create the coloring */
535: if (ctype == IS_COLORING_GLOBAL) {
536: if (!dd->localcoloring) {
537: PetscMalloc1(nc*nx*ny,&colors);
538: ii = 0;
539: for (j=ys; j<ys+ny; j++) {
540: for (i=xs; i<xs+nx; i++) {
541: for (k=0; k<nc; k++) {
542: colors[ii++] = k + nc*((3*j+i) % 5);
543: }
544: }
545: }
546: ncolors = 5*nc;
547: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
548: }
549: *coloring = dd->localcoloring;
550: } else if (ctype == IS_COLORING_LOCAL) {
551: if (!dd->ghostedcoloring) {
552: PetscMalloc1(nc*gnx*gny,&colors);
553: ii = 0;
554: for (j=gys; j<gys+gny; j++) {
555: for (i=gxs; i<gxs+gnx; i++) {
556: for (k=0; k<nc; k++) {
557: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
558: }
559: }
560: }
561: ncolors = 5*nc;
562: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
563: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
564: }
565: *coloring = dd->ghostedcoloring;
566: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
567: return(0);
568: }
570: /* =========================================================================== */
571: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
572: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
573: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
574: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
575: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
576: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
577: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
578: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
579: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
580: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
581: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM,Mat);
582: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM,Mat);
583: extern PetscErrorCode DMCreateMatrix_DA_IS(DM,Mat);
585: /*@C
586: MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix
588: Logically Collective on Mat
590: Input Parameters:
591: + mat - the matrix
592: - da - the da
594: Level: intermediate
596: @*/
597: PetscErrorCode MatSetupDM(Mat mat,DM da)
598: {
604: PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
605: return(0);
606: }
608: PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer)
609: {
610: DM da;
611: PetscErrorCode ierr;
612: const char *prefix;
613: Mat Anatural;
614: AO ao;
615: PetscInt rstart,rend,*petsc,i;
616: IS is;
617: MPI_Comm comm;
618: PetscViewerFormat format;
621: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
622: PetscViewerGetFormat(viewer,&format);
623: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
625: PetscObjectGetComm((PetscObject)A,&comm);
626: MatGetDM(A, &da);
627: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
629: DMDAGetAO(da,&ao);
630: MatGetOwnershipRange(A,&rstart,&rend);
631: PetscMalloc1(rend-rstart,&petsc);
632: for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
633: AOApplicationToPetsc(ao,rend-rstart,petsc);
634: ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);
636: /* call viewer on natural ordering */
637: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
638: ISDestroy(&is);
639: PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
640: PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
641: PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
642: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
643: MatView(Anatural,viewer);
644: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
645: MatDestroy(&Anatural);
646: return(0);
647: }
649: PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer)
650: {
651: DM da;
653: Mat Anatural,Aapp;
654: AO ao;
655: PetscInt rstart,rend,*app,i,m,n,M,N;
656: IS is;
657: MPI_Comm comm;
660: PetscObjectGetComm((PetscObject)A,&comm);
661: MatGetDM(A, &da);
662: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
664: /* Load the matrix in natural ordering */
665: MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
666: MatSetType(Anatural,((PetscObject)A)->type_name);
667: MatGetSize(A,&M,&N);
668: MatGetLocalSize(A,&m,&n);
669: MatSetSizes(Anatural,m,n,M,N);
670: MatLoad(Anatural,viewer);
672: /* Map natural ordering to application ordering and create IS */
673: DMDAGetAO(da,&ao);
674: MatGetOwnershipRange(Anatural,&rstart,&rend);
675: PetscMalloc1(rend-rstart,&app);
676: for (i=rstart; i<rend; i++) app[i-rstart] = i;
677: AOPetscToApplication(ao,rend-rstart,app);
678: ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);
680: /* Do permutation and replace header */
681: MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
682: MatHeaderReplace(A,&Aapp);
683: ISDestroy(&is);
684: MatDestroy(&Anatural);
685: return(0);
686: }
688: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
689: {
691: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
692: Mat A;
693: MPI_Comm comm;
694: MatType Atype;
695: PetscSection section, sectionGlobal;
696: void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL,(*sell)(void)=NULL,(*is)(void)=NULL;
697: MatType mtype;
698: PetscMPIInt size;
699: DM_DA *dd = (DM_DA*)da->data;
702: MatInitializePackage();
703: mtype = da->mattype;
705: DMGetSection(da, §ion);
706: if (section) {
707: PetscInt bs = -1;
708: PetscInt localSize;
709: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
711: DMGetGlobalSection(da, §ionGlobal);
712: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
713: MatCreate(PetscObjectComm((PetscObject)da),&A);
714: MatSetSizes(A,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
715: MatSetType(A,mtype);
716: PetscStrcmp(mtype,MATSHELL,&isShell);
717: PetscStrcmp(mtype,MATBAIJ,&isBlock);
718: PetscStrcmp(mtype,MATSEQBAIJ,&isSeqBlock);
719: PetscStrcmp(mtype,MATMPIBAIJ,&isMPIBlock);
720: PetscStrcmp(mtype,MATSBAIJ,&isSymBlock);
721: PetscStrcmp(mtype,MATSEQSBAIJ,&isSymSeqBlock);
722: PetscStrcmp(mtype,MATMPISBAIJ,&isSymMPIBlock);
723: /* Check for symmetric storage */
724: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
725: if (isSymmetric) {
726: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
727: }
728: if (!isShell) {
729: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
731: if (bs < 0) {
732: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
733: PetscInt pStart, pEnd, p, dof;
735: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
736: for (p = pStart; p < pEnd; ++p) {
737: PetscSectionGetDof(sectionGlobal, p, &dof);
738: if (dof) {
739: bs = dof;
740: break;
741: }
742: }
743: } else {
744: bs = 1;
745: }
746: /* Must have same blocksize on all procs (some might have no points) */
747: bsLocal = bs;
748: MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
749: }
750: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
751: /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
752: PetscFree4(dnz, onz, dnzu, onzu);
753: }
754: }
755: /*
756: m
757: ------------------------------------------------------
758: | |
759: | |
760: | ---------------------- |
761: | | | |
762: n | ny | | |
763: | | | |
764: | .--------------------- |
765: | (xs,ys) nx |
766: | . |
767: | (gxs,gys) |
768: | |
769: -----------------------------------------------------
770: */
772: /*
773: nc - number of components per grid point
774: col - number of colors needed in one direction for single component problem
776: */
777: M = dd->M;
778: N = dd->N;
779: P = dd->P;
780: dim = da->dim;
781: dof = dd->w;
782: /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
783: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
784: PetscObjectGetComm((PetscObject)da,&comm);
785: MatCreate(comm,&A);
786: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
787: MatSetType(A,mtype);
788: MatSetFromOptions(A);
789: MatSetDM(A,da);
790: if (da->structure_only) {
791: MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);
792: }
793: MatGetType(A,&Atype);
794: /*
795: We do not provide a getmatrix function in the DMDA operations because
796: the basic DMDA does not know about matrices. We think of DMDA as being more
797: more low-level than matrices. This is kind of cheating but, cause sometimes
798: we think of DMDA has higher level than matrices.
800: We could switch based on Atype (or mtype), but we do not since the
801: specialized setting routines depend only on the particular preallocation
802: details of the matrix, not the type itself.
803: */
804: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
805: if (!aij) {
806: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
807: }
808: if (!aij) {
809: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
810: if (!baij) {
811: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
812: }
813: if (!baij) {
814: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
815: if (!sbaij) {
816: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
817: }
818: if (!sbaij) {
819: PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);
820: if (!sell) {
821: PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);
822: }
823: }
824: if (!sell) {
825: PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
826: }
827: }
828: }
829: if (aij) {
830: if (dim == 1) {
831: if (dd->ofill) {
832: DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
833: } else {
834: DMCreateMatrix_DA_1d_MPIAIJ(da,A);
835: }
836: } else if (dim == 2) {
837: if (dd->ofill) {
838: DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
839: } else {
840: DMCreateMatrix_DA_2d_MPIAIJ(da,A);
841: }
842: } else if (dim == 3) {
843: if (dd->ofill) {
844: DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
845: } else {
846: DMCreateMatrix_DA_3d_MPIAIJ(da,A);
847: }
848: }
849: } else if (baij) {
850: if (dim == 2) {
851: DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
852: } else if (dim == 3) {
853: DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
854: } else SETERRQ3(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim);
855: } else if (sbaij) {
856: if (dim == 2) {
857: DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
858: } else if (dim == 3) {
859: DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
860: } else SETERRQ3(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim);
861: } else if (sell) {
862: if (dim == 2) {
863: DMCreateMatrix_DA_2d_MPISELL(da,A);
864: } else if (dim == 3) {
865: DMCreateMatrix_DA_3d_MPISELL(da,A);
866: } else SETERRQ3(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim);
867: } else if (is) {
868: DMCreateMatrix_DA_IS(da,A);
869: } else {
870: ISLocalToGlobalMapping ltog;
872: MatSetBlockSize(A,dof);
873: MatSetUp(A);
874: DMGetLocalToGlobalMapping(da,<og);
875: MatSetLocalToGlobalMapping(A,ltog,ltog);
876: }
877: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
878: MatSetStencil(A,dim,dims,starts,dof);
879: MatSetDM(A,da);
880: MPI_Comm_size(comm,&size);
881: if (size > 1) {
882: /* change viewer to display matrix in natural ordering */
883: MatSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
884: MatSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
885: }
886: *J = A;
887: return(0);
888: }
890: /* ---------------------------------------------------------------------------------*/
891: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]);
893: PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J)
894: {
895: DM_DA *da = (DM_DA*)dm->data;
896: Mat lJ;
897: ISLocalToGlobalMapping ltog;
898: IS is_loc_filt, is_glob;
899: const PetscInt *e_loc,*idx;
900: PetscInt nel,nen,nv,dof,dim,*gidx,nb;
901: PetscBool flg;
902: PetscErrorCode ierr;
904: /* The l2g map of DMDA has all ghosted nodes, and e_loc is a subset of all the local nodes (including the ghosted)
905: We need to filter the local indices that are represented through the DMDAGetElements decomposition
906: This is because the size of the local matrices in MATIS is the local size of the l2g map */
908: dof = da->w;
909: dim = dm->dim;
911: MatSetBlockSize(J,dof);
913: /* get local elements indices in local DMDA numbering */
914: DMDAGetElements(dm,&nel,&nen,&e_loc); /* this will throw an error if the stencil type is not DMDA_STENCIL_BOX */
915: ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nel*nen,e_loc,PETSC_COPY_VALUES,&is_loc_filt);
916: DMDARestoreElements(dm,&nel,&nen,&e_loc);
918: /* obtain a consistent local ordering for MATIS */
919: ISSortRemoveDups(is_loc_filt);
920: ISBlockGetLocalSize(is_loc_filt,&nb);
921: DMGetLocalToGlobalMapping(dm,<og);
922: ISLocalToGlobalMappingGetSize(ltog,&nv);
923: PetscMalloc1(PetscMax(nb,nv/dof),&gidx);
924: ISBlockGetIndices(is_loc_filt,&idx);
925: ISLocalToGlobalMappingApplyBlock(ltog,nb,idx,gidx);
926: ISBlockRestoreIndices(is_loc_filt,&idx);
927: ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nb,gidx,PETSC_USE_POINTER,&is_glob);
928: ISLocalToGlobalMappingCreateIS(is_glob,<og);
929: ISDestroy(&is_glob);
930: MatSetLocalToGlobalMapping(J,ltog,ltog);
931: ISLocalToGlobalMappingDestroy(<og);
933: /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */
934: MatISGetLocalMat(J,&lJ);
935: ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);
936: ISDestroy(&is_loc_filt);
937: ISCreateStride(PetscObjectComm((PetscObject)lJ),nv/dof,0,1,&is_glob);
938: ISGetIndices(is_glob,&idx);
939: ISGlobalToLocalMappingApplyBlock(ltog,IS_GTOLM_MASK,nv/dof,idx,&nb,gidx);
940: ISRestoreIndices(is_glob,&idx);
941: ISDestroy(&is_glob);
942: ISLocalToGlobalMappingDestroy(<og);
943: ISCreateBlock(PETSC_COMM_SELF,dof,nb,gidx,PETSC_USE_POINTER,&is_loc_filt);
944: ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);
945: ISDestroy(&is_loc_filt);
946: MatSetLocalToGlobalMapping(lJ,ltog,ltog);
947: ISLocalToGlobalMappingDestroy(<og);
948: PetscFree(gidx);
950: /* Preallocation (not exact): we reuse the preallocation routines of the assembled version */
951: flg = dm->prealloc_only;
952: dm->prealloc_only = PETSC_TRUE;
953: switch (dim) {
954: case 1:
955: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
956: DMCreateMatrix_DA_1d_MPIAIJ(dm,J);
957: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
958: break;
959: case 2:
960: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
961: DMCreateMatrix_DA_2d_MPIAIJ(dm,J);
962: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
963: break;
964: case 3:
965: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",MatISSetPreallocation_IS);
966: DMCreateMatrix_DA_3d_MPIAIJ(dm,J);
967: PetscObjectComposeFunction((PetscObject)J,"MatMPIAIJSetPreallocation_C",NULL);
968: break;
969: default:
970: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled dimension %d",dim);
971: break;
972: }
973: dm->prealloc_only = flg;
974: return(0);
975: }
977: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J)
978: {
979: PetscErrorCode ierr;
980: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p;
981: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
982: MPI_Comm comm;
983: PetscScalar *values;
984: DMBoundaryType bx,by;
985: ISLocalToGlobalMapping ltog;
986: DMDAStencilType st;
989: /*
990: nc - number of components per grid point
991: col - number of colors needed in one direction for single component problem
993: */
994: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
995: col = 2*s + 1;
996: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
997: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
998: PetscObjectGetComm((PetscObject)da,&comm);
1000: PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1001: DMGetLocalToGlobalMapping(da,<og);
1003: MatSetBlockSize(J,nc);
1004: /* determine the matrix preallocation information */
1005: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1006: for (i=xs; i<xs+nx; i++) {
1008: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1009: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1011: for (j=ys; j<ys+ny; j++) {
1012: slot = i - gxs + gnx*(j - gys);
1014: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1015: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1017: cnt = 0;
1018: for (k=0; k<nc; k++) {
1019: for (l=lstart; l<lend+1; l++) {
1020: for (p=pstart; p<pend+1; p++) {
1021: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1022: cols[cnt++] = k + nc*(slot + gnx*l + p);
1023: }
1024: }
1025: }
1026: rows[k] = k + nc*(slot);
1027: }
1028: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1029: }
1030: }
1031: MatSetBlockSize(J,nc);
1032: MatSeqSELLSetPreallocation(J,0,dnz);
1033: MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1034: MatPreallocateFinalize(dnz,onz);
1036: MatSetLocalToGlobalMapping(J,ltog,ltog);
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: PetscCalloc1(col*col*nc*nc,&values);
1045: for (i=xs; i<xs+nx; i++) {
1047: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1048: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1050: for (j=ys; j<ys+ny; j++) {
1051: slot = i - gxs + gnx*(j - gys);
1053: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1054: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1056: cnt = 0;
1057: for (k=0; k<nc; k++) {
1058: for (l=lstart; l<lend+1; l++) {
1059: for (p=pstart; p<pend+1; p++) {
1060: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1061: cols[cnt++] = k + nc*(slot + gnx*l + p);
1062: }
1063: }
1064: }
1065: rows[k] = k + nc*(slot);
1066: }
1067: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1068: }
1069: }
1070: PetscFree(values);
1071: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1072: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1073: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1074: }
1075: PetscFree2(rows,cols);
1076: return(0);
1077: }
1079: PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da,Mat J)
1080: {
1081: PetscErrorCode ierr;
1082: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1083: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1084: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1085: MPI_Comm comm;
1086: PetscScalar *values;
1087: DMBoundaryType bx,by,bz;
1088: ISLocalToGlobalMapping ltog;
1089: DMDAStencilType st;
1092: /*
1093: nc - number of components per grid point
1094: col - number of colors needed in one direction for single component problem
1096: */
1097: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1098: col = 2*s + 1;
1099: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1100: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1101: PetscObjectGetComm((PetscObject)da,&comm);
1103: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1104: DMGetLocalToGlobalMapping(da,<og);
1106: MatSetBlockSize(J,nc);
1107: /* determine the matrix preallocation information */
1108: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1109: for (i=xs; i<xs+nx; i++) {
1110: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1111: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1112: for (j=ys; j<ys+ny; j++) {
1113: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1114: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1115: for (k=zs; k<zs+nz; k++) {
1116: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1117: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1119: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1121: cnt = 0;
1122: for (l=0; l<nc; l++) {
1123: for (ii=istart; ii<iend+1; ii++) {
1124: for (jj=jstart; jj<jend+1; jj++) {
1125: for (kk=kstart; kk<kend+1; kk++) {
1126: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1127: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1128: }
1129: }
1130: }
1131: }
1132: rows[l] = l + nc*(slot);
1133: }
1134: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1135: }
1136: }
1137: }
1138: MatSetBlockSize(J,nc);
1139: MatSeqSELLSetPreallocation(J,0,dnz);
1140: MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1141: MatPreallocateFinalize(dnz,onz);
1142: MatSetLocalToGlobalMapping(J,ltog,ltog);
1144: /*
1145: For each node in the grid: we get the neighbors in the local (on processor ordering
1146: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1147: PETSc ordering.
1148: */
1149: if (!da->prealloc_only) {
1150: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1151: for (i=xs; i<xs+nx; i++) {
1152: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1153: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1154: for (j=ys; j<ys+ny; j++) {
1155: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1156: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1157: for (k=zs; k<zs+nz; k++) {
1158: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1159: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1161: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1163: cnt = 0;
1164: for (l=0; l<nc; l++) {
1165: for (ii=istart; ii<iend+1; ii++) {
1166: for (jj=jstart; jj<jend+1; jj++) {
1167: for (kk=kstart; kk<kend+1; kk++) {
1168: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1169: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1170: }
1171: }
1172: }
1173: }
1174: rows[l] = l + nc*(slot);
1175: }
1176: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1177: }
1178: }
1179: }
1180: PetscFree(values);
1181: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1182: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1183: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1184: }
1185: PetscFree2(rows,cols);
1186: return(0);
1187: }
1189: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
1190: {
1191: PetscErrorCode ierr;
1192: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,M,N;
1193: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
1194: MPI_Comm comm;
1195: PetscScalar *values;
1196: DMBoundaryType bx,by;
1197: ISLocalToGlobalMapping ltog,mltog;
1198: DMDAStencilType st;
1199: PetscBool removedups = PETSC_FALSE;
1202: /*
1203: nc - number of components per grid point
1204: col - number of colors needed in one direction for single component problem
1206: */
1207: DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1208: col = 2*s + 1;
1209: /*
1210: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1211: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1212: */
1213: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1214: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1215: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1216: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1217: PetscObjectGetComm((PetscObject)da,&comm);
1219: PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1220: DMGetLocalToGlobalMapping(da,<og);
1222: MatSetBlockSize(J,nc);
1223: /* determine the matrix preallocation information */
1224: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1225: for (i=xs; i<xs+nx; i++) {
1227: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1228: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1230: for (j=ys; j<ys+ny; j++) {
1231: slot = i - gxs + gnx*(j - gys);
1233: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1234: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1236: cnt = 0;
1237: for (k=0; k<nc; k++) {
1238: for (l=lstart; l<lend+1; l++) {
1239: for (p=pstart; p<pend+1; p++) {
1240: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1241: cols[cnt++] = k + nc*(slot + gnx*l + p);
1242: }
1243: }
1244: }
1245: rows[k] = k + nc*(slot);
1246: }
1247: if (removedups) {
1248: MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1249: } else {
1250: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1251: }
1252: }
1253: }
1254: MatSetBlockSize(J,nc);
1255: MatSeqAIJSetPreallocation(J,0,dnz);
1256: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1257: MatPreallocateFinalize(dnz,onz);
1258: MatGetLocalToGlobalMapping(J,&mltog,NULL);
1259: if (!mltog) {
1260: MatSetLocalToGlobalMapping(J,ltog,ltog);
1261: }
1263: /*
1264: For each node in the grid: we get the neighbors in the local (on processor ordering
1265: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1266: PETSc ordering.
1267: */
1268: if (!da->prealloc_only) {
1269: PetscCalloc1(col*col*nc*nc,&values);
1270: for (i=xs; i<xs+nx; i++) {
1272: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1273: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1275: for (j=ys; j<ys+ny; j++) {
1276: slot = i - gxs + gnx*(j - gys);
1278: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1279: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1281: cnt = 0;
1282: for (k=0; k<nc; k++) {
1283: for (l=lstart; l<lend+1; l++) {
1284: for (p=pstart; p<pend+1; p++) {
1285: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1286: cols[cnt++] = k + nc*(slot + gnx*l + p);
1287: }
1288: }
1289: }
1290: rows[k] = k + nc*(slot);
1291: }
1292: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1293: }
1294: }
1295: PetscFree(values);
1296: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1297: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1298: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1299: }
1300: PetscFree2(rows,cols);
1301: return(0);
1302: }
1304: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
1305: {
1306: PetscErrorCode ierr;
1307: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1308: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
1309: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
1310: DM_DA *dd = (DM_DA*)da->data;
1311: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
1312: MPI_Comm comm;
1313: PetscScalar *values;
1314: DMBoundaryType bx,by;
1315: ISLocalToGlobalMapping ltog;
1316: DMDAStencilType st;
1317: PetscBool removedups = PETSC_FALSE;
1320: /*
1321: nc - number of components per grid point
1322: col - number of colors needed in one direction for single component problem
1324: */
1325: DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1326: col = 2*s + 1;
1327: /*
1328: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1329: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1330: */
1331: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1332: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1333: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1334: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1335: PetscObjectGetComm((PetscObject)da,&comm);
1337: PetscMalloc1(col*col*nc,&cols);
1338: DMGetLocalToGlobalMapping(da,<og);
1340: MatSetBlockSize(J,nc);
1341: /* determine the matrix preallocation information */
1342: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1343: for (i=xs; i<xs+nx; i++) {
1345: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1346: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1348: for (j=ys; j<ys+ny; j++) {
1349: slot = i - gxs + gnx*(j - gys);
1351: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1352: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1354: for (k=0; k<nc; k++) {
1355: cnt = 0;
1356: for (l=lstart; l<lend+1; l++) {
1357: for (p=pstart; p<pend+1; p++) {
1358: if (l || p) {
1359: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1360: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1361: }
1362: } else {
1363: if (dfill) {
1364: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1365: } else {
1366: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1367: }
1368: }
1369: }
1370: }
1371: row = k + nc*(slot);
1372: maxcnt = PetscMax(maxcnt,cnt);
1373: if (removedups) {
1374: MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1375: } else {
1376: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1377: }
1378: }
1379: }
1380: }
1381: MatSeqAIJSetPreallocation(J,0,dnz);
1382: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1383: MatPreallocateFinalize(dnz,onz);
1384: MatSetLocalToGlobalMapping(J,ltog,ltog);
1386: /*
1387: For each node in the grid: we get the neighbors in the local (on processor ordering
1388: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1389: PETSc ordering.
1390: */
1391: if (!da->prealloc_only) {
1392: PetscCalloc1(maxcnt,&values);
1393: for (i=xs; i<xs+nx; i++) {
1395: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1396: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1398: for (j=ys; j<ys+ny; j++) {
1399: slot = i - gxs + gnx*(j - gys);
1401: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1402: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1404: for (k=0; k<nc; k++) {
1405: cnt = 0;
1406: for (l=lstart; l<lend+1; l++) {
1407: for (p=pstart; p<pend+1; p++) {
1408: if (l || p) {
1409: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1410: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1411: }
1412: } else {
1413: if (dfill) {
1414: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1415: } else {
1416: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1417: }
1418: }
1419: }
1420: }
1421: row = k + nc*(slot);
1422: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1423: }
1424: }
1425: }
1426: PetscFree(values);
1427: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1428: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1429: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1430: }
1431: PetscFree(cols);
1432: return(0);
1433: }
1435: /* ---------------------------------------------------------------------------------*/
1437: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1438: {
1439: PetscErrorCode ierr;
1440: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1441: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1442: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1443: MPI_Comm comm;
1444: PetscScalar *values;
1445: DMBoundaryType bx,by,bz;
1446: ISLocalToGlobalMapping ltog,mltog;
1447: DMDAStencilType st;
1448: PetscBool removedups = PETSC_FALSE;
1451: /*
1452: nc - number of components per grid point
1453: col - number of colors needed in one direction for single component problem
1455: */
1456: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1457: col = 2*s + 1;
1459: /*
1460: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1461: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1462: */
1463: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1464: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1465: if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;
1467: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1468: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1469: PetscObjectGetComm((PetscObject)da,&comm);
1471: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1472: DMGetLocalToGlobalMapping(da,<og);
1474: MatSetBlockSize(J,nc);
1475: /* determine the matrix preallocation information */
1476: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1477: for (i=xs; i<xs+nx; i++) {
1478: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1479: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1480: for (j=ys; j<ys+ny; j++) {
1481: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1482: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1483: for (k=zs; k<zs+nz; k++) {
1484: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1485: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1487: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1489: cnt = 0;
1490: for (l=0; l<nc; l++) {
1491: for (ii=istart; ii<iend+1; ii++) {
1492: for (jj=jstart; jj<jend+1; jj++) {
1493: for (kk=kstart; kk<kend+1; kk++) {
1494: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1495: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1496: }
1497: }
1498: }
1499: }
1500: rows[l] = l + nc*(slot);
1501: }
1502: if (removedups) {
1503: MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1504: } else {
1505: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1506: }
1507: }
1508: }
1509: }
1510: MatSetBlockSize(J,nc);
1511: MatSeqAIJSetPreallocation(J,0,dnz);
1512: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1513: MatPreallocateFinalize(dnz,onz);
1514: MatGetLocalToGlobalMapping(J,&mltog,NULL);
1515: if (!mltog) {
1516: MatSetLocalToGlobalMapping(J,ltog,ltog);
1517: }
1519: /*
1520: For each node in the grid: we get the neighbors in the local (on processor ordering
1521: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1522: PETSc ordering.
1523: */
1524: if (!da->prealloc_only) {
1525: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1526: for (i=xs; i<xs+nx; i++) {
1527: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1528: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1529: for (j=ys; j<ys+ny; j++) {
1530: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1531: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1532: for (k=zs; k<zs+nz; k++) {
1533: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1534: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1536: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1538: cnt = 0;
1539: for (l=0; l<nc; l++) {
1540: for (ii=istart; ii<iend+1; ii++) {
1541: for (jj=jstart; jj<jend+1; jj++) {
1542: for (kk=kstart; kk<kend+1; kk++) {
1543: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1544: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1545: }
1546: }
1547: }
1548: }
1549: rows[l] = l + nc*(slot);
1550: }
1551: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1552: }
1553: }
1554: }
1555: PetscFree(values);
1556: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1557: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1558: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1559: }
1560: PetscFree2(rows,cols);
1561: return(0);
1562: }
1564: /* ---------------------------------------------------------------------------------*/
1566: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1567: {
1568: PetscErrorCode ierr;
1569: DM_DA *dd = (DM_DA*)da->data;
1570: PetscInt xs,nx,i,j,gxs,gnx,row,k,l;
1571: PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1572: PetscInt *ofill = dd->ofill,*dfill = dd->dfill;
1573: PetscScalar *values;
1574: DMBoundaryType bx;
1575: ISLocalToGlobalMapping ltog;
1576: PetscMPIInt rank,size;
1579: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1580: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1582: /*
1583: nc - number of components per grid point
1585: */
1586: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1587: if (s > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Matrix creation for 1d not implemented correctly for stencil width larger than 1");
1588: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1589: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1591: MatSetBlockSize(J,nc);
1592: PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);
1594: /*
1595: note should be smaller for first and last process with no periodic
1596: does not handle dfill
1597: */
1598: cnt = 0;
1599: /* coupling with process to the left */
1600: for (i=0; i<s; i++) {
1601: for (j=0; j<nc; j++) {
1602: ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1603: cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1604: if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1605: if (size > 1) ocols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1606: else cols[cnt] += (s - i)*(ofill[j+1] - ofill[j]);
1607: }
1608: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1609: cnt++;
1610: }
1611: }
1612: for (i=s; i<nx-s; i++) {
1613: for (j=0; j<nc; j++) {
1614: cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1615: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1616: cnt++;
1617: }
1618: }
1619: /* coupling with process to the right */
1620: for (i=nx-s; i<nx; i++) {
1621: for (j=0; j<nc; j++) {
1622: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1623: cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1624: if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1625: if (size > 1) ocols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1626: else cols[cnt] += (i - nx + s + 1)*(ofill[j+1] - ofill[j]);
1627: }
1628: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1629: cnt++;
1630: }
1631: }
1633: MatSeqAIJSetPreallocation(J,0,cols);
1634: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1635: PetscFree2(cols,ocols);
1637: DMGetLocalToGlobalMapping(da,<og);
1638: MatSetLocalToGlobalMapping(J,ltog,ltog);
1640: /*
1641: For each node in the grid: we get the neighbors in the local (on processor ordering
1642: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1643: PETSc ordering.
1644: */
1645: if (!da->prealloc_only) {
1646: PetscCalloc2(maxcnt,&values,maxcnt,&cols);
1648: row = xs*nc;
1649: /* coupling with process to the left */
1650: for (i=xs; i<xs+s; i++) {
1651: for (j=0; j<nc; j++) {
1652: cnt = 0;
1653: if (rank) {
1654: for (l=0; l<s; l++) {
1655: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1656: }
1657: }
1658: if (!rank && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1659: for (l=0; l<s; l++) {
1660: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (m + i - s - l)*nc + ofill[k];
1661: }
1662: }
1663: if (dfill) {
1664: for (k=dfill[j]; k<dfill[j+1]; k++) {
1665: cols[cnt++] = i*nc + dfill[k];
1666: }
1667: } else {
1668: for (k=0; k<nc; k++) {
1669: cols[cnt++] = i*nc + k;
1670: }
1671: }
1672: for (l=0; l<s; l++) {
1673: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1674: }
1675: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1676: row++;
1677: }
1678: }
1679: for (i=xs+s; i<xs+nx-s; i++) {
1680: for (j=0; j<nc; j++) {
1681: cnt = 0;
1682: for (l=0; l<s; l++) {
1683: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1684: }
1685: if (dfill) {
1686: for (k=dfill[j]; k<dfill[j+1]; k++) {
1687: cols[cnt++] = i*nc + dfill[k];
1688: }
1689: } else {
1690: for (k=0; k<nc; k++) {
1691: cols[cnt++] = i*nc + k;
1692: }
1693: }
1694: for (l=0; l<s; l++) {
1695: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1696: }
1697: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1698: row++;
1699: }
1700: }
1701: /* coupling with process to the right */
1702: for (i=xs+nx-s; i<xs+nx; i++) {
1703: for (j=0; j<nc; j++) {
1704: cnt = 0;
1705: for (l=0; l<s; l++) {
1706: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1707: }
1708: if (dfill) {
1709: for (k=dfill[j]; k<dfill[j+1]; k++) {
1710: cols[cnt++] = i*nc + dfill[k];
1711: }
1712: } else {
1713: for (k=0; k<nc; k++) {
1714: cols[cnt++] = i*nc + k;
1715: }
1716: }
1717: if (rank < size-1) {
1718: for (l=0; l<s; l++) {
1719: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1720: }
1721: }
1722: if ((rank == size-1) && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1723: for (l=0; l<s; l++) {
1724: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s - l - m + 2)*nc + ofill[k];
1725: }
1726: }
1727: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1728: row++;
1729: }
1730: }
1731: PetscFree2(values,cols);
1732: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1733: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1734: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1735: }
1736: return(0);
1737: }
1739: /* ---------------------------------------------------------------------------------*/
1741: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1742: {
1743: PetscErrorCode ierr;
1744: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1745: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1746: PetscInt istart,iend;
1747: PetscScalar *values;
1748: DMBoundaryType bx;
1749: ISLocalToGlobalMapping ltog,mltog;
1752: /*
1753: nc - number of components per grid point
1754: col - number of colors needed in one direction for single component problem
1756: */
1757: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1758: col = 2*s + 1;
1760: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1761: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1763: MatSetBlockSize(J,nc);
1764: MatSeqAIJSetPreallocation(J,col*nc,0);
1765: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1767: DMGetLocalToGlobalMapping(da,<og);
1768: MatGetLocalToGlobalMapping(J,&mltog,NULL);
1769: if (!mltog) {
1770: MatSetLocalToGlobalMapping(J,ltog,ltog);
1771: }
1773: /*
1774: For each node in the grid: we get the neighbors in the local (on processor ordering
1775: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1776: PETSc ordering.
1777: */
1778: if (!da->prealloc_only) {
1779: PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1780: PetscCalloc1(col*nc*nc,&values);
1781: for (i=xs; i<xs+nx; i++) {
1782: istart = PetscMax(-s,gxs - i);
1783: iend = PetscMin(s,gxs + gnx - i - 1);
1784: slot = i - gxs;
1786: cnt = 0;
1787: for (l=0; l<nc; l++) {
1788: for (i1=istart; i1<iend+1; i1++) {
1789: cols[cnt++] = l + nc*(slot + i1);
1790: }
1791: rows[l] = l + nc*(slot);
1792: }
1793: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1794: }
1795: PetscFree(values);
1796: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1797: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1798: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1799: PetscFree2(rows,cols);
1800: }
1801: return(0);
1802: }
1804: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1805: {
1806: PetscErrorCode ierr;
1807: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1808: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1809: PetscInt istart,iend,jstart,jend,ii,jj;
1810: MPI_Comm comm;
1811: PetscScalar *values;
1812: DMBoundaryType bx,by;
1813: DMDAStencilType st;
1814: ISLocalToGlobalMapping ltog;
1817: /*
1818: nc - number of components per grid point
1819: col - number of colors needed in one direction for single component problem
1820: */
1821: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1822: col = 2*s + 1;
1824: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1825: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1826: PetscObjectGetComm((PetscObject)da,&comm);
1828: PetscMalloc1(col*col*nc*nc,&cols);
1830: DMGetLocalToGlobalMapping(da,<og);
1832: /* determine the matrix preallocation information */
1833: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1834: for (i=xs; i<xs+nx; i++) {
1835: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1836: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1837: for (j=ys; j<ys+ny; j++) {
1838: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1839: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1840: slot = i - gxs + gnx*(j - gys);
1842: /* Find block columns in block row */
1843: cnt = 0;
1844: for (ii=istart; ii<iend+1; ii++) {
1845: for (jj=jstart; jj<jend+1; jj++) {
1846: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1847: cols[cnt++] = slot + ii + gnx*jj;
1848: }
1849: }
1850: }
1851: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1852: }
1853: }
1854: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1855: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1856: MatPreallocateFinalize(dnz,onz);
1858: MatSetLocalToGlobalMapping(J,ltog,ltog);
1860: /*
1861: For each node in the grid: we get the neighbors in the local (on processor ordering
1862: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1863: PETSc ordering.
1864: */
1865: if (!da->prealloc_only) {
1866: PetscCalloc1(col*col*nc*nc,&values);
1867: for (i=xs; i<xs+nx; i++) {
1868: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1869: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1870: for (j=ys; j<ys+ny; j++) {
1871: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1872: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1873: slot = i - gxs + gnx*(j - gys);
1874: cnt = 0;
1875: for (ii=istart; ii<iend+1; ii++) {
1876: for (jj=jstart; jj<jend+1; jj++) {
1877: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1878: cols[cnt++] = slot + ii + gnx*jj;
1879: }
1880: }
1881: }
1882: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1883: }
1884: }
1885: PetscFree(values);
1886: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1887: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1888: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1889: }
1890: PetscFree(cols);
1891: return(0);
1892: }
1894: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1895: {
1896: PetscErrorCode ierr;
1897: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1898: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1899: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1900: MPI_Comm comm;
1901: PetscScalar *values;
1902: DMBoundaryType bx,by,bz;
1903: DMDAStencilType st;
1904: ISLocalToGlobalMapping ltog;
1907: /*
1908: nc - number of components per grid point
1909: col - number of colors needed in one direction for single component problem
1911: */
1912: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1913: col = 2*s + 1;
1915: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1916: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1917: PetscObjectGetComm((PetscObject)da,&comm);
1919: PetscMalloc1(col*col*col,&cols);
1921: DMGetLocalToGlobalMapping(da,<og);
1923: /* determine the matrix preallocation information */
1924: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1925: for (i=xs; i<xs+nx; i++) {
1926: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1927: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1928: for (j=ys; j<ys+ny; j++) {
1929: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1930: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1931: for (k=zs; k<zs+nz; k++) {
1932: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1933: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1935: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1937: /* Find block columns in block row */
1938: cnt = 0;
1939: for (ii=istart; ii<iend+1; ii++) {
1940: for (jj=jstart; jj<jend+1; jj++) {
1941: for (kk=kstart; kk<kend+1; kk++) {
1942: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1943: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1944: }
1945: }
1946: }
1947: }
1948: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1949: }
1950: }
1951: }
1952: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1953: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1954: MatPreallocateFinalize(dnz,onz);
1956: MatSetLocalToGlobalMapping(J,ltog,ltog);
1958: /*
1959: For each node in the grid: we get the neighbors in the local (on processor ordering
1960: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1961: PETSc ordering.
1962: */
1963: if (!da->prealloc_only) {
1964: PetscCalloc1(col*col*col*nc*nc,&values);
1965: for (i=xs; i<xs+nx; i++) {
1966: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1967: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1968: for (j=ys; j<ys+ny; j++) {
1969: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1970: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1971: for (k=zs; k<zs+nz; k++) {
1972: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1973: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1975: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1977: cnt = 0;
1978: for (ii=istart; ii<iend+1; ii++) {
1979: for (jj=jstart; jj<jend+1; jj++) {
1980: for (kk=kstart; kk<kend+1; kk++) {
1981: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1982: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1983: }
1984: }
1985: }
1986: }
1987: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1988: }
1989: }
1990: }
1991: PetscFree(values);
1992: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1993: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1994: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1995: }
1996: PetscFree(cols);
1997: return(0);
1998: }
2000: /*
2001: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
2002: identify in the local ordering with periodic domain.
2003: */
2004: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
2005: {
2007: PetscInt i,n;
2010: ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
2011: ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
2012: for (i=0,n=0; i<*cnt; i++) {
2013: if (col[i] >= *row) col[n++] = col[i];
2014: }
2015: *cnt = n;
2016: return(0);
2017: }
2019: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
2020: {
2021: PetscErrorCode ierr;
2022: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2023: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
2024: PetscInt istart,iend,jstart,jend,ii,jj;
2025: MPI_Comm comm;
2026: PetscScalar *values;
2027: DMBoundaryType bx,by;
2028: DMDAStencilType st;
2029: ISLocalToGlobalMapping ltog;
2032: /*
2033: nc - number of components per grid point
2034: col - number of colors needed in one direction for single component problem
2035: */
2036: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
2037: col = 2*s + 1;
2039: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
2040: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
2041: PetscObjectGetComm((PetscObject)da,&comm);
2043: PetscMalloc1(col*col*nc*nc,&cols);
2045: DMGetLocalToGlobalMapping(da,<og);
2047: /* determine the matrix preallocation information */
2048: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
2049: for (i=xs; i<xs+nx; i++) {
2050: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2051: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2052: for (j=ys; j<ys+ny; j++) {
2053: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2054: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2055: slot = i - gxs + gnx*(j - gys);
2057: /* Find block columns in block row */
2058: cnt = 0;
2059: for (ii=istart; ii<iend+1; ii++) {
2060: for (jj=jstart; jj<jend+1; jj++) {
2061: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2062: cols[cnt++] = slot + ii + gnx*jj;
2063: }
2064: }
2065: }
2066: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2067: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2068: }
2069: }
2070: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2071: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2072: MatPreallocateFinalize(dnz,onz);
2074: MatSetLocalToGlobalMapping(J,ltog,ltog);
2076: /*
2077: For each node in the grid: we get the neighbors in the local (on processor ordering
2078: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2079: PETSc ordering.
2080: */
2081: if (!da->prealloc_only) {
2082: PetscCalloc1(col*col*nc*nc,&values);
2083: for (i=xs; i<xs+nx; i++) {
2084: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2085: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2086: for (j=ys; j<ys+ny; j++) {
2087: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2088: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2089: slot = i - gxs + gnx*(j - gys);
2091: /* Find block columns in block row */
2092: cnt = 0;
2093: for (ii=istart; ii<iend+1; ii++) {
2094: for (jj=jstart; jj<jend+1; jj++) {
2095: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2096: cols[cnt++] = slot + ii + gnx*jj;
2097: }
2098: }
2099: }
2100: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2101: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2102: }
2103: }
2104: PetscFree(values);
2105: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2106: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2107: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2108: }
2109: PetscFree(cols);
2110: return(0);
2111: }
2113: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
2114: {
2115: PetscErrorCode ierr;
2116: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2117: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
2118: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
2119: MPI_Comm comm;
2120: PetscScalar *values;
2121: DMBoundaryType bx,by,bz;
2122: DMDAStencilType st;
2123: ISLocalToGlobalMapping ltog;
2126: /*
2127: nc - number of components per grid point
2128: col - number of colors needed in one direction for single component problem
2129: */
2130: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2131: col = 2*s + 1;
2133: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2134: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2135: PetscObjectGetComm((PetscObject)da,&comm);
2137: /* create the matrix */
2138: PetscMalloc1(col*col*col,&cols);
2140: DMGetLocalToGlobalMapping(da,<og);
2142: /* determine the matrix preallocation information */
2143: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2144: for (i=xs; i<xs+nx; i++) {
2145: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2146: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2147: for (j=ys; j<ys+ny; j++) {
2148: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2149: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2150: for (k=zs; k<zs+nz; k++) {
2151: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2152: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2154: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2156: /* Find block columns in block row */
2157: cnt = 0;
2158: for (ii=istart; ii<iend+1; ii++) {
2159: for (jj=jstart; jj<jend+1; jj++) {
2160: for (kk=kstart; kk<kend+1; kk++) {
2161: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2162: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2163: }
2164: }
2165: }
2166: }
2167: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2168: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2169: }
2170: }
2171: }
2172: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2173: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2174: MatPreallocateFinalize(dnz,onz);
2176: MatSetLocalToGlobalMapping(J,ltog,ltog);
2178: /*
2179: For each node in the grid: we get the neighbors in the local (on processor ordering
2180: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2181: PETSc ordering.
2182: */
2183: if (!da->prealloc_only) {
2184: PetscCalloc1(col*col*col*nc*nc,&values);
2185: for (i=xs; i<xs+nx; i++) {
2186: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2187: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2188: for (j=ys; j<ys+ny; j++) {
2189: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2190: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2191: for (k=zs; k<zs+nz; k++) {
2192: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2193: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2195: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2197: cnt = 0;
2198: for (ii=istart; ii<iend+1; ii++) {
2199: for (jj=jstart; jj<jend+1; jj++) {
2200: for (kk=kstart; kk<kend+1; kk++) {
2201: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2202: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2203: }
2204: }
2205: }
2206: }
2207: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2208: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2209: }
2210: }
2211: }
2212: PetscFree(values);
2213: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2214: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2215: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2216: }
2217: PetscFree(cols);
2218: return(0);
2219: }
2221: /* ---------------------------------------------------------------------------------*/
2223: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2224: {
2225: PetscErrorCode ierr;
2226: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2227: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2228: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2229: DM_DA *dd = (DM_DA*)da->data;
2230: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2231: MPI_Comm comm;
2232: PetscScalar *values;
2233: DMBoundaryType bx,by,bz;
2234: ISLocalToGlobalMapping ltog;
2235: DMDAStencilType st;
2236: PetscBool removedups = PETSC_FALSE;
2239: /*
2240: nc - number of components per grid point
2241: col - number of colors needed in one direction for single component problem
2243: */
2244: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2245: col = 2*s + 1;
2246: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
2247: by 2*stencil_width + 1\n");
2248: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
2249: by 2*stencil_width + 1\n");
2250: if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
2251: by 2*stencil_width + 1\n");
2253: /*
2254: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2255: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2256: */
2257: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2258: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2259: if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;
2261: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2262: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2263: PetscObjectGetComm((PetscObject)da,&comm);
2265: PetscMalloc1(col*col*col*nc,&cols);
2266: DMGetLocalToGlobalMapping(da,<og);
2268: /* determine the matrix preallocation information */
2269: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
2271: MatSetBlockSize(J,nc);
2272: for (i=xs; i<xs+nx; i++) {
2273: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2274: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2275: for (j=ys; j<ys+ny; j++) {
2276: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2277: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2278: for (k=zs; k<zs+nz; k++) {
2279: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2280: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2282: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2284: for (l=0; l<nc; l++) {
2285: cnt = 0;
2286: for (ii=istart; ii<iend+1; ii++) {
2287: for (jj=jstart; jj<jend+1; jj++) {
2288: for (kk=kstart; kk<kend+1; kk++) {
2289: if (ii || jj || kk) {
2290: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2291: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2292: }
2293: } else {
2294: if (dfill) {
2295: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2296: } else {
2297: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2298: }
2299: }
2300: }
2301: }
2302: }
2303: row = l + nc*(slot);
2304: maxcnt = PetscMax(maxcnt,cnt);
2305: if (removedups) {
2306: MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2307: } else {
2308: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2309: }
2310: }
2311: }
2312: }
2313: }
2314: MatSeqAIJSetPreallocation(J,0,dnz);
2315: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2316: MatPreallocateFinalize(dnz,onz);
2317: MatSetLocalToGlobalMapping(J,ltog,ltog);
2319: /*
2320: For each node in the grid: we get the neighbors in the local (on processor ordering
2321: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2322: PETSc ordering.
2323: */
2324: if (!da->prealloc_only) {
2325: PetscCalloc1(maxcnt,&values);
2326: for (i=xs; i<xs+nx; i++) {
2327: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2328: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2329: for (j=ys; j<ys+ny; j++) {
2330: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2331: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2332: for (k=zs; k<zs+nz; k++) {
2333: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2334: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2336: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2338: for (l=0; l<nc; l++) {
2339: cnt = 0;
2340: for (ii=istart; ii<iend+1; ii++) {
2341: for (jj=jstart; jj<jend+1; jj++) {
2342: for (kk=kstart; kk<kend+1; kk++) {
2343: if (ii || jj || kk) {
2344: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2345: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2346: }
2347: } else {
2348: if (dfill) {
2349: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2350: } else {
2351: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2352: }
2353: }
2354: }
2355: }
2356: }
2357: row = l + nc*(slot);
2358: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2359: }
2360: }
2361: }
2362: }
2363: PetscFree(values);
2364: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2365: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2366: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2367: }
2368: PetscFree(cols);
2369: return(0);
2370: }