Actual source code: fdda.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
3: #include <petscmat.h>
4: #include <petsc-private/matimpl.h>
6: extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring*);
7: extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring*);
8: extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring*);
9: extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring*);
11: /*
12: For ghost i that may be negative or greater than the upper bound this
13: maps it into the 0:m-1 range using periodicity
14: */
15: #define SetInRange(i,m) ((i < 0) ? m+i : ((i >= m) ? i-m : i))
19: static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill)
20: {
22: PetscInt i,j,nz,*fill;
25: if (!dfill) return(0);
27: /* count number nonzeros */
28: nz = 0;
29: for (i=0; i<w; i++) {
30: for (j=0; j<w; j++) {
31: if (dfill[w*i+j]) nz++;
32: }
33: }
34: PetscMalloc((nz + w + 1)*sizeof(PetscInt),&fill);
35: /* construct modified CSR storage of nonzero structure */
36: /* fill[0 -- w] marks starts of each row of column indices (and end of last row)
37: so fill[1] - fill[0] gives number of nonzeros in first row etc */
38: nz = w + 1;
39: for (i=0; i<w; i++) {
40: fill[i] = nz;
41: for (j=0; j<w; j++) {
42: if (dfill[w*i+j]) {
43: fill[nz] = j;
44: nz++;
45: }
46: }
47: }
48: fill[w] = nz;
50: *rfill = fill;
51: return(0);
52: }
56: /*@
57: DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
58: of the matrix returned by DMCreateMatrix().
60: Logically Collective on DMDA
62: Input Parameter:
63: + da - the distributed array
64: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
65: - ofill - the fill pattern in the off-diagonal blocks
68: Level: developer
70: Notes: This only makes sense when you are doing multicomponent problems but using the
71: MPIAIJ matrix format
73: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
74: representing coupling and 0 entries for missing coupling. For example
75: $ dfill[9] = {1, 0, 0,
76: $ 1, 1, 0,
77: $ 0, 1, 1}
78: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
79: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
80: diagonal block).
82: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
83: can be represented in the dfill, ofill format
85: Contributed by Glenn Hammond
87: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
89: @*/
90: PetscErrorCode DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
91: {
92: DM_DA *dd = (DM_DA*)da->data;
94: PetscInt i,k,cnt = 1;
97: DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
98: DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);
100: /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of
101: columns to the left with any nonzeros in them plus 1 */
102: PetscMalloc(dd->w*sizeof(PetscBool),&dd->ofillcols);
103: PetscMemzero(dd->ofillcols,dd->w*sizeof(PetscBool));
104: for (i=0; i<dd->w; i++) {
105: for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
106: }
107: for (i=0; i<dd->w; i++) {
108: if (dd->ofillcols[i]) {
109: dd->ofillcols[i] = cnt++;
110: }
111: }
112: return(0);
113: }
118: PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,MatType mtype,ISColoring *coloring)
119: {
120: PetscErrorCode ierr;
121: PetscInt dim,m,n,p,nc;
122: DMDABoundaryType bx,by,bz;
123: MPI_Comm comm;
124: PetscMPIInt size;
125: PetscBool isBAIJ;
126: DM_DA *dd = (DM_DA*)da->data;
129: /*
130: m
131: ------------------------------------------------------
132: | |
133: | |
134: | ---------------------- |
135: | | | |
136: n | yn | | |
137: | | | |
138: | .--------------------- |
139: | (xs,ys) xn |
140: | . |
141: | (gxs,gys) |
142: | |
143: -----------------------------------------------------
144: */
146: /*
147: nc - number of components per grid point
148: col - number of colors needed in one direction for single component problem
150: */
151: DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);
153: PetscObjectGetComm((PetscObject)da,&comm);
154: MPI_Comm_size(comm,&size);
155: if (ctype == IS_COLORING_GHOSTED) {
156: if (size == 1) {
157: ctype = IS_COLORING_GLOBAL;
158: } else if (dim > 1) {
159: if ((m==1 && bx == DMDA_BOUNDARY_PERIODIC) || (n==1 && by == DMDA_BOUNDARY_PERIODIC) || (p==1 && bz == DMDA_BOUNDARY_PERIODIC)) {
160: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_GHOSTED cannot be used for periodic boundary condition having both ends of the domain on the same process");
161: }
162: }
163: }
165: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
166: matrices is for the blocks, not the individual matrix elements */
167: PetscStrcmp(mtype,MATBAIJ,&isBAIJ);
168: if (!isBAIJ) {PetscStrcmp(mtype,MATMPIBAIJ,&isBAIJ);}
169: if (!isBAIJ) {PetscStrcmp(mtype,MATSEQBAIJ,&isBAIJ);}
170: if (isBAIJ) {
171: dd->w = 1;
172: dd->xs = dd->xs/nc;
173: dd->xe = dd->xe/nc;
174: dd->Xs = dd->Xs/nc;
175: dd->Xe = dd->Xe/nc;
176: }
178: /*
179: We do not provide a getcoloring function in the DMDA operations because
180: the basic DMDA does not know about matrices. We think of DMDA as being more
181: more low-level then matrices.
182: */
183: if (dim == 1) {
184: DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
185: } else if (dim == 2) {
186: DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
187: } else if (dim == 3) {
188: DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
189: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
190: if (isBAIJ) {
191: dd->w = nc;
192: dd->xs = dd->xs*nc;
193: dd->xe = dd->xe*nc;
194: dd->Xs = dd->Xs*nc;
195: dd->Xe = dd->Xe*nc;
196: }
197: return(0);
198: }
200: /* ---------------------------------------------------------------------------------*/
204: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
205: {
206: PetscErrorCode ierr;
207: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
208: PetscInt ncolors;
209: MPI_Comm comm;
210: DMDABoundaryType bx,by;
211: DMDAStencilType st;
212: ISColoringValue *colors;
213: DM_DA *dd = (DM_DA*)da->data;
216: /*
217: nc - number of components per grid point
218: col - number of colors needed in one direction for single component problem
220: */
221: DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
222: col = 2*s + 1;
223: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
224: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
225: PetscObjectGetComm((PetscObject)da,&comm);
227: /* special case as taught to us by Paul Hovland */
228: if (st == DMDA_STENCIL_STAR && s == 1) {
229: DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
230: } else {
232: if (bx == DMDA_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\
233: by 2*stencil_width + 1 (%d)\n", m, col);
234: if (by == DMDA_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\
235: by 2*stencil_width + 1 (%d)\n", n, col);
236: if (ctype == IS_COLORING_GLOBAL) {
237: if (!dd->localcoloring) {
238: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
239: ii = 0;
240: for (j=ys; j<ys+ny; j++) {
241: for (i=xs; i<xs+nx; i++) {
242: for (k=0; k<nc; k++) {
243: colors[ii++] = k + nc*((i % col) + col*(j % col));
244: }
245: }
246: }
247: ncolors = nc + nc*(col-1 + col*(col-1));
248: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);
249: }
250: *coloring = dd->localcoloring;
251: } else if (ctype == IS_COLORING_GHOSTED) {
252: if (!dd->ghostedcoloring) {
253: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
254: ii = 0;
255: for (j=gys; j<gys+gny; j++) {
256: for (i=gxs; i<gxs+gnx; i++) {
257: for (k=0; k<nc; k++) {
258: /* the complicated stuff is to handle periodic boundaries */
259: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
260: }
261: }
262: }
263: ncolors = nc + nc*(col - 1 + col*(col-1));
264: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);
265: /* PetscIntView(ncolors,(PetscInt*)colors,0); */
267: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
268: }
269: *coloring = dd->ghostedcoloring;
270: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
271: }
272: ISColoringReference(*coloring);
273: return(0);
274: }
276: /* ---------------------------------------------------------------------------------*/
280: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
281: {
282: PetscErrorCode ierr;
283: 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;
284: PetscInt ncolors;
285: MPI_Comm comm;
286: DMDABoundaryType bx,by,bz;
287: DMDAStencilType st;
288: ISColoringValue *colors;
289: DM_DA *dd = (DM_DA*)da->data;
292: /*
293: nc - number of components per grid point
294: col - number of colors needed in one direction for single component problem
296: */
297: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
298: col = 2*s + 1;
299: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
300: by 2*stencil_width + 1\n");
301: if (by == DMDA_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
302: by 2*stencil_width + 1\n");
303: if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
304: by 2*stencil_width + 1\n");
306: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
307: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
308: PetscObjectGetComm((PetscObject)da,&comm);
310: /* create the coloring */
311: if (ctype == IS_COLORING_GLOBAL) {
312: if (!dd->localcoloring) {
313: PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);
314: ii = 0;
315: for (k=zs; k<zs+nz; k++) {
316: for (j=ys; j<ys+ny; j++) {
317: for (i=xs; i<xs+nx; i++) {
318: for (l=0; l<nc; l++) {
319: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
320: }
321: }
322: }
323: }
324: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
325: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,&dd->localcoloring);
326: }
327: *coloring = dd->localcoloring;
328: } else if (ctype == IS_COLORING_GHOSTED) {
329: if (!dd->ghostedcoloring) {
330: PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);
331: ii = 0;
332: for (k=gzs; k<gzs+gnz; k++) {
333: for (j=gys; j<gys+gny; j++) {
334: for (i=gxs; i<gxs+gnx; i++) {
335: for (l=0; l<nc; l++) {
336: /* the complicated stuff is to handle periodic boundaries */
337: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
338: }
339: }
340: }
341: }
342: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
343: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,&dd->ghostedcoloring);
344: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
345: }
346: *coloring = dd->ghostedcoloring;
347: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
348: ISColoringReference(*coloring);
349: return(0);
350: }
352: /* ---------------------------------------------------------------------------------*/
356: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
357: {
358: PetscErrorCode ierr;
359: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
360: PetscInt ncolors;
361: MPI_Comm comm;
362: DMDABoundaryType bx;
363: ISColoringValue *colors;
364: DM_DA *dd = (DM_DA*)da->data;
367: /*
368: nc - number of components per grid point
369: col - number of colors needed in one direction for single component problem
371: */
372: DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
373: col = 2*s + 1;
375: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\
376: by 2*stencil_width + 1 %d\n",(int)m,(int)col);
378: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
379: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
380: PetscObjectGetComm((PetscObject)da,&comm);
382: /* create the coloring */
383: if (ctype == IS_COLORING_GLOBAL) {
384: if (!dd->localcoloring) {
385: PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);
386: if (dd->ofillcols) {
387: PetscInt tc = 0;
388: for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
389: i1 = 0;
390: for (i=xs; i<xs+nx; i++) {
391: for (l=0; l<nc; l++) {
392: if (dd->ofillcols[l] && (i % col)) {
393: colors[i1++] = nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
394: } else {
395: colors[i1++] = l;
396: }
397: }
398: }
399: ncolors = nc + 2*s*tc;
400: } else {
401: i1 = 0;
402: for (i=xs; i<xs+nx; i++) {
403: for (l=0; l<nc; l++) {
404: colors[i1++] = l + nc*(i % col);
405: }
406: }
407: ncolors = nc + nc*(col-1);
408: }
409: ISColoringCreate(comm,ncolors,nc*nx,colors,&dd->localcoloring);
410: }
411: *coloring = dd->localcoloring;
412: } else if (ctype == IS_COLORING_GHOSTED) {
413: if (!dd->ghostedcoloring) {
414: PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);
415: i1 = 0;
416: for (i=gxs; i<gxs+gnx; i++) {
417: for (l=0; l<nc; l++) {
418: /* the complicated stuff is to handle periodic boundaries */
419: colors[i1++] = l + nc*(SetInRange(i,m) % col);
420: }
421: }
422: ncolors = nc + nc*(col-1);
423: ISColoringCreate(comm,ncolors,nc*gnx,colors,&dd->ghostedcoloring);
424: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
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: }
434: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
435: {
436: PetscErrorCode ierr;
437: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
438: PetscInt ncolors;
439: MPI_Comm comm;
440: DMDABoundaryType bx,by;
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,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
451: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
452: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
453: PetscObjectGetComm((PetscObject)da,&comm);
455: if (bx == DMDA_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");
456: if (by == DMDA_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");
458: /* create the coloring */
459: if (ctype == IS_COLORING_GLOBAL) {
460: if (!dd->localcoloring) {
461: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
462: ii = 0;
463: for (j=ys; j<ys+ny; j++) {
464: for (i=xs; i<xs+nx; i++) {
465: for (k=0; k<nc; k++) {
466: colors[ii++] = k + nc*((3*j+i) % 5);
467: }
468: }
469: }
470: ncolors = 5*nc;
471: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);
472: }
473: *coloring = dd->localcoloring;
474: } else if (ctype == IS_COLORING_GHOSTED) {
475: if (!dd->ghostedcoloring) {
476: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
477: ii = 0;
478: for (j=gys; j<gys+gny; j++) {
479: for (i=gxs; i<gxs+gnx; i++) {
480: for (k=0; k<nc; k++) {
481: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
482: }
483: }
484: }
485: ncolors = 5*nc;
486: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);
487: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
488: }
489: *coloring = dd->ghostedcoloring;
490: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
491: return(0);
492: }
494: /* =========================================================================== */
495: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
496: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
497: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
498: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
499: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
500: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
501: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
502: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
503: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
504: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
508: /*@C
509: MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix
511: Logically Collective on Mat
513: Input Parameters:
514: + mat - the matrix
515: - da - the da
517: Level: intermediate
519: @*/
520: PetscErrorCode MatSetupDM(Mat mat,DM da)
521: {
527: PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
528: return(0);
529: }
533: PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer)
534: {
535: DM da;
536: PetscErrorCode ierr;
537: const char *prefix;
538: Mat Anatural;
539: AO ao;
540: PetscInt rstart,rend,*petsc,i;
541: IS is;
542: MPI_Comm comm;
543: PetscViewerFormat format;
546: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
547: PetscViewerGetFormat(viewer,&format);
548: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
550: PetscObjectGetComm((PetscObject)A,&comm);
551: MatGetDM(A, &da);
552: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
554: DMDAGetAO(da,&ao);
555: MatGetOwnershipRange(A,&rstart,&rend);
556: PetscMalloc((rend-rstart)*sizeof(PetscInt),&petsc);
557: for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
558: AOApplicationToPetsc(ao,rend-rstart,petsc);
559: ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);
561: /* call viewer on natural ordering */
562: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
563: ISDestroy(&is);
564: PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
565: PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
566: PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
567: MatView(Anatural,viewer);
568: MatDestroy(&Anatural);
569: return(0);
570: }
574: PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer)
575: {
576: DM da;
578: Mat Anatural,Aapp;
579: AO ao;
580: PetscInt rstart,rend,*app,i;
581: IS is;
582: MPI_Comm comm;
585: PetscObjectGetComm((PetscObject)A,&comm);
586: MatGetDM(A, &da);
587: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
589: /* Load the matrix in natural ordering */
590: MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
591: MatSetType(Anatural,((PetscObject)A)->type_name);
592: MatSetSizes(Anatural,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
593: MatLoad(Anatural,viewer);
595: /* Map natural ordering to application ordering and create IS */
596: DMDAGetAO(da,&ao);
597: MatGetOwnershipRange(Anatural,&rstart,&rend);
598: PetscMalloc((rend-rstart)*sizeof(PetscInt),&app);
599: for (i=rstart; i<rend; i++) app[i-rstart] = i;
600: AOPetscToApplication(ao,rend-rstart,app);
601: ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);
603: /* Do permutation and replace header */
604: MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
605: MatHeaderReplace(A,Aapp);
606: ISDestroy(&is);
607: MatDestroy(&Anatural);
608: return(0);
609: }
613: PetscErrorCode DMCreateMatrix_DA(DM da, MatType mtype,Mat *J)
614: {
616: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
617: Mat A;
618: MPI_Comm comm;
619: MatType Atype;
620: PetscSection section, sectionGlobal;
621: void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL;
622: MatType ttype[256];
623: PetscBool flg;
624: PetscMPIInt size;
625: DM_DA *dd = (DM_DA*)da->data;
628: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
629: MatInitializePackage();
630: #endif
631: if (!mtype) mtype = MATAIJ;
632: PetscStrcpy((char*)ttype,mtype);
633: PetscOptionsBegin(PetscObjectComm((PetscObject)da),((PetscObject)da)->prefix,"DMDA options","Mat");
634: PetscOptionsList("-dm_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)ttype,256,&flg);
635: PetscOptionsEnd();
637: DMGetDefaultSection(da, §ion);
638: if (section) {
639: PetscInt bs = -1;
640: PetscInt localSize;
641: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
643: DMGetDefaultGlobalSection(da, §ionGlobal);
644: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
645: MatCreate(PetscObjectComm((PetscObject)da), J);
646: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
647: MatSetType(*J, mtype);
648: MatSetFromOptions(*J);
649: PetscStrcmp(mtype, MATSHELL, &isShell);
650: PetscStrcmp(mtype, MATBAIJ, &isBlock);
651: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
652: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
653: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
654: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
655: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
656: /* Check for symmetric storage */
657: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
658: if (isSymmetric) {
659: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
660: }
661: if (!isShell) {
662: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
664: if (bs < 0) {
665: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
666: PetscInt pStart, pEnd, p, dof;
668: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
669: for (p = pStart; p < pEnd; ++p) {
670: PetscSectionGetDof(sectionGlobal, p, &dof);
671: if (dof) {
672: bs = dof;
673: break;
674: }
675: }
676: } else {
677: bs = 1;
678: }
679: /* Must have same blocksize on all procs (some might have no points) */
680: bsLocal = bs;
681: MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
682: }
683: PetscMalloc4(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz, localSize/bs, PetscInt, &dnzu, localSize/bs, PetscInt, &onzu);
684: PetscMemzero(dnz, localSize/bs * sizeof(PetscInt));
685: PetscMemzero(onz, localSize/bs * sizeof(PetscInt));
686: PetscMemzero(dnzu, localSize/bs * sizeof(PetscInt));
687: PetscMemzero(onzu, localSize/bs * sizeof(PetscInt));
688: /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
689: PetscFree4(dnz, onz, dnzu, onzu);
690: }
691: }
692: /*
693: m
694: ------------------------------------------------------
695: | |
696: | |
697: | ---------------------- |
698: | | | |
699: n | ny | | |
700: | | | |
701: | .--------------------- |
702: | (xs,ys) nx |
703: | . |
704: | (gxs,gys) |
705: | |
706: -----------------------------------------------------
707: */
709: /*
710: nc - number of components per grid point
711: col - number of colors needed in one direction for single component problem
713: */
714: M = dd->M;
715: N = dd->N;
716: P = dd->P;
717: dim = dd->dim;
718: dof = dd->w;
719: /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
720: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
721: PetscObjectGetComm((PetscObject)da,&comm);
722: MatCreate(comm,&A);
723: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
724: MatSetType(A,(MatType)ttype);
725: MatSetDM(A,da);
726: MatSetFromOptions(A);
727: MatGetType(A,&Atype);
728: /*
729: We do not provide a getmatrix function in the DMDA operations because
730: the basic DMDA does not know about matrices. We think of DMDA as being more
731: more low-level than matrices. This is kind of cheating but, cause sometimes
732: we think of DMDA has higher level than matrices.
734: We could switch based on Atype (or mtype), but we do not since the
735: specialized setting routines depend only the particular preallocation
736: details of the matrix, not the type itself.
737: */
738: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
739: if (!aij) {
740: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
741: }
742: if (!aij) {
743: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
744: if (!baij) {
745: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
746: }
747: if (!baij) {
748: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
749: if (!sbaij) {
750: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
751: }
752: }
753: }
754: if (aij) {
755: if (dim == 1) {
756: if (dd->ofill) {
757: DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
758: } else {
759: DMCreateMatrix_DA_1d_MPIAIJ(da,A);
760: }
761: } else if (dim == 2) {
762: if (dd->ofill) {
763: DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
764: } else {
765: DMCreateMatrix_DA_2d_MPIAIJ(da,A);
766: }
767: } else if (dim == 3) {
768: if (dd->ofill) {
769: DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
770: } else {
771: DMCreateMatrix_DA_3d_MPIAIJ(da,A);
772: }
773: }
774: } else if (baij) {
775: if (dim == 2) {
776: DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
777: } else if (dim == 3) {
778: DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
779: } 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);
780: } else if (sbaij) {
781: if (dim == 2) {
782: DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
783: } else if (dim == 3) {
784: DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
785: } 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);
786: } else {
787: ISLocalToGlobalMapping ltog,ltogb;
788: DMGetLocalToGlobalMapping(da,<og);
789: DMGetLocalToGlobalMappingBlock(da,<ogb);
790: MatSetUp(A);
791: MatSetLocalToGlobalMapping(A,ltog,ltog);
792: MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);
793: }
794: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
795: MatSetStencil(A,dim,dims,starts,dof);
796: MatSetDM(A,da);
797: MPI_Comm_size(comm,&size);
798: if (size > 1) {
799: /* change viewer to display matrix in natural ordering */
800: MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
801: MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
802: }
803: *J = A;
804: return(0);
805: }
807: /* ---------------------------------------------------------------------------------*/
810: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
811: {
812: PetscErrorCode ierr;
813: 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;
814: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
815: MPI_Comm comm;
816: PetscScalar *values;
817: DMDABoundaryType bx,by;
818: ISLocalToGlobalMapping ltog,ltogb;
819: DMDAStencilType st;
822: /*
823: nc - number of components per grid point
824: col - number of colors needed in one direction for single component problem
826: */
827: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
828: col = 2*s + 1;
829: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
830: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
831: PetscObjectGetComm((PetscObject)da,&comm);
833: PetscMalloc2(nc,PetscInt,&rows,col*col*nc*nc,PetscInt,&cols);
834: DMGetLocalToGlobalMapping(da,<og);
835: DMGetLocalToGlobalMappingBlock(da,<ogb);
837: /* determine the matrix preallocation information */
838: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
839: for (i=xs; i<xs+nx; i++) {
841: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
842: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
844: for (j=ys; j<ys+ny; j++) {
845: slot = i - gxs + gnx*(j - gys);
847: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
848: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
850: cnt = 0;
851: for (k=0; k<nc; k++) {
852: for (l=lstart; l<lend+1; l++) {
853: for (p=pstart; p<pend+1; p++) {
854: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
855: cols[cnt++] = k + nc*(slot + gnx*l + p);
856: }
857: }
858: }
859: rows[k] = k + nc*(slot);
860: }
861: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
862: }
863: }
864: MatSetBlockSize(J,nc);
865: MatSeqAIJSetPreallocation(J,0,dnz);
866: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
867: MatPreallocateFinalize(dnz,onz);
869: MatSetLocalToGlobalMapping(J,ltog,ltog);
870: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
872: /*
873: For each node in the grid: we get the neighbors in the local (on processor ordering
874: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
875: PETSc ordering.
876: */
877: if (!da->prealloc_only) {
878: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
879: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
880: for (i=xs; i<xs+nx; i++) {
882: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
883: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
885: for (j=ys; j<ys+ny; j++) {
886: slot = i - gxs + gnx*(j - gys);
888: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
889: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
891: cnt = 0;
892: for (k=0; k<nc; k++) {
893: for (l=lstart; l<lend+1; l++) {
894: for (p=pstart; p<pend+1; p++) {
895: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
896: cols[cnt++] = k + nc*(slot + gnx*l + p);
897: }
898: }
899: }
900: rows[k] = k + nc*(slot);
901: }
902: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
903: }
904: }
905: PetscFree(values);
906: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
907: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
908: }
909: PetscFree2(rows,cols);
910: return(0);
911: }
915: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
916: {
917: PetscErrorCode ierr;
918: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
919: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p;
920: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
921: DM_DA *dd = (DM_DA*)da->data;
922: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
923: MPI_Comm comm;
924: PetscScalar *values;
925: DMDABoundaryType bx,by;
926: ISLocalToGlobalMapping ltog,ltogb;
927: DMDAStencilType st;
930: /*
931: nc - number of components per grid point
932: col - number of colors needed in one direction for single component problem
934: */
935: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
936: col = 2*s + 1;
937: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
938: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
939: PetscObjectGetComm((PetscObject)da,&comm);
941: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
942: DMGetLocalToGlobalMapping(da,<og);
943: DMGetLocalToGlobalMappingBlock(da,<ogb);
945: /* determine the matrix preallocation information */
946: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
947: for (i=xs; i<xs+nx; i++) {
949: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
950: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
952: for (j=ys; j<ys+ny; j++) {
953: slot = i - gxs + gnx*(j - gys);
955: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
956: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
958: for (k=0; k<nc; k++) {
959: cnt = 0;
960: for (l=lstart; l<lend+1; l++) {
961: for (p=pstart; p<pend+1; p++) {
962: if (l || p) {
963: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
964: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
965: }
966: } else {
967: if (dfill) {
968: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
969: } else {
970: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
971: }
972: }
973: }
974: }
975: row = k + nc*(slot);
976: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
977: }
978: }
979: }
980: MatSeqAIJSetPreallocation(J,0,dnz);
981: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
982: MatPreallocateFinalize(dnz,onz);
983: MatSetLocalToGlobalMapping(J,ltog,ltog);
984: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
986: /*
987: For each node in the grid: we get the neighbors in the local (on processor ordering
988: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
989: PETSc ordering.
990: */
991: if (!da->prealloc_only) {
992: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
993: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
994: for (i=xs; i<xs+nx; i++) {
996: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
997: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
999: for (j=ys; j<ys+ny; j++) {
1000: slot = i - gxs + gnx*(j - gys);
1002: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1003: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1005: for (k=0; k<nc; k++) {
1006: cnt = 0;
1007: for (l=lstart; l<lend+1; l++) {
1008: for (p=pstart; p<pend+1; p++) {
1009: if (l || p) {
1010: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1011: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1012: }
1013: } else {
1014: if (dfill) {
1015: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1016: } else {
1017: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1018: }
1019: }
1020: }
1021: }
1022: row = k + nc*(slot);
1023: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1024: }
1025: }
1026: }
1027: PetscFree(values);
1028: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1029: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1030: }
1031: PetscFree(cols);
1032: return(0);
1033: }
1035: /* ---------------------------------------------------------------------------------*/
1039: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1040: {
1041: PetscErrorCode ierr;
1042: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1043: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1044: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1045: MPI_Comm comm;
1046: PetscScalar *values;
1047: DMDABoundaryType bx,by,bz;
1048: ISLocalToGlobalMapping ltog,ltogb;
1049: DMDAStencilType st;
1052: /*
1053: nc - number of components per grid point
1054: col - number of colors needed in one direction for single component problem
1056: */
1057: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1058: col = 2*s + 1;
1060: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1061: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1062: PetscObjectGetComm((PetscObject)da,&comm);
1064: PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);
1065: DMGetLocalToGlobalMapping(da,<og);
1066: DMGetLocalToGlobalMappingBlock(da,<ogb);
1068: /* determine the matrix preallocation information */
1069: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1070: for (i=xs; i<xs+nx; i++) {
1071: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1072: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1073: for (j=ys; j<ys+ny; j++) {
1074: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1075: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1076: for (k=zs; k<zs+nz; k++) {
1077: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1078: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1080: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1082: cnt = 0;
1083: for (l=0; l<nc; l++) {
1084: for (ii=istart; ii<iend+1; ii++) {
1085: for (jj=jstart; jj<jend+1; jj++) {
1086: for (kk=kstart; kk<kend+1; kk++) {
1087: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1088: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1089: }
1090: }
1091: }
1092: }
1093: rows[l] = l + nc*(slot);
1094: }
1095: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1096: }
1097: }
1098: }
1099: MatSetBlockSize(J,nc);
1100: MatSeqAIJSetPreallocation(J,0,dnz);
1101: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1102: MatPreallocateFinalize(dnz,onz);
1103: MatSetLocalToGlobalMapping(J,ltog,ltog);
1104: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1106: /*
1107: For each node in the grid: we get the neighbors in the local (on processor ordering
1108: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1109: PETSc ordering.
1110: */
1111: if (!da->prealloc_only) {
1112: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1113: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1114: for (i=xs; i<xs+nx; i++) {
1115: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1116: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1117: for (j=ys; j<ys+ny; j++) {
1118: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1119: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1120: for (k=zs; k<zs+nz; k++) {
1121: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1122: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1124: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1126: cnt = 0;
1127: for (l=0; l<nc; l++) {
1128: for (ii=istart; ii<iend+1; ii++) {
1129: for (jj=jstart; jj<jend+1; jj++) {
1130: for (kk=kstart; kk<kend+1; kk++) {
1131: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1132: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1133: }
1134: }
1135: }
1136: }
1137: rows[l] = l + nc*(slot);
1138: }
1139: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1140: }
1141: }
1142: }
1143: PetscFree(values);
1144: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1145: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1146: }
1147: PetscFree2(rows,cols);
1148: return(0);
1149: }
1151: /* ---------------------------------------------------------------------------------*/
1155: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1156: {
1157: PetscErrorCode ierr;
1158: DM_DA *dd = (DM_DA*)da->data;
1159: PetscInt xs,nx,i,j,gxs,gnx,row,k,l;
1160: PetscInt m,dim,s,*cols = NULL,nc,col,cnt,*ocols;
1161: PetscInt *ofill = dd->ofill;
1162: PetscScalar *values;
1163: DMDABoundaryType bx;
1164: ISLocalToGlobalMapping ltog,ltogb;
1165: PetscMPIInt rank,size;
1168: if (dd->bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1169: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1170: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1172: /*
1173: nc - number of components per grid point
1174: col - number of colors needed in one direction for single component problem
1176: */
1177: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1178: col = 2*s + 1;
1180: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1181: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1183: MatSetBlockSize(J,nc);
1184: PetscMalloc2(nx*nc,PetscInt,&cols,nx*nc,PetscInt,&ocols);
1185: PetscMemzero(cols,nx*nc*sizeof(PetscInt));
1186: PetscMemzero(ocols,nx*nc*sizeof(PetscInt));
1188: /*
1189: note should be smaller for first and last process with no periodic
1190: does not handle dfill
1191: */
1192: cnt = 0;
1193: /* coupling with process to the left */
1194: for (i=0; i<s; i++) {
1195: for (j=0; j<nc; j++) {
1196: ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1197: cols[cnt] = nc + (s + i)*(ofill[j+1] - ofill[j]);
1198: cnt++;
1199: }
1200: }
1201: for (i=s; i<nx-s; i++) {
1202: for (j=0; j<nc; j++) {
1203: cols[cnt] = nc + 2*s*(ofill[j+1] - ofill[j]);
1204: cnt++;
1205: }
1206: }
1207: /* coupling with process to the right */
1208: for (i=nx-s; i<nx; i++) {
1209: for (j=0; j<nc; j++) {
1210: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1211: cols[cnt] = nc + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1212: cnt++;
1213: }
1214: }
1216: MatSeqAIJSetPreallocation(J,0,cols);
1217: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1218: PetscFree2(cols,ocols);
1220: DMGetLocalToGlobalMapping(da,<og);
1221: DMGetLocalToGlobalMappingBlock(da,<ogb);
1222: MatSetLocalToGlobalMapping(J,ltog,ltog);
1223: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1225: /*
1226: For each node in the grid: we get the neighbors in the local (on processor ordering
1227: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1228: PETSc ordering.
1229: */
1230: if (!da->prealloc_only) {
1231: PetscMalloc(col*nc*nc*sizeof(PetscInt),&cols);
1232: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
1233: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
1235: row = xs*nc;
1236: /* coupling with process to the left */
1237: for (i=xs; i<xs+s; i++) {
1238: for (j=0; j<nc; j++) {
1239: cnt = 0;
1240: if (rank) {
1241: for (l=0; l<s; l++) {
1242: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1243: }
1244: }
1245: for (k=0; k<nc; k++) {
1246: cols[cnt++] = i*nc + k;
1247: }
1248: for (l=0; l<s; l++) {
1249: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1250: }
1251: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1252: row++;
1253: }
1254: }
1255: for (i=xs+s; i<xs+nx-s; i++) {
1256: for (j=0; j<nc; j++) {
1257: cnt = 0;
1258: for (l=0; l<s; l++) {
1259: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1260: }
1261: for (k=0; k<nc; k++) {
1262: cols[cnt++] = i*nc + k;
1263: }
1264: for (l=0; l<s; l++) {
1265: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1266: }
1267: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1268: row++;
1269: }
1270: }
1271: /* coupling with process to the right */
1272: for (i=xs+nx-s; i<xs+nx; i++) {
1273: for (j=0; j<nc; j++) {
1274: cnt = 0;
1275: for (l=0; l<s; l++) {
1276: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1277: }
1278: for (k=0; k<nc; k++) {
1279: cols[cnt++] = i*nc + k;
1280: }
1281: if (rank < size-1) {
1282: for (l=0; l<s; l++) {
1283: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1284: }
1285: }
1286: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1287: row++;
1288: }
1289: }
1290: PetscFree(values);
1291: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1292: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1293: PetscFree(cols);
1294: }
1295: return(0);
1296: }
1298: /* ---------------------------------------------------------------------------------*/
1302: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1303: {
1304: PetscErrorCode ierr;
1305: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1306: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1307: PetscInt istart,iend;
1308: PetscScalar *values;
1309: DMDABoundaryType bx;
1310: ISLocalToGlobalMapping ltog,ltogb;
1313: /*
1314: nc - number of components per grid point
1315: col - number of colors needed in one direction for single component problem
1317: */
1318: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1319: col = 2*s + 1;
1321: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1322: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1324: MatSetBlockSize(J,nc);
1325: MatSeqAIJSetPreallocation(J,col*nc,0);
1326: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1328: DMGetLocalToGlobalMapping(da,<og);
1329: DMGetLocalToGlobalMappingBlock(da,<ogb);
1330: MatSetLocalToGlobalMapping(J,ltog,ltog);
1331: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1333: /*
1334: For each node in the grid: we get the neighbors in the local (on processor ordering
1335: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1336: PETSc ordering.
1337: */
1338: if (!da->prealloc_only) {
1339: PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);
1340: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
1341: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
1342: for (i=xs; i<xs+nx; i++) {
1343: istart = PetscMax(-s,gxs - i);
1344: iend = PetscMin(s,gxs + gnx - i - 1);
1345: slot = i - gxs;
1347: cnt = 0;
1348: for (l=0; l<nc; l++) {
1349: for (i1=istart; i1<iend+1; i1++) {
1350: cols[cnt++] = l + nc*(slot + i1);
1351: }
1352: rows[l] = l + nc*(slot);
1353: }
1354: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1355: }
1356: PetscFree(values);
1357: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1358: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1359: PetscFree2(rows,cols);
1360: }
1361: return(0);
1362: }
1366: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1367: {
1368: PetscErrorCode ierr;
1369: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1370: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1371: PetscInt istart,iend,jstart,jend,ii,jj;
1372: MPI_Comm comm;
1373: PetscScalar *values;
1374: DMDABoundaryType bx,by;
1375: DMDAStencilType st;
1376: ISLocalToGlobalMapping ltog,ltogb;
1379: /*
1380: nc - number of components per grid point
1381: col - number of colors needed in one direction for single component problem
1382: */
1383: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1384: col = 2*s + 1;
1386: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1387: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1388: PetscObjectGetComm((PetscObject)da,&comm);
1390: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1392: DMGetLocalToGlobalMapping(da,<og);
1393: DMGetLocalToGlobalMappingBlock(da,<ogb);
1395: /* determine the matrix preallocation information */
1396: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1397: for (i=xs; i<xs+nx; i++) {
1398: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1399: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1400: for (j=ys; j<ys+ny; j++) {
1401: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1402: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1403: slot = i - gxs + gnx*(j - gys);
1405: /* Find block columns in block row */
1406: cnt = 0;
1407: for (ii=istart; ii<iend+1; ii++) {
1408: for (jj=jstart; jj<jend+1; jj++) {
1409: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1410: cols[cnt++] = slot + ii + gnx*jj;
1411: }
1412: }
1413: }
1414: MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);
1415: }
1416: }
1417: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1418: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1419: MatPreallocateFinalize(dnz,onz);
1421: MatSetLocalToGlobalMapping(J,ltog,ltog);
1422: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1424: /*
1425: For each node in the grid: we get the neighbors in the local (on processor ordering
1426: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1427: PETSc ordering.
1428: */
1429: if (!da->prealloc_only) {
1430: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1431: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1432: for (i=xs; i<xs+nx; i++) {
1433: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1434: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1435: for (j=ys; j<ys+ny; j++) {
1436: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1437: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1438: slot = i - gxs + gnx*(j - gys);
1439: cnt = 0;
1440: for (ii=istart; ii<iend+1; ii++) {
1441: for (jj=jstart; jj<jend+1; jj++) {
1442: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1443: cols[cnt++] = slot + ii + gnx*jj;
1444: }
1445: }
1446: }
1447: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1448: }
1449: }
1450: PetscFree(values);
1451: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1452: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1453: }
1454: PetscFree(cols);
1455: return(0);
1456: }
1460: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1461: {
1462: PetscErrorCode ierr;
1463: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1464: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1465: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1466: MPI_Comm comm;
1467: PetscScalar *values;
1468: DMDABoundaryType bx,by,bz;
1469: DMDAStencilType st;
1470: ISLocalToGlobalMapping ltog,ltogb;
1473: /*
1474: nc - number of components per grid point
1475: col - number of colors needed in one direction for single component problem
1477: */
1478: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1479: col = 2*s + 1;
1481: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1482: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1483: PetscObjectGetComm((PetscObject)da,&comm);
1485: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1487: DMGetLocalToGlobalMapping(da,<og);
1488: DMGetLocalToGlobalMappingBlock(da,<ogb);
1490: /* determine the matrix preallocation information */
1491: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1492: for (i=xs; i<xs+nx; i++) {
1493: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1494: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1495: for (j=ys; j<ys+ny; j++) {
1496: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1497: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1498: for (k=zs; k<zs+nz; k++) {
1499: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1500: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1502: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1504: /* Find block columns in block row */
1505: cnt = 0;
1506: for (ii=istart; ii<iend+1; ii++) {
1507: for (jj=jstart; jj<jend+1; jj++) {
1508: for (kk=kstart; kk<kend+1; kk++) {
1509: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1510: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1511: }
1512: }
1513: }
1514: }
1515: MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);
1516: }
1517: }
1518: }
1519: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1520: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1521: MatPreallocateFinalize(dnz,onz);
1523: MatSetLocalToGlobalMapping(J,ltog,ltog);
1524: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1526: /*
1527: For each node in the grid: we get the neighbors in the local (on processor ordering
1528: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1529: PETSc ordering.
1530: */
1531: if (!da->prealloc_only) {
1532: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1533: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1534: for (i=xs; i<xs+nx; i++) {
1535: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1536: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1537: for (j=ys; j<ys+ny; j++) {
1538: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1539: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1540: for (k=zs; k<zs+nz; k++) {
1541: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1542: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1544: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1546: cnt = 0;
1547: for (ii=istart; ii<iend+1; ii++) {
1548: for (jj=jstart; jj<jend+1; jj++) {
1549: for (kk=kstart; kk<kend+1; kk++) {
1550: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1551: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1552: }
1553: }
1554: }
1555: }
1556: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1557: }
1558: }
1559: }
1560: PetscFree(values);
1561: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1562: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1563: }
1564: PetscFree(cols);
1565: return(0);
1566: }
1570: /*
1571: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1572: identify in the local ordering with periodic domain.
1573: */
1574: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1575: {
1577: PetscInt i,n;
1580: ISLocalToGlobalMappingApply(ltog,1,row,row);
1581: ISLocalToGlobalMappingApply(ltog,*cnt,col,col);
1582: for (i=0,n=0; i<*cnt; i++) {
1583: if (col[i] >= *row) col[n++] = col[i];
1584: }
1585: *cnt = n;
1586: return(0);
1587: }
1591: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1592: {
1593: PetscErrorCode ierr;
1594: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1595: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1596: PetscInt istart,iend,jstart,jend,ii,jj;
1597: MPI_Comm comm;
1598: PetscScalar *values;
1599: DMDABoundaryType bx,by;
1600: DMDAStencilType st;
1601: ISLocalToGlobalMapping ltog,ltogb;
1604: /*
1605: nc - number of components per grid point
1606: col - number of colors needed in one direction for single component problem
1607: */
1608: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1609: col = 2*s + 1;
1611: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1612: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1613: PetscObjectGetComm((PetscObject)da,&comm);
1615: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1617: DMGetLocalToGlobalMapping(da,<og);
1618: DMGetLocalToGlobalMappingBlock(da,<ogb);
1620: /* determine the matrix preallocation information */
1621: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1622: for (i=xs; i<xs+nx; i++) {
1623: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1624: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1625: for (j=ys; j<ys+ny; j++) {
1626: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1627: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1628: slot = i - gxs + gnx*(j - gys);
1630: /* Find block columns in block row */
1631: cnt = 0;
1632: for (ii=istart; ii<iend+1; ii++) {
1633: for (jj=jstart; jj<jend+1; jj++) {
1634: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1635: cols[cnt++] = slot + ii + gnx*jj;
1636: }
1637: }
1638: }
1639: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1640: MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);
1641: }
1642: }
1643: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1644: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1645: MatPreallocateFinalize(dnz,onz);
1647: MatSetLocalToGlobalMapping(J,ltog,ltog);
1648: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1650: /*
1651: For each node in the grid: we get the neighbors in the local (on processor ordering
1652: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1653: PETSc ordering.
1654: */
1655: if (!da->prealloc_only) {
1656: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1657: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1658: for (i=xs; i<xs+nx; i++) {
1659: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1660: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1661: for (j=ys; j<ys+ny; j++) {
1662: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1663: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1664: slot = i - gxs + gnx*(j - gys);
1666: /* Find block columns in block row */
1667: cnt = 0;
1668: for (ii=istart; ii<iend+1; ii++) {
1669: for (jj=jstart; jj<jend+1; jj++) {
1670: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1671: cols[cnt++] = slot + ii + gnx*jj;
1672: }
1673: }
1674: }
1675: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1676: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1677: }
1678: }
1679: PetscFree(values);
1680: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1681: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1682: }
1683: PetscFree(cols);
1684: return(0);
1685: }
1689: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1690: {
1691: PetscErrorCode ierr;
1692: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1693: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1694: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1695: MPI_Comm comm;
1696: PetscScalar *values;
1697: DMDABoundaryType bx,by,bz;
1698: DMDAStencilType st;
1699: ISLocalToGlobalMapping ltog,ltogb;
1702: /*
1703: nc - number of components per grid point
1704: col - number of colors needed in one direction for single component problem
1705: */
1706: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1707: col = 2*s + 1;
1709: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1710: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1711: PetscObjectGetComm((PetscObject)da,&comm);
1713: /* create the matrix */
1714: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1716: DMGetLocalToGlobalMapping(da,<og);
1717: DMGetLocalToGlobalMappingBlock(da,<ogb);
1719: /* determine the matrix preallocation information */
1720: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1721: for (i=xs; i<xs+nx; i++) {
1722: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1723: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1724: for (j=ys; j<ys+ny; j++) {
1725: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1726: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1727: for (k=zs; k<zs+nz; k++) {
1728: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1729: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1731: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1733: /* Find block columns in block row */
1734: cnt = 0;
1735: for (ii=istart; ii<iend+1; ii++) {
1736: for (jj=jstart; jj<jend+1; jj++) {
1737: for (kk=kstart; kk<kend+1; kk++) {
1738: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1739: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1740: }
1741: }
1742: }
1743: }
1744: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1745: MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);
1746: }
1747: }
1748: }
1749: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1750: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1751: MatPreallocateFinalize(dnz,onz);
1753: MatSetLocalToGlobalMapping(J,ltog,ltog);
1754: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1756: /*
1757: For each node in the grid: we get the neighbors in the local (on processor ordering
1758: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1759: PETSc ordering.
1760: */
1761: if (!da->prealloc_only) {
1762: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1763: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1764: for (i=xs; i<xs+nx; i++) {
1765: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1766: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1767: for (j=ys; j<ys+ny; j++) {
1768: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1769: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1770: for (k=zs; k<zs+nz; k++) {
1771: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1772: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1774: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1776: cnt = 0;
1777: for (ii=istart; ii<iend+1; ii++) {
1778: for (jj=jstart; jj<jend+1; jj++) {
1779: for (kk=kstart; kk<kend+1; kk++) {
1780: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1781: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1782: }
1783: }
1784: }
1785: }
1786: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1787: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1788: }
1789: }
1790: }
1791: PetscFree(values);
1792: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1793: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1794: }
1795: PetscFree(cols);
1796: return(0);
1797: }
1799: /* ---------------------------------------------------------------------------------*/
1803: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
1804: {
1805: PetscErrorCode ierr;
1806: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1807: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz;
1808: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1809: DM_DA *dd = (DM_DA*)da->data;
1810: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
1811: MPI_Comm comm;
1812: PetscScalar *values;
1813: DMDABoundaryType bx,by,bz;
1814: ISLocalToGlobalMapping ltog,ltogb;
1815: DMDAStencilType st;
1818: /*
1819: nc - number of components per grid point
1820: col - number of colors needed in one direction for single component problem
1822: */
1823: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1824: col = 2*s + 1;
1825: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1826: by 2*stencil_width + 1\n");
1827: if (by == DMDA_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1828: by 2*stencil_width + 1\n");
1829: if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1830: by 2*stencil_width + 1\n");
1832: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1833: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1834: PetscObjectGetComm((PetscObject)da,&comm);
1836: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
1837: DMGetLocalToGlobalMapping(da,<og);
1838: DMGetLocalToGlobalMappingBlock(da,<ogb);
1840: /* determine the matrix preallocation information */
1841: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1844: for (i=xs; i<xs+nx; i++) {
1845: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1846: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1847: for (j=ys; j<ys+ny; j++) {
1848: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1849: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1850: for (k=zs; k<zs+nz; k++) {
1851: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1852: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1854: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1856: for (l=0; l<nc; l++) {
1857: cnt = 0;
1858: for (ii=istart; ii<iend+1; ii++) {
1859: for (jj=jstart; jj<jend+1; jj++) {
1860: for (kk=kstart; kk<kend+1; kk++) {
1861: if (ii || jj || kk) {
1862: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1863: 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);
1864: }
1865: } else {
1866: if (dfill) {
1867: 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);
1868: } else {
1869: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1870: }
1871: }
1872: }
1873: }
1874: }
1875: row = l + nc*(slot);
1876: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1877: }
1878: }
1879: }
1880: }
1881: MatSeqAIJSetPreallocation(J,0,dnz);
1882: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1883: MatPreallocateFinalize(dnz,onz);
1884: MatSetLocalToGlobalMapping(J,ltog,ltog);
1885: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1887: /*
1888: For each node in the grid: we get the neighbors in the local (on processor ordering
1889: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1890: PETSc ordering.
1891: */
1892: if (!da->prealloc_only) {
1893: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1894: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1895: for (i=xs; i<xs+nx; i++) {
1896: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1897: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1898: for (j=ys; j<ys+ny; j++) {
1899: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1900: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1901: for (k=zs; k<zs+nz; k++) {
1902: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1903: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1905: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1907: for (l=0; l<nc; l++) {
1908: cnt = 0;
1909: for (ii=istart; ii<iend+1; ii++) {
1910: for (jj=jstart; jj<jend+1; jj++) {
1911: for (kk=kstart; kk<kend+1; kk++) {
1912: if (ii || jj || kk) {
1913: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1914: 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);
1915: }
1916: } else {
1917: if (dfill) {
1918: 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);
1919: } else {
1920: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1921: }
1922: }
1923: }
1924: }
1925: }
1926: row = l + nc*(slot);
1927: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1928: }
1929: }
1930: }
1931: }
1932: PetscFree(values);
1933: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1934: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1935: }
1936: PetscFree(cols);
1937: return(0);
1938: }