Actual source code: fdda.c
petsc-3.3-p7 2013-05-11
1:
2: #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/
3: #include <petscmat.h> /*I "petscmat.h" I*/
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(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: nz = w + 1;
37: for (i=0; i<w; i++) {
38: fill[i] = nz;
39: for (j=0; j<w; j++) {
40: if (dfill[w*i+j]) {
41: fill[nz] = j;
42: nz++;
43: }
44: }
45: }
46: fill[w] = nz;
47:
48: *rfill = fill;
49: return(0);
50: }
54: /*@
55: DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
56: of the matrix returned by DMCreateMatrix().
58: Logically Collective on DMDA
60: Input Parameter:
61: + da - the distributed array
62: . dfill - the fill pattern in the diagonal block (may be PETSC_NULL, means use dense block)
63: - ofill - the fill pattern in the off-diagonal blocks
66: Level: developer
68: Notes: This only makes sense when you are doing multicomponent problems but using the
69: MPIAIJ matrix format
71: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
72: representing coupling and 0 entries for missing coupling. For example
73: $ dfill[9] = {1, 0, 0,
74: $ 1, 1, 0,
75: $ 0, 1, 1}
76: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
77: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
78: diagonal block).
80: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
81: can be represented in the dfill, ofill format
83: Contributed by Glenn Hammond
85: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
87: @*/
88: PetscErrorCode DMDASetBlockFills(DM da,PetscInt *dfill,PetscInt *ofill)
89: {
90: DM_DA *dd = (DM_DA*)da->data;
94: DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
95: DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);
96: return(0);
97: }
102: PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,const MatType mtype,ISColoring *coloring)
103: {
104: PetscErrorCode ierr;
105: PetscInt dim,m,n,p,nc;
106: DMDABoundaryType bx,by,bz;
107: MPI_Comm comm;
108: PetscMPIInt size;
109: PetscBool isBAIJ;
110: DM_DA *dd = (DM_DA*)da->data;
113: /*
114: m
115: ------------------------------------------------------
116: | |
117: | |
118: | ---------------------- |
119: | | | |
120: n | yn | | |
121: | | | |
122: | .--------------------- |
123: | (xs,ys) xn |
124: | . |
125: | (gxs,gys) |
126: | |
127: -----------------------------------------------------
128: */
130: /*
131: nc - number of components per grid point
132: col - number of colors needed in one direction for single component problem
133:
134: */
135: DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);
137: PetscObjectGetComm((PetscObject)da,&comm);
138: MPI_Comm_size(comm,&size);
139: if (ctype == IS_COLORING_GHOSTED){
140: if (size == 1) {
141: ctype = IS_COLORING_GLOBAL;
142: } else if (dim > 1){
143: if ((m==1 && bx == DMDA_BOUNDARY_PERIODIC) || (n==1 && by == DMDA_BOUNDARY_PERIODIC) || (p==1 && bz == DMDA_BOUNDARY_PERIODIC)){
144: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"IS_COLORING_GHOSTED cannot be used for periodic boundary condition having both ends of the domain on the same process");
145: }
146: }
147: }
149: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
150: matrices is for the blocks, not the individual matrix elements */
151: PetscStrcmp(mtype,MATBAIJ,&isBAIJ);
152: if (!isBAIJ) {PetscStrcmp(mtype,MATMPIBAIJ,&isBAIJ);}
153: if (!isBAIJ) {PetscStrcmp(mtype,MATSEQBAIJ,&isBAIJ);}
154: if (isBAIJ) {
155: dd->w = 1;
156: dd->xs = dd->xs/nc;
157: dd->xe = dd->xe/nc;
158: dd->Xs = dd->Xs/nc;
159: dd->Xe = dd->Xe/nc;
160: }
162: /*
163: We do not provide a getcoloring function in the DMDA operations because
164: the basic DMDA does not know about matrices. We think of DMDA as being more
165: more low-level then matrices.
166: */
167: if (dim == 1) {
168: DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
169: } else if (dim == 2) {
170: DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
171: } else if (dim == 3) {
172: DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
173: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
174: if (isBAIJ) {
175: dd->w = nc;
176: dd->xs = dd->xs*nc;
177: dd->xe = dd->xe*nc;
178: dd->Xs = dd->Xs*nc;
179: dd->Xe = dd->Xe*nc;
180: }
181: return(0);
182: }
184: /* ---------------------------------------------------------------------------------*/
188: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
189: {
190: PetscErrorCode ierr;
191: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
192: PetscInt ncolors;
193: MPI_Comm comm;
194: DMDABoundaryType bx,by;
195: DMDAStencilType st;
196: ISColoringValue *colors;
197: DM_DA *dd = (DM_DA*)da->data;
200: /*
201: nc - number of components per grid point
202: col - number of colors needed in one direction for single component problem
203:
204: */
205: DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
206: col = 2*s + 1;
207: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
208: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
209: PetscObjectGetComm((PetscObject)da,&comm);
211: /* special case as taught to us by Paul Hovland */
212: if (st == DMDA_STENCIL_STAR && s == 1) {
213: DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
214: } else {
216: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)){
217: SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\
218: by 2*stencil_width + 1 (%d)\n", m, col);
219: }
220: if (by == DMDA_BOUNDARY_PERIODIC && (n % col)){
221: SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\
222: by 2*stencil_width + 1 (%d)\n", n, col);
223: }
224: if (ctype == IS_COLORING_GLOBAL) {
225: if (!dd->localcoloring) {
226: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
227: ii = 0;
228: for (j=ys; j<ys+ny; j++) {
229: for (i=xs; i<xs+nx; i++) {
230: for (k=0; k<nc; k++) {
231: colors[ii++] = k + nc*((i % col) + col*(j % col));
232: }
233: }
234: }
235: ncolors = nc + nc*(col-1 + col*(col-1));
236: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);
237: }
238: *coloring = dd->localcoloring;
239: } else if (ctype == IS_COLORING_GHOSTED) {
240: if (!dd->ghostedcoloring) {
241: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
242: ii = 0;
243: for (j=gys; j<gys+gny; j++) {
244: for (i=gxs; i<gxs+gnx; i++) {
245: for (k=0; k<nc; k++) {
246: /* the complicated stuff is to handle periodic boundaries */
247: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
248: }
249: }
250: }
251: ncolors = nc + nc*(col - 1 + col*(col-1));
252: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);
253: /* PetscIntView(ncolors,(PetscInt *)colors,0); */
255: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
256: }
257: *coloring = dd->ghostedcoloring;
258: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
259: }
260: ISColoringReference(*coloring);
261: return(0);
262: }
264: /* ---------------------------------------------------------------------------------*/
268: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
269: {
270: PetscErrorCode ierr;
271: 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;
272: PetscInt ncolors;
273: MPI_Comm comm;
274: DMDABoundaryType bx,by,bz;
275: DMDAStencilType st;
276: ISColoringValue *colors;
277: DM_DA *dd = (DM_DA*)da->data;
280: /*
281: nc - number of components per grid point
282: col - number of colors needed in one direction for single component problem
283:
284: */
285: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
286: col = 2*s + 1;
287: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)){
288: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
289: by 2*stencil_width + 1\n");
290: }
291: if (by == DMDA_BOUNDARY_PERIODIC && (n % col)){
292: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
293: by 2*stencil_width + 1\n");
294: }
295: if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)){
296: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
297: by 2*stencil_width + 1\n");
298: }
300: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
301: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
302: PetscObjectGetComm((PetscObject)da,&comm);
304: /* create the coloring */
305: if (ctype == IS_COLORING_GLOBAL) {
306: if (!dd->localcoloring) {
307: PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);
308: ii = 0;
309: for (k=zs; k<zs+nz; k++) {
310: for (j=ys; j<ys+ny; j++) {
311: for (i=xs; i<xs+nx; i++) {
312: for (l=0; l<nc; l++) {
313: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
314: }
315: }
316: }
317: }
318: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
319: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,&dd->localcoloring);
320: }
321: *coloring = dd->localcoloring;
322: } else if (ctype == IS_COLORING_GHOSTED) {
323: if (!dd->ghostedcoloring) {
324: PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);
325: ii = 0;
326: for (k=gzs; k<gzs+gnz; k++) {
327: for (j=gys; j<gys+gny; j++) {
328: for (i=gxs; i<gxs+gnx; i++) {
329: for (l=0; l<nc; l++) {
330: /* the complicated stuff is to handle periodic boundaries */
331: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
332: }
333: }
334: }
335: }
336: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
337: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,&dd->ghostedcoloring);
338: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
339: }
340: *coloring = dd->ghostedcoloring;
341: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
342: ISColoringReference(*coloring);
343: return(0);
344: }
346: /* ---------------------------------------------------------------------------------*/
350: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
351: {
352: PetscErrorCode ierr;
353: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
354: PetscInt ncolors;
355: MPI_Comm comm;
356: DMDABoundaryType bx;
357: ISColoringValue *colors;
358: DM_DA *dd = (DM_DA*)da->data;
361: /*
362: nc - number of components per grid point
363: col - number of colors needed in one direction for single component problem
364:
365: */
366: DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
367: col = 2*s + 1;
369: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)) {
370: SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\
371: by 2*stencil_width + 1 %d\n",(int)m,(int)col);
372: }
374: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
375: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
376: PetscObjectGetComm((PetscObject)da,&comm);
378: /* create the coloring */
379: if (ctype == IS_COLORING_GLOBAL) {
380: if (!dd->localcoloring) {
381: PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);
382: i1 = 0;
383: for (i=xs; i<xs+nx; i++) {
384: for (l=0; l<nc; l++) {
385: colors[i1++] = l + nc*(i % col);
386: }
387: }
388: ncolors = nc + nc*(col-1);
389: ISColoringCreate(comm,ncolors,nc*nx,colors,&dd->localcoloring);
390: }
391: *coloring = dd->localcoloring;
392: } else if (ctype == IS_COLORING_GHOSTED) {
393: if (!dd->ghostedcoloring) {
394: PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);
395: i1 = 0;
396: for (i=gxs; i<gxs+gnx; i++) {
397: for (l=0; l<nc; l++) {
398: /* the complicated stuff is to handle periodic boundaries */
399: colors[i1++] = l + nc*(SetInRange(i,m) % col);
400: }
401: }
402: ncolors = nc + nc*(col-1);
403: ISColoringCreate(comm,ncolors,nc*gnx,colors,&dd->ghostedcoloring);
404: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
405: }
406: *coloring = dd->ghostedcoloring;
407: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
408: ISColoringReference(*coloring);
409: return(0);
410: }
414: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
415: {
416: PetscErrorCode ierr;
417: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
418: PetscInt ncolors;
419: MPI_Comm comm;
420: DMDABoundaryType bx,by;
421: ISColoringValue *colors;
422: DM_DA *dd = (DM_DA*)da->data;
425: /*
426: nc - number of components per grid point
427: col - number of colors needed in one direction for single component problem
428:
429: */
430: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
431: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
432: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
433: PetscObjectGetComm((PetscObject)da,&comm);
435: if (bx == DMDA_BOUNDARY_PERIODIC && (m % 5)){
436: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
437: by 5\n");
438: }
439: if (by == DMDA_BOUNDARY_PERIODIC && (n % 5)){
440: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
441: by 5\n");
442: }
444: /* create the coloring */
445: if (ctype == IS_COLORING_GLOBAL) {
446: if (!dd->localcoloring) {
447: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
448: ii = 0;
449: for (j=ys; j<ys+ny; j++) {
450: for (i=xs; i<xs+nx; i++) {
451: for (k=0; k<nc; k++) {
452: colors[ii++] = k + nc*((3*j+i) % 5);
453: }
454: }
455: }
456: ncolors = 5*nc;
457: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,&dd->localcoloring);
458: }
459: *coloring = dd->localcoloring;
460: } else if (ctype == IS_COLORING_GHOSTED) {
461: if (!dd->ghostedcoloring) {
462: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
463: ii = 0;
464: for (j=gys; j<gys+gny; j++) {
465: for (i=gxs; i<gxs+gnx; i++) {
466: for (k=0; k<nc; k++) {
467: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
468: }
469: }
470: }
471: ncolors = 5*nc;
472: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,&dd->ghostedcoloring);
473: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
474: }
475: *coloring = dd->ghostedcoloring;
476: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
477: return(0);
478: }
480: /* =========================================================================== */
481: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
482: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
483: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
484: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
485: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
486: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
487: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
488: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
489: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
493: /*@
494: MatSetDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix
496: Logically Collective on Mat
498: Input Parameters:
499: + mat - the matrix
500: - da - the da
502: Level: intermediate
504: @*/
505: PetscErrorCode MatSetDM(Mat mat,DM da)
506: {
512: PetscTryMethod(mat,"MatSetDM_C",(Mat,DM),(mat,da));
513: return(0);
514: }
516: EXTERN_C_BEGIN
519: PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer)
520: {
521: DM da;
523: const char *prefix;
524: Mat Anatural;
525: AO ao;
526: PetscInt rstart,rend,*petsc,i;
527: IS is;
528: MPI_Comm comm;
529: PetscViewerFormat format;
532: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
533: PetscViewerGetFormat(viewer,&format);
534: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
536: PetscObjectGetComm((PetscObject)A,&comm);
537: PetscObjectQuery((PetscObject)A,"DM",(PetscObject*)&da);
538: if (!da) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
540: DMDAGetAO(da,&ao);
541: MatGetOwnershipRange(A,&rstart,&rend);
542: PetscMalloc((rend-rstart)*sizeof(PetscInt),&petsc);
543: for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
544: AOApplicationToPetsc(ao,rend-rstart,petsc);
545: ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);
547: /* call viewer on natural ordering */
548: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
549: ISDestroy(&is);
550: PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
551: PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
552: PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
553: MatView(Anatural,viewer);
554: MatDestroy(&Anatural);
555: return(0);
556: }
557: EXTERN_C_END
559: EXTERN_C_BEGIN
562: PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer)
563: {
564: DM da;
566: Mat Anatural,Aapp;
567: AO ao;
568: PetscInt rstart,rend,*app,i;
569: IS is;
570: MPI_Comm comm;
573: PetscObjectGetComm((PetscObject)A,&comm);
574: PetscObjectQuery((PetscObject)A,"DM",(PetscObject*)&da);
575: if (!da) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
577: /* Load the matrix in natural ordering */
578: MatCreate(((PetscObject)A)->comm,&Anatural);
579: MatSetType(Anatural,((PetscObject)A)->type_name);
580: MatSetSizes(Anatural,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
581: MatLoad(Anatural,viewer);
583: /* Map natural ordering to application ordering and create IS */
584: DMDAGetAO(da,&ao);
585: MatGetOwnershipRange(Anatural,&rstart,&rend);
586: PetscMalloc((rend-rstart)*sizeof(PetscInt),&app);
587: for (i=rstart; i<rend; i++) app[i-rstart] = i;
588: AOPetscToApplication(ao,rend-rstart,app);
589: ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);
591: /* Do permutation and replace header */
592: MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
593: MatHeaderReplace(A,Aapp);
594: ISDestroy(&is);
595: MatDestroy(&Anatural);
596: return(0);
597: }
598: EXTERN_C_END
602: PetscErrorCode DMCreateMatrix_DA(DM da, const MatType mtype,Mat *J)
603: {
605: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
606: Mat A;
607: MPI_Comm comm;
608: const MatType Atype;
609: PetscSection section, sectionGlobal;
610: void (*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL;
611: MatType ttype[256];
612: PetscBool flg;
613: PetscMPIInt size;
614: DM_DA *dd = (DM_DA*)da->data;
617: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
618: MatInitializePackage(PETSC_NULL);
619: #endif
620: if (!mtype) mtype = MATAIJ;
621: PetscStrcpy((char*)ttype,mtype);
622: PetscOptionsBegin(((PetscObject)da)->comm,((PetscObject)da)->prefix,"DMDA options","Mat");
623: PetscOptionsList("-dm_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)ttype,256,&flg);
624: PetscOptionsEnd();
626: DMGetDefaultSection(da, §ion);
627: if (section) {
628: PetscInt bs = -1;
629: PetscInt localSize;
630: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
632: DMGetDefaultGlobalSection(da, §ionGlobal);
633: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
634: MatCreate(((PetscObject) da)->comm, J);
635: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
636: MatSetType(*J, mtype);
637: MatSetFromOptions(*J);
638: PetscStrcmp(mtype, MATSHELL, &isShell);
639: PetscStrcmp(mtype, MATBAIJ, &isBlock);
640: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
641: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
642: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
643: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
644: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
645: /* Check for symmetric storage */
646: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
647: if (isSymmetric) {
648: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
649: }
650: if (!isShell) {
651: /* PetscBool fillMatrix = (PetscBool) !da->prealloc_only; */
652: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
654: if (bs < 0) {
655: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
656: PetscInt pStart, pEnd, p, dof;
658: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
659: for(p = pStart; p < pEnd; ++p) {
660: PetscSectionGetDof(sectionGlobal, p, &dof);
661: if (dof) {
662: bs = dof;
663: break;
664: }
665: }
666: } else {
667: bs = 1;
668: }
669: /* Must have same blocksize on all procs (some might have no points) */
670: bsLocal = bs;
671: MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, ((PetscObject) da)->comm);
672: }
673: PetscMalloc4(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz, localSize/bs, PetscInt, &dnzu, localSize/bs, PetscInt, &onzu);
674: PetscMemzero(dnz, localSize/bs * sizeof(PetscInt));
675: PetscMemzero(onz, localSize/bs * sizeof(PetscInt));
676: PetscMemzero(dnzu, localSize/bs * sizeof(PetscInt));
677: PetscMemzero(onzu, localSize/bs * sizeof(PetscInt));
678: /* DMComplexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
679: PetscFree4(dnz, onz, dnzu, onzu);
680: }
681: }
682: /*
683: m
684: ------------------------------------------------------
685: | |
686: | |
687: | ---------------------- |
688: | | | |
689: n | ny | | |
690: | | | |
691: | .--------------------- |
692: | (xs,ys) nx |
693: | . |
694: | (gxs,gys) |
695: | |
696: -----------------------------------------------------
697: */
699: /*
700: nc - number of components per grid point
701: col - number of colors needed in one direction for single component problem
702:
703: */
704: DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0);
705: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
706: PetscObjectGetComm((PetscObject)da,&comm);
707: MatCreate(comm,&A);
708: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
709: MatSetType(A,(const MatType)ttype);
710: MatSetDM(A,da);
711: MatSetFromOptions(A);
712: MatGetType(A,&Atype);
713: /*
714: We do not provide a getmatrix function in the DMDA operations because
715: the basic DMDA does not know about matrices. We think of DMDA as being more
716: more low-level than matrices. This is kind of cheating but, cause sometimes
717: we think of DMDA has higher level than matrices.
719: We could switch based on Atype (or mtype), but we do not since the
720: specialized setting routines depend only the particular preallocation
721: details of the matrix, not the type itself.
722: */
723: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
724: if (!aij) {
725: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
726: }
727: if (!aij) {
728: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
729: if (!baij) {
730: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
731: }
732: if (!baij){
733: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
734: if (!sbaij) {
735: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
736: }
737: }
738: }
739: if (aij) {
740: if (dim == 1) {
741: DMCreateMatrix_DA_1d_MPIAIJ(da,A);
742: } else if (dim == 2) {
743: if (dd->ofill) {
744: DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
745: } else {
746: DMCreateMatrix_DA_2d_MPIAIJ(da,A);
747: }
748: } else if (dim == 3) {
749: if (dd->ofill) {
750: DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
751: } else {
752: DMCreateMatrix_DA_3d_MPIAIJ(da,A);
753: }
754: }
755: } else if (baij) {
756: if (dim == 2) {
757: DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
758: } else if (dim == 3) {
759: DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
760: } else {
761: SETERRQ3(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension!\n" \
762: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim);
763: }
764: } else if (sbaij) {
765: if (dim == 2) {
766: DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
767: } else if (dim == 3) {
768: DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
769: } else {
770: SETERRQ3(((PetscObject)da)->comm,PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension!\n" \
771: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim);
772: }
773: } else {
774: ISLocalToGlobalMapping ltog,ltogb;
775: DMGetLocalToGlobalMapping(da,<og);
776: DMGetLocalToGlobalMappingBlock(da,<ogb);
777: MatSetUp(A);
778: MatSetLocalToGlobalMapping(A,ltog,ltog);
779: MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);
780: }
781: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
782: MatSetStencil(A,dim,dims,starts,dof);
783: PetscObjectCompose((PetscObject)A,"DM",(PetscObject)da);
784: MPI_Comm_size(comm,&size);
785: if (size > 1) {
786: /* change viewer to display matrix in natural ordering */
787: MatShellSetOperation(A, MATOP_VIEW, (void (*)(void)) MatView_MPI_DA);
788: MatShellSetOperation(A, MATOP_LOAD, (void (*)(void)) MatLoad_MPI_DA);
789: }
790: *J = A;
791: return(0);
792: }
794: /* ---------------------------------------------------------------------------------*/
797: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
798: {
799: PetscErrorCode ierr;
800: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p;
801: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
802: MPI_Comm comm;
803: PetscScalar *values;
804: DMDABoundaryType bx,by;
805: ISLocalToGlobalMapping ltog,ltogb;
806: DMDAStencilType st;
809: /*
810: nc - number of components per grid point
811: col - number of colors needed in one direction for single component problem
812:
813: */
814: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
815: col = 2*s + 1;
816: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
817: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
818: PetscObjectGetComm((PetscObject)da,&comm);
820: PetscMalloc2(nc,PetscInt,&rows,col*col*nc*nc,PetscInt,&cols);
821: DMGetLocalToGlobalMapping(da,<og);
822: DMGetLocalToGlobalMappingBlock(da,<ogb);
823:
824: /* determine the matrix preallocation information */
825: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
826: for (i=xs; i<xs+nx; i++) {
828: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
829: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
831: for (j=ys; j<ys+ny; j++) {
832: slot = i - gxs + gnx*(j - gys);
834: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
835: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
837: cnt = 0;
838: for (k=0; k<nc; k++) {
839: for (l=lstart; l<lend+1; l++) {
840: for (p=pstart; p<pend+1; p++) {
841: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
842: cols[cnt++] = k + nc*(slot + gnx*l + p);
843: }
844: }
845: }
846: rows[k] = k + nc*(slot);
847: }
848: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
849: }
850: }
851: MatSetBlockSize(J,nc);
852: MatSeqAIJSetPreallocation(J,0,dnz);
853: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
854: MatPreallocateFinalize(dnz,onz);
856: MatSetLocalToGlobalMapping(J,ltog,ltog);
857: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
859: /*
860: For each node in the grid: we get the neighbors in the local (on processor ordering
861: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
862: PETSc ordering.
863: */
864: if (!da->prealloc_only) {
865: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
866: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
867: for (i=xs; i<xs+nx; i++) {
868:
869: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
870: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
871:
872: for (j=ys; j<ys+ny; j++) {
873: slot = i - gxs + gnx*(j - gys);
874:
875: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
876: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
878: cnt = 0;
879: for (k=0; k<nc; k++) {
880: for (l=lstart; l<lend+1; l++) {
881: for (p=pstart; p<pend+1; p++) {
882: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
883: cols[cnt++] = k + nc*(slot + gnx*l + p);
884: }
885: }
886: }
887: rows[k] = k + nc*(slot);
888: }
889: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
890: }
891: }
892: PetscFree(values);
893: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
894: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
895: }
896: PetscFree2(rows,cols);
897: return(0);
898: }
902: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
903: {
904: PetscErrorCode ierr;
905: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
906: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p;
907: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
908: DM_DA *dd = (DM_DA*)da->data;
909: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
910: MPI_Comm comm;
911: PetscScalar *values;
912: DMDABoundaryType bx,by;
913: ISLocalToGlobalMapping ltog,ltogb;
914: DMDAStencilType st;
917: /*
918: nc - number of components per grid point
919: col - number of colors needed in one direction for single component problem
920:
921: */
922: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
923: col = 2*s + 1;
924: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
925: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
926: PetscObjectGetComm((PetscObject)da,&comm);
928: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
929: DMGetLocalToGlobalMapping(da,<og);
930: DMGetLocalToGlobalMappingBlock(da,<ogb);
931:
932: /* determine the matrix preallocation information */
933: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
934: for (i=xs; i<xs+nx; i++) {
936: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
937: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
939: for (j=ys; j<ys+ny; j++) {
940: slot = i - gxs + gnx*(j - gys);
942: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
943: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
945: for (k=0; k<nc; k++) {
946: cnt = 0;
947: for (l=lstart; l<lend+1; l++) {
948: for (p=pstart; p<pend+1; p++) {
949: if (l || p) {
950: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
951: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
952: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
953: }
954: } else {
955: if (dfill) {
956: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
957: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
958: } else {
959: for (ifill_col=0; ifill_col<nc; ifill_col++)
960: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
961: }
962: }
963: }
964: }
965: row = k + nc*(slot);
966: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
967: }
968: }
969: }
970: MatSeqAIJSetPreallocation(J,0,dnz);
971: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
972: MatPreallocateFinalize(dnz,onz);
973: MatSetLocalToGlobalMapping(J,ltog,ltog);
974: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
976: /*
977: For each node in the grid: we get the neighbors in the local (on processor ordering
978: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
979: PETSc ordering.
980: */
981: if (!da->prealloc_only) {
982: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
983: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
984: for (i=xs; i<xs+nx; i++) {
985:
986: pstart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
987: pend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
988:
989: for (j=ys; j<ys+ny; j++) {
990: slot = i - gxs + gnx*(j - gys);
991:
992: lstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
993: lend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
995: for (k=0; k<nc; k++) {
996: cnt = 0;
997: for (l=lstart; l<lend+1; l++) {
998: for (p=pstart; p<pend+1; p++) {
999: if (l || p) {
1000: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1001: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
1002: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1003: }
1004: } else {
1005: if (dfill) {
1006: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
1007: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1008: } else {
1009: for (ifill_col=0; ifill_col<nc; ifill_col++)
1010: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1011: }
1012: }
1013: }
1014: }
1015: row = k + nc*(slot);
1016: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1017: }
1018: }
1019: }
1020: PetscFree(values);
1021: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1022: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1023: }
1024: PetscFree(cols);
1025: return(0);
1026: }
1028: /* ---------------------------------------------------------------------------------*/
1032: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1033: {
1034: PetscErrorCode ierr;
1035: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1036: PetscInt m,n,dim,s,*cols = PETSC_NULL,k,nc,*rows = PETSC_NULL,col,cnt,l,p,*dnz = PETSC_NULL,*onz = PETSC_NULL;
1037: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1038: MPI_Comm comm;
1039: PetscScalar *values;
1040: DMDABoundaryType bx,by,bz;
1041: ISLocalToGlobalMapping ltog,ltogb;
1042: DMDAStencilType st;
1045: /*
1046: nc - number of components per grid point
1047: col - number of colors needed in one direction for single component problem
1048:
1049: */
1050: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1051: col = 2*s + 1;
1053: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1054: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1055: PetscObjectGetComm((PetscObject)da,&comm);
1057: PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);
1058: DMGetLocalToGlobalMapping(da,<og);
1059: DMGetLocalToGlobalMappingBlock(da,<ogb);
1061: /* determine the matrix preallocation information */
1062: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1063: for (i=xs; i<xs+nx; i++) {
1064: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1065: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1066: for (j=ys; j<ys+ny; j++) {
1067: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1068: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1069: for (k=zs; k<zs+nz; k++) {
1070: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1071: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1072:
1073: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1074:
1075: cnt = 0;
1076: for (l=0; l<nc; l++) {
1077: for (ii=istart; ii<iend+1; ii++) {
1078: for (jj=jstart; jj<jend+1; jj++) {
1079: for (kk=kstart; kk<kend+1; kk++) {
1080: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1081: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1082: }
1083: }
1084: }
1085: }
1086: rows[l] = l + nc*(slot);
1087: }
1088: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1089: }
1090: }
1091: }
1092: MatSetBlockSize(J,nc);
1093: MatSeqAIJSetPreallocation(J,0,dnz);
1094: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1095: MatPreallocateFinalize(dnz,onz);
1096: MatSetLocalToGlobalMapping(J,ltog,ltog);
1097: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1099: /*
1100: For each node in the grid: we get the neighbors in the local (on processor ordering
1101: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1102: PETSc ordering.
1103: */
1104: if (!da->prealloc_only) {
1105: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1106: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1107: for (i=xs; i<xs+nx; i++) {
1108: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1109: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1110: for (j=ys; j<ys+ny; j++) {
1111: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1112: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1113: for (k=zs; k<zs+nz; k++) {
1114: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1115: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1116:
1117: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1118:
1119: cnt = 0;
1120: for (l=0; l<nc; l++) {
1121: for (ii=istart; ii<iend+1; ii++) {
1122: for (jj=jstart; jj<jend+1; jj++) {
1123: for (kk=kstart; kk<kend+1; kk++) {
1124: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1125: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1126: }
1127: }
1128: }
1129: }
1130: rows[l] = l + nc*(slot);
1131: }
1132: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1133: }
1134: }
1135: }
1136: PetscFree(values);
1137: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1138: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1139: }
1140: PetscFree2(rows,cols);
1141: return(0);
1142: }
1144: /* ---------------------------------------------------------------------------------*/
1148: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1149: {
1150: PetscErrorCode ierr;
1151: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1152: PetscInt m,dim,s,*cols = PETSC_NULL,nc,*rows = PETSC_NULL,col,cnt,l;
1153: PetscInt istart,iend;
1154: PetscScalar *values;
1155: DMDABoundaryType bx;
1156: ISLocalToGlobalMapping ltog,ltogb;
1159: /*
1160: nc - number of components per grid point
1161: col - number of colors needed in one direction for single component problem
1162:
1163: */
1164: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1165: col = 2*s + 1;
1167: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1168: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1170: MatSetBlockSize(J,nc);
1171: MatSeqAIJSetPreallocation(J,col*nc,0);
1172: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1173: PetscMalloc2(nc,PetscInt,&rows,col*nc*nc,PetscInt,&cols);
1174:
1175: DMGetLocalToGlobalMapping(da,<og);
1176: DMGetLocalToGlobalMappingBlock(da,<ogb);
1177: MatSetLocalToGlobalMapping(J,ltog,ltog);
1178: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1179:
1180: /*
1181: For each node in the grid: we get the neighbors in the local (on processor ordering
1182: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1183: PETSc ordering.
1184: */
1185: if (!da->prealloc_only) {
1186: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
1187: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
1188: for (i=xs; i<xs+nx; i++) {
1189: istart = PetscMax(-s,gxs - i);
1190: iend = PetscMin(s,gxs + gnx - i - 1);
1191: slot = i - gxs;
1192:
1193: cnt = 0;
1194: for (l=0; l<nc; l++) {
1195: for (i1=istart; i1<iend+1; i1++) {
1196: cols[cnt++] = l + nc*(slot + i1);
1197: }
1198: rows[l] = l + nc*(slot);
1199: }
1200: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1201: }
1202: PetscFree(values);
1203: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1204: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1205: }
1206: PetscFree2(rows,cols);
1207: return(0);
1208: }
1212: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1213: {
1214: PetscErrorCode ierr;
1215: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1216: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1217: PetscInt istart,iend,jstart,jend,ii,jj;
1218: MPI_Comm comm;
1219: PetscScalar *values;
1220: DMDABoundaryType bx,by;
1221: DMDAStencilType st;
1222: ISLocalToGlobalMapping ltog,ltogb;
1225: /*
1226: nc - number of components per grid point
1227: col - number of colors needed in one direction for single component problem
1228: */
1229: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1230: col = 2*s + 1;
1232: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1233: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1234: PetscObjectGetComm((PetscObject)da,&comm);
1236: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1238: DMGetLocalToGlobalMapping(da,<og);
1239: DMGetLocalToGlobalMappingBlock(da,<ogb);
1241: /* determine the matrix preallocation information */
1242: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1243: for (i=xs; i<xs+nx; i++) {
1244: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1245: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1246: for (j=ys; j<ys+ny; j++) {
1247: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1248: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1249: slot = i - gxs + gnx*(j - gys);
1251: /* Find block columns in block row */
1252: cnt = 0;
1253: for (ii=istart; ii<iend+1; ii++) {
1254: for (jj=jstart; jj<jend+1; jj++) {
1255: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1256: cols[cnt++] = slot + ii + gnx*jj;
1257: }
1258: }
1259: }
1260: MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);
1261: }
1262: }
1263: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1264: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1265: MatPreallocateFinalize(dnz,onz);
1267: MatSetLocalToGlobalMapping(J,ltog,ltog);
1268: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1270: /*
1271: For each node in the grid: we get the neighbors in the local (on processor ordering
1272: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1273: PETSc ordering.
1274: */
1275: if (!da->prealloc_only) {
1276: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1277: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1278: for (i=xs; i<xs+nx; i++) {
1279: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1280: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1281: for (j=ys; j<ys+ny; j++) {
1282: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1283: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1284: slot = i - gxs + gnx*(j - gys);
1285: cnt = 0;
1286: for (ii=istart; ii<iend+1; ii++) {
1287: for (jj=jstart; jj<jend+1; jj++) {
1288: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1289: cols[cnt++] = slot + ii + gnx*jj;
1290: }
1291: }
1292: }
1293: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1294: }
1295: }
1296: PetscFree(values);
1297: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1298: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1299: }
1300: PetscFree(cols);
1301: return(0);
1302: }
1306: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1307: {
1308: PetscErrorCode ierr;
1309: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1310: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1311: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1312: MPI_Comm comm;
1313: PetscScalar *values;
1314: DMDABoundaryType bx,by,bz;
1315: DMDAStencilType st;
1316: ISLocalToGlobalMapping ltog,ltogb;
1319: /*
1320: nc - number of components per grid point
1321: col - number of colors needed in one direction for single component problem
1322:
1323: */
1324: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1325: col = 2*s + 1;
1327: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1328: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1329: PetscObjectGetComm((PetscObject)da,&comm);
1331: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1333: DMGetLocalToGlobalMapping(da,<og);
1334: DMGetLocalToGlobalMappingBlock(da,<ogb);
1336: /* determine the matrix preallocation information */
1337: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1338: for (i=xs; i<xs+nx; i++) {
1339: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1340: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1341: for (j=ys; j<ys+ny; j++) {
1342: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1343: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1344: for (k=zs; k<zs+nz; k++) {
1345: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1346: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1348: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1350: /* Find block columns in block row */
1351: cnt = 0;
1352: for (ii=istart; ii<iend+1; ii++) {
1353: for (jj=jstart; jj<jend+1; jj++) {
1354: for (kk=kstart; kk<kend+1; kk++) {
1355: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1356: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1357: }
1358: }
1359: }
1360: }
1361: MatPreallocateSetLocal(ltogb,1,&slot,ltogb,cnt,cols,dnz,onz);
1362: }
1363: }
1364: }
1365: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1366: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1367: MatPreallocateFinalize(dnz,onz);
1369: MatSetLocalToGlobalMapping(J,ltog,ltog);
1370: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1372: /*
1373: For each node in the grid: we get the neighbors in the local (on processor ordering
1374: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1375: PETSc ordering.
1376: */
1377: if (!da->prealloc_only) {
1378: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1379: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1380: for (i=xs; i<xs+nx; i++) {
1381: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1382: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1383: for (j=ys; j<ys+ny; j++) {
1384: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1385: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1386: for (k=zs; k<zs+nz; k++) {
1387: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1388: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1389:
1390: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1391:
1392: cnt = 0;
1393: for (ii=istart; ii<iend+1; ii++) {
1394: for (jj=jstart; jj<jend+1; jj++) {
1395: for (kk=kstart; kk<kend+1; kk++) {
1396: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1397: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1398: }
1399: }
1400: }
1401: }
1402: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1403: }
1404: }
1405: }
1406: PetscFree(values);
1407: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1408: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1409: }
1410: PetscFree(cols);
1411: return(0);
1412: }
1416: /*
1417: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1418: identify in the local ordering with periodic domain.
1419: */
1420: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1421: {
1423: PetscInt i,n;
1426: ISLocalToGlobalMappingApply(ltog,1,row,row);
1427: ISLocalToGlobalMappingApply(ltog,*cnt,col,col);
1428: for (i=0,n=0; i<*cnt; i++) {
1429: if (col[i] >= *row) col[n++] = col[i];
1430: }
1431: *cnt = n;
1432: return(0);
1433: }
1437: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(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,nc,col,cnt,*dnz,*onz;
1442: PetscInt istart,iend,jstart,jend,ii,jj;
1443: MPI_Comm comm;
1444: PetscScalar *values;
1445: DMDABoundaryType bx,by;
1446: DMDAStencilType st;
1447: ISLocalToGlobalMapping ltog,ltogb;
1450: /*
1451: nc - number of components per grid point
1452: col - number of colors needed in one direction for single component problem
1453: */
1454: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1455: col = 2*s + 1;
1457: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1458: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1459: PetscObjectGetComm((PetscObject)da,&comm);
1461: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1463: DMGetLocalToGlobalMapping(da,<og);
1464: DMGetLocalToGlobalMappingBlock(da,<ogb);
1466: /* determine the matrix preallocation information */
1467: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1468: for (i=xs; i<xs+nx; i++) {
1469: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1470: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1471: for (j=ys; j<ys+ny; j++) {
1472: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1473: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1474: slot = i - gxs + gnx*(j - gys);
1476: /* Find block columns in block row */
1477: cnt = 0;
1478: for (ii=istart; ii<iend+1; ii++) {
1479: for (jj=jstart; jj<jend+1; jj++) {
1480: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1481: cols[cnt++] = slot + ii + gnx*jj;
1482: }
1483: }
1484: }
1485: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1486: MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);
1487: }
1488: }
1489: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1490: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1491: MatPreallocateFinalize(dnz,onz);
1493: MatSetLocalToGlobalMapping(J,ltog,ltog);
1494: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1496: /*
1497: For each node in the grid: we get the neighbors in the local (on processor ordering
1498: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1499: PETSc ordering.
1500: */
1501: if (!da->prealloc_only) {
1502: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1503: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1504: for (i=xs; i<xs+nx; i++) {
1505: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1506: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1507: for (j=ys; j<ys+ny; j++) {
1508: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1509: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1510: slot = i - gxs + gnx*(j - gys);
1512: /* Find block columns in block row */
1513: cnt = 0;
1514: for (ii=istart; ii<iend+1; ii++) {
1515: for (jj=jstart; jj<jend+1; jj++) {
1516: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1517: cols[cnt++] = slot + ii + gnx*jj;
1518: }
1519: }
1520: }
1521: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1522: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1523: }
1524: }
1525: PetscFree(values);
1526: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1527: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1528: }
1529: PetscFree(cols);
1530: return(0);
1531: }
1535: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1536: {
1537: PetscErrorCode ierr;
1538: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1539: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1540: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1541: MPI_Comm comm;
1542: PetscScalar *values;
1543: DMDABoundaryType bx,by,bz;
1544: DMDAStencilType st;
1545: ISLocalToGlobalMapping ltog,ltogb;
1548: /*
1549: nc - number of components per grid point
1550: col - number of colors needed in one direction for single component problem
1551: */
1552: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1553: col = 2*s + 1;
1555: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1556: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1557: PetscObjectGetComm((PetscObject)da,&comm);
1559: /* create the matrix */
1560: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1562: DMGetLocalToGlobalMapping(da,<og);
1563: DMGetLocalToGlobalMappingBlock(da,<ogb);
1565: /* determine the matrix preallocation information */
1566: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1567: for (i=xs; i<xs+nx; i++) {
1568: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1569: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1570: for (j=ys; j<ys+ny; j++) {
1571: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1572: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1573: for (k=zs; k<zs+nz; k++) {
1574: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1575: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1577: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1579: /* Find block columns in block row */
1580: cnt = 0;
1581: for (ii=istart; ii<iend+1; ii++) {
1582: for (jj=jstart; jj<jend+1; jj++) {
1583: for (kk=kstart; kk<kend+1; kk++) {
1584: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1585: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1586: }
1587: }
1588: }
1589: }
1590: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1591: MatPreallocateSymmetricSet(slot,cnt,cols,dnz,onz);
1592: }
1593: }
1594: }
1595: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1596: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1597: MatPreallocateFinalize(dnz,onz);
1599: MatSetLocalToGlobalMapping(J,ltog,ltog);
1600: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1602: /*
1603: For each node in the grid: we get the neighbors in the local (on processor ordering
1604: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1605: PETSc ordering.
1606: */
1607: if (!da->prealloc_only) {
1608: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1609: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1610: for (i=xs; i<xs+nx; i++) {
1611: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1612: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1613: for (j=ys; j<ys+ny; j++) {
1614: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1615: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1616: for (k=zs; k<zs+nz; k++) {
1617: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1618: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1619:
1620: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1621:
1622: cnt = 0;
1623: for (ii=istart; ii<iend+1; ii++) {
1624: for (jj=jstart; jj<jend+1; jj++) {
1625: for (kk=kstart; kk<kend+1; kk++) {
1626: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1627: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1628: }
1629: }
1630: }
1631: }
1632: L2GFilterUpperTriangular(ltogb,&slot,&cnt,cols);
1633: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1634: }
1635: }
1636: }
1637: PetscFree(values);
1638: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1639: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1640: }
1641: PetscFree(cols);
1642: return(0);
1643: }
1645: /* ---------------------------------------------------------------------------------*/
1649: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
1650: {
1651: PetscErrorCode ierr;
1652: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1653: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz;
1654: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1655: DM_DA *dd = (DM_DA*)da->data;
1656: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
1657: MPI_Comm comm;
1658: PetscScalar *values;
1659: DMDABoundaryType bx,by,bz;
1660: ISLocalToGlobalMapping ltog,ltogb;
1661: DMDAStencilType st;
1664: /*
1665: nc - number of components per grid point
1666: col - number of colors needed in one direction for single component problem
1667:
1668: */
1669: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1670: col = 2*s + 1;
1671: if (bx == DMDA_BOUNDARY_PERIODIC && (m % col)){
1672: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1673: by 2*stencil_width + 1\n");
1674: }
1675: if (by == DMDA_BOUNDARY_PERIODIC && (n % col)){
1676: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1677: by 2*stencil_width + 1\n");
1678: }
1679: if (bz == DMDA_BOUNDARY_PERIODIC && (p % col)){
1680: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1681: by 2*stencil_width + 1\n");
1682: }
1684: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1685: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1686: PetscObjectGetComm((PetscObject)da,&comm);
1688: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
1689: DMGetLocalToGlobalMapping(da,<og);
1690: DMGetLocalToGlobalMappingBlock(da,<ogb);
1692: /* determine the matrix preallocation information */
1693: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1696: for (i=xs; i<xs+nx; i++) {
1697: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1698: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1699: for (j=ys; j<ys+ny; j++) {
1700: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1701: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1702: for (k=zs; k<zs+nz; k++) {
1703: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1704: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1705:
1706: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1707:
1708: for (l=0; l<nc; l++) {
1709: cnt = 0;
1710: for (ii=istart; ii<iend+1; ii++) {
1711: for (jj=jstart; jj<jend+1; jj++) {
1712: for (kk=kstart; kk<kend+1; kk++) {
1713: if (ii || jj || kk) {
1714: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1715: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1716: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1717: }
1718: } else {
1719: if (dfill) {
1720: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1721: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1722: } else {
1723: for (ifill_col=0; ifill_col<nc; ifill_col++)
1724: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1725: }
1726: }
1727: }
1728: }
1729: }
1730: row = l + nc*(slot);
1731: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1732: }
1733: }
1734: }
1735: }
1736: MatSeqAIJSetPreallocation(J,0,dnz);
1737: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1738: MatPreallocateFinalize(dnz,onz);
1739: MatSetLocalToGlobalMapping(J,ltog,ltog);
1740: MatSetLocalToGlobalMappingBlock(J,ltogb,ltogb);
1742: /*
1743: For each node in the grid: we get the neighbors in the local (on processor ordering
1744: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1745: PETSc ordering.
1746: */
1747: if (!da->prealloc_only) {
1748: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1749: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1750: for (i=xs; i<xs+nx; i++) {
1751: istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1752: iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1753: for (j=ys; j<ys+ny; j++) {
1754: jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1755: jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1756: for (k=zs; k<zs+nz; k++) {
1757: kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1758: kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1759:
1760: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1761:
1762: for (l=0; l<nc; l++) {
1763: cnt = 0;
1764: for (ii=istart; ii<iend+1; ii++) {
1765: for (jj=jstart; jj<jend+1; jj++) {
1766: for (kk=kstart; kk<kend+1; kk++) {
1767: if (ii || jj || kk) {
1768: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1769: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1770: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1771: }
1772: } else {
1773: if (dfill) {
1774: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1775: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1776: } else {
1777: for (ifill_col=0; ifill_col<nc; ifill_col++)
1778: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1779: }
1780: }
1781: }
1782: }
1783: }
1784: row = l + nc*(slot);
1785: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1786: }
1787: }
1788: }
1789: }
1790: PetscFree(values);
1791: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1792: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1793: }
1794: PetscFree(cols);
1795: return(0);
1796: }