Actual source code: fdda.c
petsc-3.7.7 2017-09-25
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: PetscMalloc1(nz + w + 1,&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: PetscCalloc1(dd->w,&dd->ofillcols);
103: for (i=0; i<dd->w; i++) {
104: for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
105: }
106: for (i=0; i<dd->w; i++) {
107: if (dd->ofillcols[i]) {
108: dd->ofillcols[i] = cnt++;
109: }
110: }
111: return(0);
112: }
117: PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
118: {
119: PetscErrorCode ierr;
120: PetscInt dim,m,n,p,nc;
121: DMBoundaryType bx,by,bz;
122: MPI_Comm comm;
123: PetscMPIInt size;
124: PetscBool isBAIJ;
125: DM_DA *dd = (DM_DA*)da->data;
128: /*
129: m
130: ------------------------------------------------------
131: | |
132: | |
133: | ---------------------- |
134: | | | |
135: n | yn | | |
136: | | | |
137: | .--------------------- |
138: | (xs,ys) xn |
139: | . |
140: | (gxs,gys) |
141: | |
142: -----------------------------------------------------
143: */
145: /*
146: nc - number of components per grid point
147: col - number of colors needed in one direction for single component problem
149: */
150: DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);
152: PetscObjectGetComm((PetscObject)da,&comm);
153: MPI_Comm_size(comm,&size);
154: if (ctype == IS_COLORING_GHOSTED) {
155: if (size == 1) {
156: ctype = IS_COLORING_GLOBAL;
157: } else if (dim > 1) {
158: if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
159: 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");
160: }
161: }
162: }
164: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
165: matrices is for the blocks, not the individual matrix elements */
166: PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);
167: if (!isBAIJ) {PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);}
168: if (!isBAIJ) {PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);}
169: if (isBAIJ) {
170: dd->w = 1;
171: dd->xs = dd->xs/nc;
172: dd->xe = dd->xe/nc;
173: dd->Xs = dd->Xs/nc;
174: dd->Xe = dd->Xe/nc;
175: }
177: /*
178: We do not provide a getcoloring function in the DMDA operations because
179: the basic DMDA does not know about matrices. We think of DMDA as being more
180: more low-level then matrices.
181: */
182: if (dim == 1) {
183: DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
184: } else if (dim == 2) {
185: DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
186: } else if (dim == 3) {
187: DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
188: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
189: if (isBAIJ) {
190: dd->w = nc;
191: dd->xs = dd->xs*nc;
192: dd->xe = dd->xe*nc;
193: dd->Xs = dd->Xs*nc;
194: dd->Xe = dd->Xe*nc;
195: }
196: return(0);
197: }
199: /* ---------------------------------------------------------------------------------*/
203: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
204: {
205: PetscErrorCode ierr;
206: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
207: PetscInt ncolors;
208: MPI_Comm comm;
209: DMBoundaryType bx,by;
210: DMDAStencilType st;
211: ISColoringValue *colors;
212: DM_DA *dd = (DM_DA*)da->data;
215: /*
216: nc - number of components per grid point
217: col - number of colors needed in one direction for single component problem
219: */
220: DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
221: col = 2*s + 1;
222: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
223: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
224: PetscObjectGetComm((PetscObject)da,&comm);
226: /* special case as taught to us by Paul Hovland */
227: if (st == DMDA_STENCIL_STAR && s == 1) {
228: DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
229: } else {
231: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\
232: by 2*stencil_width + 1 (%d)\n", m, col);
233: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\
234: by 2*stencil_width + 1 (%d)\n", n, col);
235: if (ctype == IS_COLORING_GLOBAL) {
236: if (!dd->localcoloring) {
237: PetscMalloc1(nc*nx*ny,&colors);
238: ii = 0;
239: for (j=ys; j<ys+ny; j++) {
240: for (i=xs; i<xs+nx; i++) {
241: for (k=0; k<nc; k++) {
242: colors[ii++] = k + nc*((i % col) + col*(j % col));
243: }
244: }
245: }
246: ncolors = nc + nc*(col-1 + col*(col-1));
247: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
248: }
249: *coloring = dd->localcoloring;
250: } else if (ctype == IS_COLORING_GHOSTED) {
251: if (!dd->ghostedcoloring) {
252: PetscMalloc1(nc*gnx*gny,&colors);
253: ii = 0;
254: for (j=gys; j<gys+gny; j++) {
255: for (i=gxs; i<gxs+gnx; i++) {
256: for (k=0; k<nc; k++) {
257: /* the complicated stuff is to handle periodic boundaries */
258: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
259: }
260: }
261: }
262: ncolors = nc + nc*(col - 1 + col*(col-1));
263: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
264: /* PetscIntView(ncolors,(PetscInt*)colors,0); */
266: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
267: }
268: *coloring = dd->ghostedcoloring;
269: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
270: }
271: ISColoringReference(*coloring);
272: return(0);
273: }
275: /* ---------------------------------------------------------------------------------*/
279: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
280: {
281: PetscErrorCode ierr;
282: 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;
283: PetscInt ncolors;
284: MPI_Comm comm;
285: DMBoundaryType bx,by,bz;
286: DMDAStencilType st;
287: ISColoringValue *colors;
288: DM_DA *dd = (DM_DA*)da->data;
291: /*
292: nc - number of components per grid point
293: col - number of colors needed in one direction for single component problem
295: */
296: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
297: col = 2*s + 1;
298: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
299: by 2*stencil_width + 1\n");
300: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
301: by 2*stencil_width + 1\n");
302: if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
303: by 2*stencil_width + 1\n");
305: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
306: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
307: PetscObjectGetComm((PetscObject)da,&comm);
309: /* create the coloring */
310: if (ctype == IS_COLORING_GLOBAL) {
311: if (!dd->localcoloring) {
312: PetscMalloc1(nc*nx*ny*nz,&colors);
313: ii = 0;
314: for (k=zs; k<zs+nz; k++) {
315: for (j=ys; j<ys+ny; j++) {
316: for (i=xs; i<xs+nx; i++) {
317: for (l=0; l<nc; l++) {
318: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
319: }
320: }
321: }
322: }
323: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
324: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
325: }
326: *coloring = dd->localcoloring;
327: } else if (ctype == IS_COLORING_GHOSTED) {
328: if (!dd->ghostedcoloring) {
329: PetscMalloc1(nc*gnx*gny*gnz,&colors);
330: ii = 0;
331: for (k=gzs; k<gzs+gnz; k++) {
332: for (j=gys; j<gys+gny; j++) {
333: for (i=gxs; i<gxs+gnx; i++) {
334: for (l=0; l<nc; l++) {
335: /* the complicated stuff is to handle periodic boundaries */
336: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
337: }
338: }
339: }
340: }
341: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
342: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
343: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
344: }
345: *coloring = dd->ghostedcoloring;
346: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
347: ISColoringReference(*coloring);
348: return(0);
349: }
351: /* ---------------------------------------------------------------------------------*/
355: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
356: {
357: PetscErrorCode ierr;
358: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
359: PetscInt ncolors;
360: MPI_Comm comm;
361: DMBoundaryType bx;
362: ISColoringValue *colors;
363: DM_DA *dd = (DM_DA*)da->data;
366: /*
367: nc - number of components per grid point
368: col - number of colors needed in one direction for single component problem
370: */
371: DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
372: col = 2*s + 1;
374: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\
375: by 2*stencil_width + 1 %d\n",(int)m,(int)col);
377: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
378: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
379: PetscObjectGetComm((PetscObject)da,&comm);
381: /* create the coloring */
382: if (ctype == IS_COLORING_GLOBAL) {
383: if (!dd->localcoloring) {
384: PetscMalloc1(nc*nx,&colors);
385: if (dd->ofillcols) {
386: PetscInt tc = 0;
387: for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
388: i1 = 0;
389: for (i=xs; i<xs+nx; i++) {
390: for (l=0; l<nc; l++) {
391: if (dd->ofillcols[l] && (i % col)) {
392: colors[i1++] = nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
393: } else {
394: colors[i1++] = l;
395: }
396: }
397: }
398: ncolors = nc + 2*s*tc;
399: } else {
400: i1 = 0;
401: for (i=xs; i<xs+nx; i++) {
402: for (l=0; l<nc; l++) {
403: colors[i1++] = l + nc*(i % col);
404: }
405: }
406: ncolors = nc + nc*(col-1);
407: }
408: ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
409: }
410: *coloring = dd->localcoloring;
411: } else if (ctype == IS_COLORING_GHOSTED) {
412: if (!dd->ghostedcoloring) {
413: PetscMalloc1(nc*gnx,&colors);
414: i1 = 0;
415: for (i=gxs; i<gxs+gnx; i++) {
416: for (l=0; l<nc; l++) {
417: /* the complicated stuff is to handle periodic boundaries */
418: colors[i1++] = l + nc*(SetInRange(i,m) % col);
419: }
420: }
421: ncolors = nc + nc*(col-1);
422: ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
423: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
424: }
425: *coloring = dd->ghostedcoloring;
426: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
427: ISColoringReference(*coloring);
428: return(0);
429: }
433: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
434: {
435: PetscErrorCode ierr;
436: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
437: PetscInt ncolors;
438: MPI_Comm comm;
439: DMBoundaryType bx,by;
440: ISColoringValue *colors;
441: DM_DA *dd = (DM_DA*)da->data;
444: /*
445: nc - number of components per grid point
446: col - number of colors needed in one direction for single component problem
448: */
449: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
450: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
451: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
452: PetscObjectGetComm((PetscObject)da,&comm);
454: if (bx == DM_BOUNDARY_PERIODIC && (m % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible by 5\n");
455: if (by == DM_BOUNDARY_PERIODIC && (n % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible by 5\n");
457: /* create the coloring */
458: if (ctype == IS_COLORING_GLOBAL) {
459: if (!dd->localcoloring) {
460: PetscMalloc1(nc*nx*ny,&colors);
461: ii = 0;
462: for (j=ys; j<ys+ny; j++) {
463: for (i=xs; i<xs+nx; i++) {
464: for (k=0; k<nc; k++) {
465: colors[ii++] = k + nc*((3*j+i) % 5);
466: }
467: }
468: }
469: ncolors = 5*nc;
470: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
471: }
472: *coloring = dd->localcoloring;
473: } else if (ctype == IS_COLORING_GHOSTED) {
474: if (!dd->ghostedcoloring) {
475: PetscMalloc1(nc*gnx*gny,&colors);
476: ii = 0;
477: for (j=gys; j<gys+gny; j++) {
478: for (i=gxs; i<gxs+gnx; i++) {
479: for (k=0; k<nc; k++) {
480: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
481: }
482: }
483: }
484: ncolors = 5*nc;
485: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
486: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_GHOSTED);
487: }
488: *coloring = dd->ghostedcoloring;
489: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
490: return(0);
491: }
493: /* =========================================================================== */
494: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
495: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
496: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
497: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
498: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
499: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
500: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
501: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
502: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
503: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
507: /*@C
508: MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix
510: Logically Collective on Mat
512: Input Parameters:
513: + mat - the matrix
514: - da - the da
516: Level: intermediate
518: @*/
519: PetscErrorCode MatSetupDM(Mat mat,DM da)
520: {
526: PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
527: return(0);
528: }
532: PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer)
533: {
534: DM da;
535: PetscErrorCode ierr;
536: const char *prefix;
537: Mat Anatural;
538: AO ao;
539: PetscInt rstart,rend,*petsc,i;
540: IS is;
541: MPI_Comm comm;
542: PetscViewerFormat format;
545: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
546: PetscViewerGetFormat(viewer,&format);
547: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
549: PetscObjectGetComm((PetscObject)A,&comm);
550: MatGetDM(A, &da);
551: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
553: DMDAGetAO(da,&ao);
554: MatGetOwnershipRange(A,&rstart,&rend);
555: PetscMalloc1(rend-rstart,&petsc);
556: for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
557: AOApplicationToPetsc(ao,rend-rstart,petsc);
558: ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);
560: /* call viewer on natural ordering */
561: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
562: ISDestroy(&is);
563: PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
564: PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
565: PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
566: (Anatural->ops->view)(Anatural,viewer);
567: MatDestroy(&Anatural);
568: return(0);
569: }
573: PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer)
574: {
575: DM da;
577: Mat Anatural,Aapp;
578: AO ao;
579: PetscInt rstart,rend,*app,i;
580: IS is;
581: MPI_Comm comm;
584: PetscObjectGetComm((PetscObject)A,&comm);
585: MatGetDM(A, &da);
586: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
588: /* Load the matrix in natural ordering */
589: MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
590: MatSetType(Anatural,((PetscObject)A)->type_name);
591: MatSetSizes(Anatural,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
592: MatLoad(Anatural,viewer);
594: /* Map natural ordering to application ordering and create IS */
595: DMDAGetAO(da,&ao);
596: MatGetOwnershipRange(Anatural,&rstart,&rend);
597: PetscMalloc1(rend-rstart,&app);
598: for (i=rstart; i<rend; i++) app[i-rstart] = i;
599: AOPetscToApplication(ao,rend-rstart,app);
600: ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);
602: /* Do permutation and replace header */
603: MatGetSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
604: MatHeaderReplace(A,&Aapp);
605: ISDestroy(&is);
606: MatDestroy(&Anatural);
607: return(0);
608: }
612: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
613: {
615: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
616: Mat A;
617: MPI_Comm comm;
618: MatType Atype;
619: PetscSection section, sectionGlobal;
620: void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL;
621: MatType mtype;
622: PetscMPIInt size;
623: DM_DA *dd = (DM_DA*)da->data;
626: MatInitializePackage();
627: mtype = da->mattype;
629: DMGetDefaultSection(da, §ion);
630: if (section) {
631: PetscInt bs = -1;
632: PetscInt localSize;
633: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
635: DMGetDefaultGlobalSection(da, §ionGlobal);
636: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
637: MatCreate(PetscObjectComm((PetscObject)da), J);
638: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
639: MatSetType(*J, mtype);
640: MatSetFromOptions(*J);
641: PetscStrcmp(mtype, MATSHELL, &isShell);
642: PetscStrcmp(mtype, MATBAIJ, &isBlock);
643: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
644: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
645: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
646: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
647: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
648: /* Check for symmetric storage */
649: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
650: if (isSymmetric) {
651: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
652: }
653: if (!isShell) {
654: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
656: if (bs < 0) {
657: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
658: PetscInt pStart, pEnd, p, dof;
660: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
661: for (p = pStart; p < pEnd; ++p) {
662: PetscSectionGetDof(sectionGlobal, p, &dof);
663: if (dof) {
664: bs = dof;
665: break;
666: }
667: }
668: } else {
669: bs = 1;
670: }
671: /* Must have same blocksize on all procs (some might have no points) */
672: bsLocal = bs;
673: MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
674: }
675: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
676: /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
677: PetscFree4(dnz, onz, dnzu, onzu);
678: }
679: }
680: /*
681: m
682: ------------------------------------------------------
683: | |
684: | |
685: | ---------------------- |
686: | | | |
687: n | ny | | |
688: | | | |
689: | .--------------------- |
690: | (xs,ys) nx |
691: | . |
692: | (gxs,gys) |
693: | |
694: -----------------------------------------------------
695: */
697: /*
698: nc - number of components per grid point
699: col - number of colors needed in one direction for single component problem
701: */
702: M = dd->M;
703: N = dd->N;
704: P = dd->P;
705: dim = da->dim;
706: dof = dd->w;
707: /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
708: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
709: PetscObjectGetComm((PetscObject)da,&comm);
710: MatCreate(comm,&A);
711: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
712: MatSetType(A,mtype);
713: MatSetDM(A,da);
714: MatSetFromOptions(A);
715: MatGetType(A,&Atype);
716: /*
717: We do not provide a getmatrix function in the DMDA operations because
718: the basic DMDA does not know about matrices. We think of DMDA as being more
719: more low-level than matrices. This is kind of cheating but, cause sometimes
720: we think of DMDA has higher level than matrices.
722: We could switch based on Atype (or mtype), but we do not since the
723: specialized setting routines depend only the particular preallocation
724: details of the matrix, not the type itself.
725: */
726: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
727: if (!aij) {
728: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
729: }
730: if (!aij) {
731: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
732: if (!baij) {
733: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
734: }
735: if (!baij) {
736: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
737: if (!sbaij) {
738: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
739: }
740: }
741: }
742: if (aij) {
743: if (dim == 1) {
744: if (dd->ofill) {
745: DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
746: } else {
747: DMCreateMatrix_DA_1d_MPIAIJ(da,A);
748: }
749: } else if (dim == 2) {
750: if (dd->ofill) {
751: DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
752: } else {
753: DMCreateMatrix_DA_2d_MPIAIJ(da,A);
754: }
755: } else if (dim == 3) {
756: if (dd->ofill) {
757: DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
758: } else {
759: DMCreateMatrix_DA_3d_MPIAIJ(da,A);
760: }
761: }
762: } else if (baij) {
763: if (dim == 2) {
764: DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
765: } else if (dim == 3) {
766: DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
767: } 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);
768: } else if (sbaij) {
769: if (dim == 2) {
770: DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
771: } else if (dim == 3) {
772: DMCreateMatrix_DA_3d_MPISBAIJ(da,A);
773: } 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);
774: } else {
775: ISLocalToGlobalMapping ltog;
776: DMGetLocalToGlobalMapping(da,<og);
777: MatSetUp(A);
778: MatSetLocalToGlobalMapping(A,ltog,ltog);
779: }
780: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
781: MatSetStencil(A,dim,dims,starts,dof);
782: MatSetDM(A,da);
783: MPI_Comm_size(comm,&size);
784: if (size > 1) {
785: /* change viewer to display matrix in natural ordering */
786: MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
787: MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
788: }
789: *J = A;
790: return(0);
791: }
793: /* ---------------------------------------------------------------------------------*/
796: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
797: {
798: PetscErrorCode ierr;
799: 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;
800: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
801: MPI_Comm comm;
802: PetscScalar *values;
803: DMBoundaryType bx,by;
804: ISLocalToGlobalMapping ltog;
805: DMDAStencilType st;
808: /*
809: nc - number of components per grid point
810: col - number of colors needed in one direction for single component problem
812: */
813: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
814: col = 2*s + 1;
815: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
816: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
817: PetscObjectGetComm((PetscObject)da,&comm);
819: PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
820: DMGetLocalToGlobalMapping(da,<og);
822: MatSetBlockSize(J,nc);
823: /* determine the matrix preallocation information */
824: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
825: for (i=xs; i<xs+nx; i++) {
827: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
828: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
830: for (j=ys; j<ys+ny; j++) {
831: slot = i - gxs + gnx*(j - gys);
833: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
834: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
836: cnt = 0;
837: for (k=0; k<nc; k++) {
838: for (l=lstart; l<lend+1; l++) {
839: for (p=pstart; p<pend+1; p++) {
840: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
841: cols[cnt++] = k + nc*(slot + gnx*l + p);
842: }
843: }
844: }
845: rows[k] = k + nc*(slot);
846: }
847: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
848: }
849: }
850: MatSetBlockSize(J,nc);
851: MatSeqAIJSetPreallocation(J,0,dnz);
852: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
853: MatPreallocateFinalize(dnz,onz);
855: MatSetLocalToGlobalMapping(J,ltog,ltog);
857: /*
858: For each node in the grid: we get the neighbors in the local (on processor ordering
859: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
860: PETSc ordering.
861: */
862: if (!da->prealloc_only) {
863: PetscCalloc1(col*col*nc*nc,&values);
864: for (i=xs; i<xs+nx; i++) {
866: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
867: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
869: for (j=ys; j<ys+ny; j++) {
870: slot = i - gxs + gnx*(j - gys);
872: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
873: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
875: cnt = 0;
876: for (k=0; k<nc; k++) {
877: for (l=lstart; l<lend+1; l++) {
878: for (p=pstart; p<pend+1; p++) {
879: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
880: cols[cnt++] = k + nc*(slot + gnx*l + p);
881: }
882: }
883: }
884: rows[k] = k + nc*(slot);
885: }
886: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
887: }
888: }
889: PetscFree(values);
890: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
891: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
892: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
893: }
894: PetscFree2(rows,cols);
895: return(0);
896: }
900: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
901: {
902: PetscErrorCode ierr;
903: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
904: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p;
905: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
906: DM_DA *dd = (DM_DA*)da->data;
907: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
908: MPI_Comm comm;
909: PetscScalar *values;
910: DMBoundaryType bx,by;
911: ISLocalToGlobalMapping ltog;
912: DMDAStencilType st;
915: /*
916: nc - number of components per grid point
917: col - number of colors needed in one direction for single component problem
919: */
920: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
921: col = 2*s + 1;
922: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
923: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
924: PetscObjectGetComm((PetscObject)da,&comm);
926: PetscMalloc1(col*col*nc,&cols);
927: DMGetLocalToGlobalMapping(da,<og);
929: MatSetBlockSize(J,nc);
930: /* determine the matrix preallocation information */
931: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
932: for (i=xs; i<xs+nx; i++) {
934: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
935: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
937: for (j=ys; j<ys+ny; j++) {
938: slot = i - gxs + gnx*(j - gys);
940: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
941: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
943: for (k=0; k<nc; k++) {
944: cnt = 0;
945: for (l=lstart; l<lend+1; l++) {
946: for (p=pstart; p<pend+1; p++) {
947: if (l || p) {
948: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
949: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
950: }
951: } else {
952: if (dfill) {
953: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
954: } else {
955: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
956: }
957: }
958: }
959: }
960: row = k + nc*(slot);
961: maxcnt = PetscMax(maxcnt,cnt);
962: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
963: }
964: }
965: }
966: MatSeqAIJSetPreallocation(J,0,dnz);
967: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
968: MatPreallocateFinalize(dnz,onz);
969: MatSetLocalToGlobalMapping(J,ltog,ltog);
971: /*
972: For each node in the grid: we get the neighbors in the local (on processor ordering
973: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
974: PETSc ordering.
975: */
976: if (!da->prealloc_only) {
977: PetscCalloc1(maxcnt,&values);
978: for (i=xs; i<xs+nx; i++) {
980: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
981: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
983: for (j=ys; j<ys+ny; j++) {
984: slot = i - gxs + gnx*(j - gys);
986: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
987: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
989: for (k=0; k<nc; k++) {
990: cnt = 0;
991: for (l=lstart; l<lend+1; l++) {
992: for (p=pstart; p<pend+1; p++) {
993: if (l || p) {
994: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
995: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
996: }
997: } else {
998: if (dfill) {
999: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1000: } else {
1001: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1002: }
1003: }
1004: }
1005: }
1006: row = k + nc*(slot);
1007: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1008: }
1009: }
1010: }
1011: PetscFree(values);
1012: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1013: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1014: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1015: }
1016: PetscFree(cols);
1017: return(0);
1018: }
1020: /* ---------------------------------------------------------------------------------*/
1024: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1025: {
1026: PetscErrorCode ierr;
1027: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1028: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1029: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1030: MPI_Comm comm;
1031: PetscScalar *values;
1032: DMBoundaryType bx,by,bz;
1033: ISLocalToGlobalMapping ltog;
1034: DMDAStencilType st;
1037: /*
1038: nc - number of components per grid point
1039: col - number of colors needed in one direction for single component problem
1041: */
1042: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1043: col = 2*s + 1;
1045: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1046: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1047: PetscObjectGetComm((PetscObject)da,&comm);
1049: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1050: DMGetLocalToGlobalMapping(da,<og);
1052: MatSetBlockSize(J,nc);
1053: /* determine the matrix preallocation information */
1054: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1055: for (i=xs; i<xs+nx; i++) {
1056: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1057: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1058: for (j=ys; j<ys+ny; j++) {
1059: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1060: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1061: for (k=zs; k<zs+nz; k++) {
1062: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1063: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1065: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1067: cnt = 0;
1068: for (l=0; l<nc; l++) {
1069: for (ii=istart; ii<iend+1; ii++) {
1070: for (jj=jstart; jj<jend+1; jj++) {
1071: for (kk=kstart; kk<kend+1; kk++) {
1072: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1073: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1074: }
1075: }
1076: }
1077: }
1078: rows[l] = l + nc*(slot);
1079: }
1080: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1081: }
1082: }
1083: }
1084: MatSetBlockSize(J,nc);
1085: MatSeqAIJSetPreallocation(J,0,dnz);
1086: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1087: MatPreallocateFinalize(dnz,onz);
1088: MatSetLocalToGlobalMapping(J,ltog,ltog);
1090: /*
1091: For each node in the grid: we get the neighbors in the local (on processor ordering
1092: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1093: PETSc ordering.
1094: */
1095: if (!da->prealloc_only) {
1096: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1097: for (i=xs; i<xs+nx; i++) {
1098: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1099: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1100: for (j=ys; j<ys+ny; j++) {
1101: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1102: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1103: for (k=zs; k<zs+nz; k++) {
1104: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1105: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1107: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1109: cnt = 0;
1110: for (l=0; l<nc; l++) {
1111: for (ii=istart; ii<iend+1; ii++) {
1112: for (jj=jstart; jj<jend+1; jj++) {
1113: for (kk=kstart; kk<kend+1; kk++) {
1114: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1115: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1116: }
1117: }
1118: }
1119: }
1120: rows[l] = l + nc*(slot);
1121: }
1122: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1123: }
1124: }
1125: }
1126: PetscFree(values);
1127: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1128: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1129: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1130: }
1131: PetscFree2(rows,cols);
1132: return(0);
1133: }
1135: /* ---------------------------------------------------------------------------------*/
1139: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1140: {
1141: PetscErrorCode ierr;
1142: DM_DA *dd = (DM_DA*)da->data;
1143: PetscInt xs,nx,i,j,gxs,gnx,row,k,l;
1144: PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1145: PetscInt *ofill = dd->ofill,*dfill = dd->dfill;
1146: PetscScalar *values;
1147: DMBoundaryType bx;
1148: ISLocalToGlobalMapping ltog;
1149: PetscMPIInt rank,size;
1152: if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1153: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1154: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1156: /*
1157: nc - number of components per grid point
1159: */
1160: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1161: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1162: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1164: MatSetBlockSize(J,nc);
1165: PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);
1167: if (nx < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Need at least two grid points per process");
1168: /*
1169: note should be smaller for first and last process with no periodic
1170: does not handle dfill
1171: */
1172: cnt = 0;
1173: /* coupling with process to the left */
1174: for (i=0; i<s; i++) {
1175: for (j=0; j<nc; j++) {
1176: ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1177: cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1178: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1179: cnt++;
1180: }
1181: }
1182: for (i=s; i<nx-s; i++) {
1183: for (j=0; j<nc; j++) {
1184: cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1185: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1186: cnt++;
1187: }
1188: }
1189: /* coupling with process to the right */
1190: for (i=nx-s; i<nx; i++) {
1191: for (j=0; j<nc; j++) {
1192: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1193: cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1194: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1195: cnt++;
1196: }
1197: }
1199: MatSeqAIJSetPreallocation(J,0,cols);
1200: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1201: PetscFree2(cols,ocols);
1203: DMGetLocalToGlobalMapping(da,<og);
1204: MatSetLocalToGlobalMapping(J,ltog,ltog);
1206: /*
1207: For each node in the grid: we get the neighbors in the local (on processor ordering
1208: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1209: PETSc ordering.
1210: */
1211: if (!da->prealloc_only) {
1212: PetscCalloc2(maxcnt,&values,maxcnt,&cols);
1214: row = xs*nc;
1215: /* coupling with process to the left */
1216: for (i=xs; i<xs+s; i++) {
1217: for (j=0; j<nc; j++) {
1218: cnt = 0;
1219: if (rank) {
1220: for (l=0; l<s; l++) {
1221: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1222: }
1223: }
1224: if (dfill) {
1225: for (k=dfill[j]; k<dfill[j+1]; k++) {
1226: cols[cnt++] = i*nc + dfill[k];
1227: }
1228: } else {
1229: for (k=0; k<nc; k++) {
1230: cols[cnt++] = i*nc + k;
1231: }
1232: }
1233: for (l=0; l<s; l++) {
1234: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1235: }
1236: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1237: row++;
1238: }
1239: }
1240: for (i=xs+s; i<xs+nx-s; i++) {
1241: for (j=0; j<nc; j++) {
1242: cnt = 0;
1243: for (l=0; l<s; l++) {
1244: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1245: }
1246: if (dfill) {
1247: for (k=dfill[j]; k<dfill[j+1]; k++) {
1248: cols[cnt++] = i*nc + dfill[k];
1249: }
1250: } else {
1251: for (k=0; k<nc; k++) {
1252: cols[cnt++] = i*nc + k;
1253: }
1254: }
1255: for (l=0; l<s; l++) {
1256: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1257: }
1258: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1259: row++;
1260: }
1261: }
1262: /* coupling with process to the right */
1263: for (i=xs+nx-s; i<xs+nx; i++) {
1264: for (j=0; j<nc; j++) {
1265: cnt = 0;
1266: for (l=0; l<s; l++) {
1267: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1268: }
1269: if (dfill) {
1270: for (k=dfill[j]; k<dfill[j+1]; k++) {
1271: cols[cnt++] = i*nc + dfill[k];
1272: }
1273: } else {
1274: for (k=0; k<nc; k++) {
1275: cols[cnt++] = i*nc + k;
1276: }
1277: }
1278: if (rank < size-1) {
1279: for (l=0; l<s; l++) {
1280: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1281: }
1282: }
1283: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1284: row++;
1285: }
1286: }
1287: PetscFree2(values,cols);
1288: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1289: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1290: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1291: }
1292: return(0);
1293: }
1295: /* ---------------------------------------------------------------------------------*/
1299: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1300: {
1301: PetscErrorCode ierr;
1302: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1303: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1304: PetscInt istart,iend;
1305: PetscScalar *values;
1306: DMBoundaryType bx;
1307: ISLocalToGlobalMapping ltog;
1310: /*
1311: nc - number of components per grid point
1312: col - number of colors needed in one direction for single component problem
1314: */
1315: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1316: col = 2*s + 1;
1318: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1319: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1321: MatSetBlockSize(J,nc);
1322: MatSeqAIJSetPreallocation(J,col*nc,0);
1323: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1325: DMGetLocalToGlobalMapping(da,<og);
1326: MatSetLocalToGlobalMapping(J,ltog,ltog);
1328: /*
1329: For each node in the grid: we get the neighbors in the local (on processor ordering
1330: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1331: PETSc ordering.
1332: */
1333: if (!da->prealloc_only) {
1334: PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1335: PetscCalloc1(col*nc*nc,&values);
1336: for (i=xs; i<xs+nx; i++) {
1337: istart = PetscMax(-s,gxs - i);
1338: iend = PetscMin(s,gxs + gnx - i - 1);
1339: slot = i - gxs;
1341: cnt = 0;
1342: for (l=0; l<nc; l++) {
1343: for (i1=istart; i1<iend+1; i1++) {
1344: cols[cnt++] = l + nc*(slot + i1);
1345: }
1346: rows[l] = l + nc*(slot);
1347: }
1348: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1349: }
1350: PetscFree(values);
1351: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1352: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1353: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1354: PetscFree2(rows,cols);
1355: }
1356: return(0);
1357: }
1361: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1362: {
1363: PetscErrorCode ierr;
1364: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1365: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1366: PetscInt istart,iend,jstart,jend,ii,jj;
1367: MPI_Comm comm;
1368: PetscScalar *values;
1369: DMBoundaryType bx,by;
1370: DMDAStencilType st;
1371: ISLocalToGlobalMapping ltog;
1374: /*
1375: nc - number of components per grid point
1376: col - number of colors needed in one direction for single component problem
1377: */
1378: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1379: col = 2*s + 1;
1381: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1382: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1383: PetscObjectGetComm((PetscObject)da,&comm);
1385: PetscMalloc1(col*col*nc*nc,&cols);
1387: DMGetLocalToGlobalMapping(da,<og);
1389: /* determine the matrix preallocation information */
1390: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1391: for (i=xs; i<xs+nx; i++) {
1392: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1393: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1394: for (j=ys; j<ys+ny; j++) {
1395: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1396: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1397: slot = i - gxs + gnx*(j - gys);
1399: /* Find block columns in block row */
1400: cnt = 0;
1401: for (ii=istart; ii<iend+1; ii++) {
1402: for (jj=jstart; jj<jend+1; jj++) {
1403: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1404: cols[cnt++] = slot + ii + gnx*jj;
1405: }
1406: }
1407: }
1408: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1409: }
1410: }
1411: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1412: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1413: MatPreallocateFinalize(dnz,onz);
1415: MatSetLocalToGlobalMapping(J,ltog,ltog);
1417: /*
1418: For each node in the grid: we get the neighbors in the local (on processor ordering
1419: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1420: PETSc ordering.
1421: */
1422: if (!da->prealloc_only) {
1423: PetscCalloc1(col*col*nc*nc,&values);
1424: for (i=xs; i<xs+nx; i++) {
1425: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1426: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1427: for (j=ys; j<ys+ny; j++) {
1428: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1429: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1430: slot = i - gxs + gnx*(j - gys);
1431: cnt = 0;
1432: for (ii=istart; ii<iend+1; ii++) {
1433: for (jj=jstart; jj<jend+1; jj++) {
1434: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1435: cols[cnt++] = slot + ii + gnx*jj;
1436: }
1437: }
1438: }
1439: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1440: }
1441: }
1442: PetscFree(values);
1443: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1444: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1445: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1446: }
1447: PetscFree(cols);
1448: return(0);
1449: }
1453: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1454: {
1455: PetscErrorCode ierr;
1456: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1457: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1458: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1459: MPI_Comm comm;
1460: PetscScalar *values;
1461: DMBoundaryType bx,by,bz;
1462: DMDAStencilType st;
1463: ISLocalToGlobalMapping ltog;
1466: /*
1467: nc - number of components per grid point
1468: col - number of colors needed in one direction for single component problem
1470: */
1471: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1472: col = 2*s + 1;
1474: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1475: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1476: PetscObjectGetComm((PetscObject)da,&comm);
1478: PetscMalloc1(col*col*col,&cols);
1480: DMGetLocalToGlobalMapping(da,<og);
1482: /* determine the matrix preallocation information */
1483: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1484: for (i=xs; i<xs+nx; i++) {
1485: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1486: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1487: for (j=ys; j<ys+ny; j++) {
1488: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1489: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1490: for (k=zs; k<zs+nz; k++) {
1491: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1492: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1494: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1496: /* Find block columns in block row */
1497: cnt = 0;
1498: for (ii=istart; ii<iend+1; ii++) {
1499: for (jj=jstart; jj<jend+1; jj++) {
1500: for (kk=kstart; kk<kend+1; kk++) {
1501: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1502: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1503: }
1504: }
1505: }
1506: }
1507: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1508: }
1509: }
1510: }
1511: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1512: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1513: MatPreallocateFinalize(dnz,onz);
1515: MatSetLocalToGlobalMapping(J,ltog,ltog);
1517: /*
1518: For each node in the grid: we get the neighbors in the local (on processor ordering
1519: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1520: PETSc ordering.
1521: */
1522: if (!da->prealloc_only) {
1523: PetscCalloc1(col*col*col*nc*nc,&values);
1524: for (i=xs; i<xs+nx; i++) {
1525: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1526: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1527: for (j=ys; j<ys+ny; j++) {
1528: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1529: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1530: for (k=zs; k<zs+nz; k++) {
1531: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1532: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1534: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1536: cnt = 0;
1537: for (ii=istart; ii<iend+1; ii++) {
1538: for (jj=jstart; jj<jend+1; jj++) {
1539: for (kk=kstart; kk<kend+1; kk++) {
1540: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1541: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1542: }
1543: }
1544: }
1545: }
1546: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1547: }
1548: }
1549: }
1550: PetscFree(values);
1551: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1552: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1553: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1554: }
1555: PetscFree(cols);
1556: return(0);
1557: }
1561: /*
1562: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1563: identify in the local ordering with periodic domain.
1564: */
1565: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1566: {
1568: PetscInt i,n;
1571: ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1572: ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1573: for (i=0,n=0; i<*cnt; i++) {
1574: if (col[i] >= *row) col[n++] = col[i];
1575: }
1576: *cnt = n;
1577: return(0);
1578: }
1582: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1583: {
1584: PetscErrorCode ierr;
1585: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1586: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1587: PetscInt istart,iend,jstart,jend,ii,jj;
1588: MPI_Comm comm;
1589: PetscScalar *values;
1590: DMBoundaryType bx,by;
1591: DMDAStencilType st;
1592: ISLocalToGlobalMapping ltog;
1595: /*
1596: nc - number of components per grid point
1597: col - number of colors needed in one direction for single component problem
1598: */
1599: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1600: col = 2*s + 1;
1602: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1603: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1604: PetscObjectGetComm((PetscObject)da,&comm);
1606: PetscMalloc1(col*col*nc*nc,&cols);
1608: DMGetLocalToGlobalMapping(da,<og);
1610: /* determine the matrix preallocation information */
1611: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1612: for (i=xs; i<xs+nx; i++) {
1613: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1614: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1615: for (j=ys; j<ys+ny; j++) {
1616: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1617: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1618: slot = i - gxs + gnx*(j - gys);
1620: /* Find block columns in block row */
1621: cnt = 0;
1622: for (ii=istart; ii<iend+1; ii++) {
1623: for (jj=jstart; jj<jend+1; jj++) {
1624: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1625: cols[cnt++] = slot + ii + gnx*jj;
1626: }
1627: }
1628: }
1629: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1630: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1631: }
1632: }
1633: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1634: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1635: MatPreallocateFinalize(dnz,onz);
1637: MatSetLocalToGlobalMapping(J,ltog,ltog);
1639: /*
1640: For each node in the grid: we get the neighbors in the local (on processor ordering
1641: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1642: PETSc ordering.
1643: */
1644: if (!da->prealloc_only) {
1645: PetscCalloc1(col*col*nc*nc,&values);
1646: for (i=xs; i<xs+nx; i++) {
1647: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1648: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1649: for (j=ys; j<ys+ny; j++) {
1650: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1651: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1652: slot = i - gxs + gnx*(j - gys);
1654: /* Find block columns in block row */
1655: cnt = 0;
1656: for (ii=istart; ii<iend+1; ii++) {
1657: for (jj=jstart; jj<jend+1; jj++) {
1658: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1659: cols[cnt++] = slot + ii + gnx*jj;
1660: }
1661: }
1662: }
1663: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1664: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1665: }
1666: }
1667: PetscFree(values);
1668: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1669: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1670: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1671: }
1672: PetscFree(cols);
1673: return(0);
1674: }
1678: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1679: {
1680: PetscErrorCode ierr;
1681: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1682: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1683: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1684: MPI_Comm comm;
1685: PetscScalar *values;
1686: DMBoundaryType bx,by,bz;
1687: DMDAStencilType st;
1688: ISLocalToGlobalMapping ltog;
1691: /*
1692: nc - number of components per grid point
1693: col - number of colors needed in one direction for single component problem
1694: */
1695: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1696: col = 2*s + 1;
1698: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1699: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1700: PetscObjectGetComm((PetscObject)da,&comm);
1702: /* create the matrix */
1703: PetscMalloc1(col*col*col,&cols);
1705: DMGetLocalToGlobalMapping(da,<og);
1707: /* determine the matrix preallocation information */
1708: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1709: for (i=xs; i<xs+nx; i++) {
1710: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1711: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1712: for (j=ys; j<ys+ny; j++) {
1713: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1714: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1715: for (k=zs; k<zs+nz; k++) {
1716: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1717: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1719: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1721: /* Find block columns in block row */
1722: cnt = 0;
1723: for (ii=istart; ii<iend+1; ii++) {
1724: for (jj=jstart; jj<jend+1; jj++) {
1725: for (kk=kstart; kk<kend+1; kk++) {
1726: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1727: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1728: }
1729: }
1730: }
1731: }
1732: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1733: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1734: }
1735: }
1736: }
1737: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1738: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1739: MatPreallocateFinalize(dnz,onz);
1741: MatSetLocalToGlobalMapping(J,ltog,ltog);
1743: /*
1744: For each node in the grid: we get the neighbors in the local (on processor ordering
1745: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1746: PETSc ordering.
1747: */
1748: if (!da->prealloc_only) {
1749: PetscCalloc1(col*col*col*nc*nc,&values);
1750: for (i=xs; i<xs+nx; i++) {
1751: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1752: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1753: for (j=ys; j<ys+ny; j++) {
1754: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1755: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1756: for (k=zs; k<zs+nz; k++) {
1757: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1758: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1760: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1762: cnt = 0;
1763: for (ii=istart; ii<iend+1; ii++) {
1764: for (jj=jstart; jj<jend+1; jj++) {
1765: for (kk=kstart; kk<kend+1; kk++) {
1766: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1767: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1768: }
1769: }
1770: }
1771: }
1772: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1773: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1774: }
1775: }
1776: }
1777: PetscFree(values);
1778: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1779: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1780: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1781: }
1782: PetscFree(cols);
1783: return(0);
1784: }
1786: /* ---------------------------------------------------------------------------------*/
1790: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
1791: {
1792: PetscErrorCode ierr;
1793: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1794: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
1795: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1796: DM_DA *dd = (DM_DA*)da->data;
1797: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
1798: MPI_Comm comm;
1799: PetscScalar *values;
1800: DMBoundaryType bx,by,bz;
1801: ISLocalToGlobalMapping ltog;
1802: DMDAStencilType st;
1805: /*
1806: nc - number of components per grid point
1807: col - number of colors needed in one direction for single component problem
1809: */
1810: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1811: col = 2*s + 1;
1812: if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1813: by 2*stencil_width + 1\n");
1814: if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1815: by 2*stencil_width + 1\n");
1816: if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1817: by 2*stencil_width + 1\n");
1819: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1820: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1821: PetscObjectGetComm((PetscObject)da,&comm);
1823: PetscMalloc1(col*col*col*nc,&cols);
1824: DMGetLocalToGlobalMapping(da,<og);
1826: /* determine the matrix preallocation information */
1827: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1829: MatSetBlockSize(J,nc);
1830: for (i=xs; i<xs+nx; i++) {
1831: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1832: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1833: for (j=ys; j<ys+ny; j++) {
1834: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1835: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1836: for (k=zs; k<zs+nz; k++) {
1837: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1838: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1840: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1842: for (l=0; l<nc; l++) {
1843: cnt = 0;
1844: for (ii=istart; ii<iend+1; ii++) {
1845: for (jj=jstart; jj<jend+1; jj++) {
1846: for (kk=kstart; kk<kend+1; kk++) {
1847: if (ii || jj || kk) {
1848: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1849: 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);
1850: }
1851: } else {
1852: if (dfill) {
1853: 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);
1854: } else {
1855: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1856: }
1857: }
1858: }
1859: }
1860: }
1861: row = l + nc*(slot);
1862: maxcnt = PetscMax(maxcnt,cnt);
1863: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1864: }
1865: }
1866: }
1867: }
1868: MatSeqAIJSetPreallocation(J,0,dnz);
1869: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1870: MatPreallocateFinalize(dnz,onz);
1871: MatSetLocalToGlobalMapping(J,ltog,ltog);
1873: /*
1874: For each node in the grid: we get the neighbors in the local (on processor ordering
1875: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1876: PETSc ordering.
1877: */
1878: if (!da->prealloc_only) {
1879: PetscCalloc1(maxcnt,&values);
1880: for (i=xs; i<xs+nx; i++) {
1881: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1882: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1883: for (j=ys; j<ys+ny; j++) {
1884: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1885: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1886: for (k=zs; k<zs+nz; k++) {
1887: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1888: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1890: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1892: for (l=0; l<nc; l++) {
1893: cnt = 0;
1894: for (ii=istart; ii<iend+1; ii++) {
1895: for (jj=jstart; jj<jend+1; jj++) {
1896: for (kk=kstart; kk<kend+1; kk++) {
1897: if (ii || jj || kk) {
1898: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1899: 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);
1900: }
1901: } else {
1902: if (dfill) {
1903: 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);
1904: } else {
1905: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1906: }
1907: }
1908: }
1909: }
1910: }
1911: row = l + nc*(slot);
1912: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1913: }
1914: }
1915: }
1916: }
1917: PetscFree(values);
1918: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1919: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1920: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1921: }
1922: PetscFree(cols);
1923: return(0);
1924: }