Actual source code: fdda.c
petsc-3.10.5 2019-03-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: if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1580: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1581: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1583: /*
1584: nc - number of components per grid point
1586: */
1587: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
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: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1605: cnt++;
1606: }
1607: }
1608: for (i=s; i<nx-s; i++) {
1609: for (j=0; j<nc; j++) {
1610: cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1611: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1612: cnt++;
1613: }
1614: }
1615: /* coupling with process to the right */
1616: for (i=nx-s; i<nx; i++) {
1617: for (j=0; j<nc; j++) {
1618: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1619: cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1620: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1621: cnt++;
1622: }
1623: }
1625: MatSeqAIJSetPreallocation(J,0,cols);
1626: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1627: PetscFree2(cols,ocols);
1629: DMGetLocalToGlobalMapping(da,<og);
1630: MatSetLocalToGlobalMapping(J,ltog,ltog);
1632: /*
1633: For each node in the grid: we get the neighbors in the local (on processor ordering
1634: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1635: PETSc ordering.
1636: */
1637: if (!da->prealloc_only) {
1638: PetscCalloc2(maxcnt,&values,maxcnt,&cols);
1640: row = xs*nc;
1641: /* coupling with process to the left */
1642: for (i=xs; i<xs+s; i++) {
1643: for (j=0; j<nc; j++) {
1644: cnt = 0;
1645: if (rank) {
1646: for (l=0; l<s; l++) {
1647: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1648: }
1649: }
1650: if (dfill) {
1651: for (k=dfill[j]; k<dfill[j+1]; k++) {
1652: cols[cnt++] = i*nc + dfill[k];
1653: }
1654: } else {
1655: for (k=0; k<nc; k++) {
1656: cols[cnt++] = i*nc + k;
1657: }
1658: }
1659: for (l=0; l<s; l++) {
1660: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1661: }
1662: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1663: row++;
1664: }
1665: }
1666: for (i=xs+s; i<xs+nx-s; i++) {
1667: for (j=0; j<nc; j++) {
1668: cnt = 0;
1669: for (l=0; l<s; l++) {
1670: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1671: }
1672: if (dfill) {
1673: for (k=dfill[j]; k<dfill[j+1]; k++) {
1674: cols[cnt++] = i*nc + dfill[k];
1675: }
1676: } else {
1677: for (k=0; k<nc; k++) {
1678: cols[cnt++] = i*nc + k;
1679: }
1680: }
1681: for (l=0; l<s; l++) {
1682: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1683: }
1684: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1685: row++;
1686: }
1687: }
1688: /* coupling with process to the right */
1689: for (i=xs+nx-s; i<xs+nx; i++) {
1690: for (j=0; j<nc; j++) {
1691: cnt = 0;
1692: for (l=0; l<s; l++) {
1693: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1694: }
1695: if (dfill) {
1696: for (k=dfill[j]; k<dfill[j+1]; k++) {
1697: cols[cnt++] = i*nc + dfill[k];
1698: }
1699: } else {
1700: for (k=0; k<nc; k++) {
1701: cols[cnt++] = i*nc + k;
1702: }
1703: }
1704: if (rank < size-1) {
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: }
1709: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1710: row++;
1711: }
1712: }
1713: PetscFree2(values,cols);
1714: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1715: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1716: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1717: }
1718: return(0);
1719: }
1721: /* ---------------------------------------------------------------------------------*/
1723: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1724: {
1725: PetscErrorCode ierr;
1726: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1727: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1728: PetscInt istart,iend;
1729: PetscScalar *values;
1730: DMBoundaryType bx;
1731: ISLocalToGlobalMapping ltog,mltog;
1734: /*
1735: nc - number of components per grid point
1736: col - number of colors needed in one direction for single component problem
1738: */
1739: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1740: col = 2*s + 1;
1742: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1743: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1745: MatSetBlockSize(J,nc);
1746: MatSeqAIJSetPreallocation(J,col*nc,0);
1747: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1749: DMGetLocalToGlobalMapping(da,<og);
1750: MatGetLocalToGlobalMapping(J,&mltog,NULL);
1751: if (!mltog) {
1752: MatSetLocalToGlobalMapping(J,ltog,ltog);
1753: }
1755: /*
1756: For each node in the grid: we get the neighbors in the local (on processor ordering
1757: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1758: PETSc ordering.
1759: */
1760: if (!da->prealloc_only) {
1761: PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1762: PetscCalloc1(col*nc*nc,&values);
1763: for (i=xs; i<xs+nx; i++) {
1764: istart = PetscMax(-s,gxs - i);
1765: iend = PetscMin(s,gxs + gnx - i - 1);
1766: slot = i - gxs;
1768: cnt = 0;
1769: for (l=0; l<nc; l++) {
1770: for (i1=istart; i1<iend+1; i1++) {
1771: cols[cnt++] = l + nc*(slot + i1);
1772: }
1773: rows[l] = l + nc*(slot);
1774: }
1775: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1776: }
1777: PetscFree(values);
1778: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1779: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1780: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1781: PetscFree2(rows,cols);
1782: }
1783: return(0);
1784: }
1786: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1787: {
1788: PetscErrorCode ierr;
1789: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1790: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1791: PetscInt istart,iend,jstart,jend,ii,jj;
1792: MPI_Comm comm;
1793: PetscScalar *values;
1794: DMBoundaryType bx,by;
1795: DMDAStencilType st;
1796: ISLocalToGlobalMapping ltog;
1799: /*
1800: nc - number of components per grid point
1801: col - number of colors needed in one direction for single component problem
1802: */
1803: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1804: col = 2*s + 1;
1806: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1807: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1808: PetscObjectGetComm((PetscObject)da,&comm);
1810: PetscMalloc1(col*col*nc*nc,&cols);
1812: DMGetLocalToGlobalMapping(da,<og);
1814: /* determine the matrix preallocation information */
1815: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1816: for (i=xs; i<xs+nx; i++) {
1817: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1818: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1819: for (j=ys; j<ys+ny; j++) {
1820: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1821: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1822: slot = i - gxs + gnx*(j - gys);
1824: /* Find block columns in block row */
1825: cnt = 0;
1826: for (ii=istart; ii<iend+1; ii++) {
1827: for (jj=jstart; jj<jend+1; jj++) {
1828: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1829: cols[cnt++] = slot + ii + gnx*jj;
1830: }
1831: }
1832: }
1833: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1834: }
1835: }
1836: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1837: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1838: MatPreallocateFinalize(dnz,onz);
1840: MatSetLocalToGlobalMapping(J,ltog,ltog);
1842: /*
1843: For each node in the grid: we get the neighbors in the local (on processor ordering
1844: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1845: PETSc ordering.
1846: */
1847: if (!da->prealloc_only) {
1848: PetscCalloc1(col*col*nc*nc,&values);
1849: for (i=xs; i<xs+nx; i++) {
1850: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1851: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1852: for (j=ys; j<ys+ny; j++) {
1853: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1854: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1855: slot = i - gxs + gnx*(j - gys);
1856: cnt = 0;
1857: for (ii=istart; ii<iend+1; ii++) {
1858: for (jj=jstart; jj<jend+1; jj++) {
1859: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1860: cols[cnt++] = slot + ii + gnx*jj;
1861: }
1862: }
1863: }
1864: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1865: }
1866: }
1867: PetscFree(values);
1868: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1869: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1870: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1871: }
1872: PetscFree(cols);
1873: return(0);
1874: }
1876: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1877: {
1878: PetscErrorCode ierr;
1879: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1880: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1881: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1882: MPI_Comm comm;
1883: PetscScalar *values;
1884: DMBoundaryType bx,by,bz;
1885: DMDAStencilType st;
1886: ISLocalToGlobalMapping ltog;
1889: /*
1890: nc - number of components per grid point
1891: col - number of colors needed in one direction for single component problem
1893: */
1894: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1895: col = 2*s + 1;
1897: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1898: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1899: PetscObjectGetComm((PetscObject)da,&comm);
1901: PetscMalloc1(col*col*col,&cols);
1903: DMGetLocalToGlobalMapping(da,<og);
1905: /* determine the matrix preallocation information */
1906: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1907: for (i=xs; i<xs+nx; i++) {
1908: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1909: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1910: for (j=ys; j<ys+ny; j++) {
1911: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1912: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1913: for (k=zs; k<zs+nz; k++) {
1914: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1915: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1917: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1919: /* Find block columns in block row */
1920: cnt = 0;
1921: for (ii=istart; ii<iend+1; ii++) {
1922: for (jj=jstart; jj<jend+1; jj++) {
1923: for (kk=kstart; kk<kend+1; kk++) {
1924: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1925: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1926: }
1927: }
1928: }
1929: }
1930: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1931: }
1932: }
1933: }
1934: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1935: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1936: MatPreallocateFinalize(dnz,onz);
1938: MatSetLocalToGlobalMapping(J,ltog,ltog);
1940: /*
1941: For each node in the grid: we get the neighbors in the local (on processor ordering
1942: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1943: PETSc ordering.
1944: */
1945: if (!da->prealloc_only) {
1946: PetscCalloc1(col*col*col*nc*nc,&values);
1947: for (i=xs; i<xs+nx; i++) {
1948: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1949: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1950: for (j=ys; j<ys+ny; j++) {
1951: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1952: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1953: for (k=zs; k<zs+nz; k++) {
1954: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1955: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1957: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1959: cnt = 0;
1960: for (ii=istart; ii<iend+1; ii++) {
1961: for (jj=jstart; jj<jend+1; jj++) {
1962: for (kk=kstart; kk<kend+1; kk++) {
1963: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1964: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1965: }
1966: }
1967: }
1968: }
1969: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1970: }
1971: }
1972: }
1973: PetscFree(values);
1974: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1975: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1976: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1977: }
1978: PetscFree(cols);
1979: return(0);
1980: }
1982: /*
1983: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1984: identify in the local ordering with periodic domain.
1985: */
1986: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1987: {
1989: PetscInt i,n;
1992: ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1993: ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1994: for (i=0,n=0; i<*cnt; i++) {
1995: if (col[i] >= *row) col[n++] = col[i];
1996: }
1997: *cnt = n;
1998: return(0);
1999: }
2001: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
2002: {
2003: PetscErrorCode ierr;
2004: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2005: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
2006: PetscInt istart,iend,jstart,jend,ii,jj;
2007: MPI_Comm comm;
2008: PetscScalar *values;
2009: DMBoundaryType bx,by;
2010: DMDAStencilType st;
2011: ISLocalToGlobalMapping ltog;
2014: /*
2015: nc - number of components per grid point
2016: col - number of colors needed in one direction for single component problem
2017: */
2018: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
2019: col = 2*s + 1;
2021: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
2022: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
2023: PetscObjectGetComm((PetscObject)da,&comm);
2025: PetscMalloc1(col*col*nc*nc,&cols);
2027: DMGetLocalToGlobalMapping(da,<og);
2029: /* determine the matrix preallocation information */
2030: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
2031: for (i=xs; i<xs+nx; i++) {
2032: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2033: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2034: for (j=ys; j<ys+ny; j++) {
2035: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2036: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2037: slot = i - gxs + gnx*(j - gys);
2039: /* Find block columns in block row */
2040: cnt = 0;
2041: for (ii=istart; ii<iend+1; ii++) {
2042: for (jj=jstart; jj<jend+1; jj++) {
2043: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2044: cols[cnt++] = slot + ii + gnx*jj;
2045: }
2046: }
2047: }
2048: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2049: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2050: }
2051: }
2052: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2053: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2054: MatPreallocateFinalize(dnz,onz);
2056: MatSetLocalToGlobalMapping(J,ltog,ltog);
2058: /*
2059: For each node in the grid: we get the neighbors in the local (on processor ordering
2060: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2061: PETSc ordering.
2062: */
2063: if (!da->prealloc_only) {
2064: PetscCalloc1(col*col*nc*nc,&values);
2065: for (i=xs; i<xs+nx; i++) {
2066: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2067: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2068: for (j=ys; j<ys+ny; j++) {
2069: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2070: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2071: slot = i - gxs + gnx*(j - gys);
2073: /* Find block columns in block row */
2074: cnt = 0;
2075: for (ii=istart; ii<iend+1; ii++) {
2076: for (jj=jstart; jj<jend+1; jj++) {
2077: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
2078: cols[cnt++] = slot + ii + gnx*jj;
2079: }
2080: }
2081: }
2082: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2083: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2084: }
2085: }
2086: PetscFree(values);
2087: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2088: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2089: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2090: }
2091: PetscFree(cols);
2092: return(0);
2093: }
2095: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
2096: {
2097: PetscErrorCode ierr;
2098: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2099: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
2100: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
2101: MPI_Comm comm;
2102: PetscScalar *values;
2103: DMBoundaryType bx,by,bz;
2104: DMDAStencilType st;
2105: ISLocalToGlobalMapping ltog;
2108: /*
2109: nc - number of components per grid point
2110: col - number of colors needed in one direction for single component problem
2111: */
2112: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2113: col = 2*s + 1;
2115: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2116: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2117: PetscObjectGetComm((PetscObject)da,&comm);
2119: /* create the matrix */
2120: PetscMalloc1(col*col*col,&cols);
2122: DMGetLocalToGlobalMapping(da,<og);
2124: /* determine the matrix preallocation information */
2125: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2126: for (i=xs; i<xs+nx; i++) {
2127: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2128: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2129: for (j=ys; j<ys+ny; j++) {
2130: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2131: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2132: for (k=zs; k<zs+nz; k++) {
2133: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2134: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2136: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2138: /* Find block columns in block row */
2139: cnt = 0;
2140: for (ii=istart; ii<iend+1; ii++) {
2141: for (jj=jstart; jj<jend+1; jj++) {
2142: for (kk=kstart; kk<kend+1; kk++) {
2143: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2144: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2145: }
2146: }
2147: }
2148: }
2149: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2150: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2151: }
2152: }
2153: }
2154: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2155: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2156: MatPreallocateFinalize(dnz,onz);
2158: MatSetLocalToGlobalMapping(J,ltog,ltog);
2160: /*
2161: For each node in the grid: we get the neighbors in the local (on processor ordering
2162: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2163: PETSc ordering.
2164: */
2165: if (!da->prealloc_only) {
2166: PetscCalloc1(col*col*col*nc*nc,&values);
2167: for (i=xs; i<xs+nx; i++) {
2168: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2169: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2170: for (j=ys; j<ys+ny; j++) {
2171: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2172: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2173: for (k=zs; k<zs+nz; k++) {
2174: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2175: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2177: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2179: cnt = 0;
2180: for (ii=istart; ii<iend+1; ii++) {
2181: for (jj=jstart; jj<jend+1; jj++) {
2182: for (kk=kstart; kk<kend+1; kk++) {
2183: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2184: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2185: }
2186: }
2187: }
2188: }
2189: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2190: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2191: }
2192: }
2193: }
2194: PetscFree(values);
2195: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2196: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2197: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2198: }
2199: PetscFree(cols);
2200: return(0);
2201: }
2203: /* ---------------------------------------------------------------------------------*/
2205: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2206: {
2207: PetscErrorCode ierr;
2208: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2209: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2210: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2211: DM_DA *dd = (DM_DA*)da->data;
2212: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2213: MPI_Comm comm;
2214: PetscScalar *values;
2215: DMBoundaryType bx,by,bz;
2216: ISLocalToGlobalMapping ltog;
2217: DMDAStencilType st;
2218: PetscBool removedups = PETSC_FALSE;
2221: /*
2222: nc - number of components per grid point
2223: col - number of colors needed in one direction for single component problem
2225: */
2226: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2227: col = 2*s + 1;
2228: 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\
2229: by 2*stencil_width + 1\n");
2230: 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\
2231: by 2*stencil_width + 1\n");
2232: 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\
2233: by 2*stencil_width + 1\n");
2235: /*
2236: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2237: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2238: */
2239: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2240: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2241: if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;
2243: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2244: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2245: PetscObjectGetComm((PetscObject)da,&comm);
2247: PetscMalloc1(col*col*col*nc,&cols);
2248: DMGetLocalToGlobalMapping(da,<og);
2250: /* determine the matrix preallocation information */
2251: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
2253: MatSetBlockSize(J,nc);
2254: for (i=xs; i<xs+nx; i++) {
2255: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2256: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2257: for (j=ys; j<ys+ny; j++) {
2258: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2259: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2260: for (k=zs; k<zs+nz; k++) {
2261: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2262: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2264: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2266: for (l=0; l<nc; l++) {
2267: cnt = 0;
2268: for (ii=istart; ii<iend+1; ii++) {
2269: for (jj=jstart; jj<jend+1; jj++) {
2270: for (kk=kstart; kk<kend+1; kk++) {
2271: if (ii || jj || kk) {
2272: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2273: 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);
2274: }
2275: } else {
2276: if (dfill) {
2277: 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);
2278: } else {
2279: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2280: }
2281: }
2282: }
2283: }
2284: }
2285: row = l + nc*(slot);
2286: maxcnt = PetscMax(maxcnt,cnt);
2287: if (removedups) {
2288: MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2289: } else {
2290: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2291: }
2292: }
2293: }
2294: }
2295: }
2296: MatSeqAIJSetPreallocation(J,0,dnz);
2297: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2298: MatPreallocateFinalize(dnz,onz);
2299: MatSetLocalToGlobalMapping(J,ltog,ltog);
2301: /*
2302: For each node in the grid: we get the neighbors in the local (on processor ordering
2303: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2304: PETSc ordering.
2305: */
2306: if (!da->prealloc_only) {
2307: PetscCalloc1(maxcnt,&values);
2308: for (i=xs; i<xs+nx; i++) {
2309: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2310: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2311: for (j=ys; j<ys+ny; j++) {
2312: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2313: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2314: for (k=zs; k<zs+nz; k++) {
2315: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2316: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2318: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2320: for (l=0; l<nc; l++) {
2321: cnt = 0;
2322: for (ii=istart; ii<iend+1; ii++) {
2323: for (jj=jstart; jj<jend+1; jj++) {
2324: for (kk=kstart; kk<kend+1; kk++) {
2325: if (ii || jj || kk) {
2326: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2327: 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);
2328: }
2329: } else {
2330: if (dfill) {
2331: 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);
2332: } else {
2333: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2334: }
2335: }
2336: }
2337: }
2338: }
2339: row = l + nc*(slot);
2340: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2341: }
2342: }
2343: }
2344: }
2345: PetscFree(values);
2346: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2347: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2348: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2349: }
2350: PetscFree(cols);
2351: return(0);
2352: }