Actual source code: fdda.c
petsc-3.9.4 2018-09-11
2: #include <petsc/private/dmdaimpl.h>
3: #include <petscmat.h>
5: extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring*);
6: extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring*);
7: extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring*);
8: extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring*);
10: /*
11: For ghost i that may be negative or greater than the upper bound this
12: maps it into the 0:m-1 range using periodicity
13: */
14: #define SetInRange(i,m) ((i < 0) ? m+i : ((i >= m) ? i-m : i))
16: static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill)
17: {
19: PetscInt i,j,nz,*fill;
22: if (!dfill) return(0);
24: /* count number nonzeros */
25: nz = 0;
26: for (i=0; i<w; i++) {
27: for (j=0; j<w; j++) {
28: if (dfill[w*i+j]) nz++;
29: }
30: }
31: PetscMalloc1(nz + w + 1,&fill);
32: /* construct modified CSR storage of nonzero structure */
33: /* fill[0 -- w] marks starts of each row of column indices (and end of last row)
34: so fill[1] - fill[0] gives number of nonzeros in first row etc */
35: nz = w + 1;
36: for (i=0; i<w; i++) {
37: fill[i] = nz;
38: for (j=0; j<w; j++) {
39: if (dfill[w*i+j]) {
40: fill[nz] = j;
41: nz++;
42: }
43: }
44: }
45: fill[w] = nz;
47: *rfill = fill;
48: return(0);
49: }
51: /*@
52: DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
53: of the matrix returned by DMCreateMatrix().
55: Logically Collective on DMDA
57: Input Parameter:
58: + da - the distributed array
59: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
60: - ofill - the fill pattern in the off-diagonal blocks
63: Level: developer
65: Notes: This only makes sense when you are doing multicomponent problems but using the
66: MPIAIJ matrix format
68: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
69: representing coupling and 0 entries for missing coupling. For example
70: $ dfill[9] = {1, 0, 0,
71: $ 1, 1, 0,
72: $ 0, 1, 1}
73: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
74: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
75: diagonal block).
77: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
78: can be represented in the dfill, ofill format
80: Contributed by Glenn Hammond
82: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
84: @*/
85: PetscErrorCode DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill)
86: {
87: DM_DA *dd = (DM_DA*)da->data;
89: PetscInt i,k,cnt = 1;
92: DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);
93: DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);
95: /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of
96: columns to the left with any nonzeros in them plus 1 */
97: PetscCalloc1(dd->w,&dd->ofillcols);
98: for (i=0; i<dd->w; i++) {
99: for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
100: }
101: for (i=0; i<dd->w; i++) {
102: if (dd->ofillcols[i]) {
103: dd->ofillcols[i] = cnt++;
104: }
105: }
106: return(0);
107: }
110: PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring)
111: {
112: PetscErrorCode ierr;
113: PetscInt dim,m,n,p,nc;
114: DMBoundaryType bx,by,bz;
115: MPI_Comm comm;
116: PetscMPIInt size;
117: PetscBool isBAIJ;
118: DM_DA *dd = (DM_DA*)da->data;
121: /*
122: m
123: ------------------------------------------------------
124: | |
125: | |
126: | ---------------------- |
127: | | | |
128: n | yn | | |
129: | | | |
130: | .--------------------- |
131: | (xs,ys) xn |
132: | . |
133: | (gxs,gys) |
134: | |
135: -----------------------------------------------------
136: */
138: /*
139: nc - number of components per grid point
140: col - number of colors needed in one direction for single component problem
142: */
143: DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);
145: PetscObjectGetComm((PetscObject)da,&comm);
146: MPI_Comm_size(comm,&size);
147: if (ctype == IS_COLORING_LOCAL) {
148: if (size == 1) {
149: ctype = IS_COLORING_GLOBAL;
150: } else if (dim > 1) {
151: if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) {
152: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_LOCAL cannot be used for periodic boundary condition having both ends of the domain on the same process");
153: }
154: }
155: }
157: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
158: matrices is for the blocks, not the individual matrix elements */
159: PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);
160: if (!isBAIJ) {PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);}
161: if (!isBAIJ) {PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);}
162: if (isBAIJ) {
163: dd->w = 1;
164: dd->xs = dd->xs/nc;
165: dd->xe = dd->xe/nc;
166: dd->Xs = dd->Xs/nc;
167: dd->Xe = dd->Xe/nc;
168: }
170: /*
171: We do not provide a getcoloring function in the DMDA operations because
172: the basic DMDA does not know about matrices. We think of DMDA as being more
173: more low-level then matrices.
174: */
175: if (dim == 1) {
176: DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);
177: } else if (dim == 2) {
178: DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);
179: } else if (dim == 3) {
180: DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);
181: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
182: if (isBAIJ) {
183: dd->w = nc;
184: dd->xs = dd->xs*nc;
185: dd->xe = dd->xe*nc;
186: dd->Xs = dd->Xs*nc;
187: dd->Xe = dd->Xe*nc;
188: }
189: return(0);
190: }
192: /* ---------------------------------------------------------------------------------*/
194: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
195: {
196: PetscErrorCode ierr;
197: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
198: PetscInt ncolors;
199: MPI_Comm comm;
200: DMBoundaryType bx,by;
201: DMDAStencilType st;
202: ISColoringValue *colors;
203: DM_DA *dd = (DM_DA*)da->data;
206: /*
207: nc - number of components per grid point
208: col - number of colors needed in one direction for single component problem
210: */
211: DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);
212: col = 2*s + 1;
213: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
214: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
215: PetscObjectGetComm((PetscObject)da,&comm);
217: /* special case as taught to us by Paul Hovland */
218: if (st == DMDA_STENCIL_STAR && s == 1) {
219: DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);
220: } else {
222: 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\
223: by 2*stencil_width + 1 (%d)\n", m, col);
224: 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\
225: by 2*stencil_width + 1 (%d)\n", n, col);
226: if (ctype == IS_COLORING_GLOBAL) {
227: if (!dd->localcoloring) {
228: PetscMalloc1(nc*nx*ny,&colors);
229: ii = 0;
230: for (j=ys; j<ys+ny; j++) {
231: for (i=xs; i<xs+nx; i++) {
232: for (k=0; k<nc; k++) {
233: colors[ii++] = k + nc*((i % col) + col*(j % col));
234: }
235: }
236: }
237: ncolors = nc + nc*(col-1 + col*(col-1));
238: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
239: }
240: *coloring = dd->localcoloring;
241: } else if (ctype == IS_COLORING_LOCAL) {
242: if (!dd->ghostedcoloring) {
243: PetscMalloc1(nc*gnx*gny,&colors);
244: ii = 0;
245: for (j=gys; j<gys+gny; j++) {
246: for (i=gxs; i<gxs+gnx; i++) {
247: for (k=0; k<nc; k++) {
248: /* the complicated stuff is to handle periodic boundaries */
249: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
250: }
251: }
252: }
253: ncolors = nc + nc*(col - 1 + col*(col-1));
254: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
255: /* PetscIntView(ncolors,(PetscInt*)colors,0); */
257: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
258: }
259: *coloring = dd->ghostedcoloring;
260: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
261: }
262: ISColoringReference(*coloring);
263: return(0);
264: }
266: /* ---------------------------------------------------------------------------------*/
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: DMBoundaryType 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
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 == 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\
288: by 2*stencil_width + 1\n");
289: 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\
290: by 2*stencil_width + 1\n");
291: 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\
292: by 2*stencil_width + 1\n");
294: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
295: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
296: PetscObjectGetComm((PetscObject)da,&comm);
298: /* create the coloring */
299: if (ctype == IS_COLORING_GLOBAL) {
300: if (!dd->localcoloring) {
301: PetscMalloc1(nc*nx*ny*nz,&colors);
302: ii = 0;
303: for (k=zs; k<zs+nz; k++) {
304: for (j=ys; j<ys+ny; j++) {
305: for (i=xs; i<xs+nx; i++) {
306: for (l=0; l<nc; l++) {
307: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
308: }
309: }
310: }
311: }
312: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
313: ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);
314: }
315: *coloring = dd->localcoloring;
316: } else if (ctype == IS_COLORING_LOCAL) {
317: if (!dd->ghostedcoloring) {
318: PetscMalloc1(nc*gnx*gny*gnz,&colors);
319: ii = 0;
320: for (k=gzs; k<gzs+gnz; k++) {
321: for (j=gys; j<gys+gny; j++) {
322: for (i=gxs; i<gxs+gnx; i++) {
323: for (l=0; l<nc; l++) {
324: /* the complicated stuff is to handle periodic boundaries */
325: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
326: }
327: }
328: }
329: }
330: ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1));
331: ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
332: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
333: }
334: *coloring = dd->ghostedcoloring;
335: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
336: ISColoringReference(*coloring);
337: return(0);
338: }
340: /* ---------------------------------------------------------------------------------*/
342: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
343: {
344: PetscErrorCode ierr;
345: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
346: PetscInt ncolors;
347: MPI_Comm comm;
348: DMBoundaryType bx;
349: ISColoringValue *colors;
350: DM_DA *dd = (DM_DA*)da->data;
353: /*
354: nc - number of components per grid point
355: col - number of colors needed in one direction for single component problem
357: */
358: DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);
359: col = 2*s + 1;
361: 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\
362: by 2*stencil_width + 1 %d\n",(int)m,(int)col);
364: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
365: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
366: PetscObjectGetComm((PetscObject)da,&comm);
368: /* create the coloring */
369: if (ctype == IS_COLORING_GLOBAL) {
370: if (!dd->localcoloring) {
371: PetscMalloc1(nc*nx,&colors);
372: if (dd->ofillcols) {
373: PetscInt tc = 0;
374: for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0);
375: i1 = 0;
376: for (i=xs; i<xs+nx; i++) {
377: for (l=0; l<nc; l++) {
378: if (dd->ofillcols[l] && (i % col)) {
379: colors[i1++] = nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l];
380: } else {
381: colors[i1++] = l;
382: }
383: }
384: }
385: ncolors = nc + 2*s*tc;
386: } else {
387: i1 = 0;
388: for (i=xs; i<xs+nx; i++) {
389: for (l=0; l<nc; l++) {
390: colors[i1++] = l + nc*(i % col);
391: }
392: }
393: ncolors = nc + nc*(col-1);
394: }
395: ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);
396: }
397: *coloring = dd->localcoloring;
398: } else if (ctype == IS_COLORING_LOCAL) {
399: if (!dd->ghostedcoloring) {
400: PetscMalloc1(nc*gnx,&colors);
401: i1 = 0;
402: for (i=gxs; i<gxs+gnx; i++) {
403: for (l=0; l<nc; l++) {
404: /* the complicated stuff is to handle periodic boundaries */
405: colors[i1++] = l + nc*(SetInRange(i,m) % col);
406: }
407: }
408: ncolors = nc + nc*(col-1);
409: ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
410: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
411: }
412: *coloring = dd->ghostedcoloring;
413: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
414: ISColoringReference(*coloring);
415: return(0);
416: }
418: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring)
419: {
420: PetscErrorCode ierr;
421: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
422: PetscInt ncolors;
423: MPI_Comm comm;
424: DMBoundaryType bx,by;
425: ISColoringValue *colors;
426: DM_DA *dd = (DM_DA*)da->data;
429: /*
430: nc - number of components per grid point
431: col - number of colors needed in one direction for single component problem
433: */
434: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);
435: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
436: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
437: PetscObjectGetComm((PetscObject)da,&comm);
439: 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");
440: 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");
442: /* create the coloring */
443: if (ctype == IS_COLORING_GLOBAL) {
444: if (!dd->localcoloring) {
445: PetscMalloc1(nc*nx*ny,&colors);
446: ii = 0;
447: for (j=ys; j<ys+ny; j++) {
448: for (i=xs; i<xs+nx; i++) {
449: for (k=0; k<nc; k++) {
450: colors[ii++] = k + nc*((3*j+i) % 5);
451: }
452: }
453: }
454: ncolors = 5*nc;
455: ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);
456: }
457: *coloring = dd->localcoloring;
458: } else if (ctype == IS_COLORING_LOCAL) {
459: if (!dd->ghostedcoloring) {
460: PetscMalloc1(nc*gnx*gny,&colors);
461: ii = 0;
462: for (j=gys; j<gys+gny; j++) {
463: for (i=gxs; i<gxs+gnx; i++) {
464: for (k=0; k<nc; k++) {
465: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
466: }
467: }
468: }
469: ncolors = 5*nc;
470: ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);
471: ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);
472: }
473: *coloring = dd->ghostedcoloring;
474: } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
475: return(0);
476: }
478: /* =========================================================================== */
479: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat);
480: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat);
481: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat);
482: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat);
483: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat);
484: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat);
485: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat);
486: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat);
487: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat);
488: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat);
489: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM,Mat);
490: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM,Mat);
491: extern PetscErrorCode DMCreateMatrix_DA_IS(DM,Mat);
493: /*@C
494: MatSetupDM - 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 MatSetupDM(Mat mat,DM da)
506: {
512: PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));
513: return(0);
514: }
516: PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer)
517: {
518: DM da;
519: PetscErrorCode ierr;
520: const char *prefix;
521: Mat Anatural;
522: AO ao;
523: PetscInt rstart,rend,*petsc,i;
524: IS is;
525: MPI_Comm comm;
526: PetscViewerFormat format;
529: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
530: PetscViewerGetFormat(viewer,&format);
531: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
533: PetscObjectGetComm((PetscObject)A,&comm);
534: MatGetDM(A, &da);
535: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
537: DMDAGetAO(da,&ao);
538: MatGetOwnershipRange(A,&rstart,&rend);
539: PetscMalloc1(rend-rstart,&petsc);
540: for (i=rstart; i<rend; i++) petsc[i-rstart] = i;
541: AOApplicationToPetsc(ao,rend-rstart,petsc);
542: ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);
544: /* call viewer on natural ordering */
545: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);
546: ISDestroy(&is);
547: PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);
548: PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);
549: PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);
550: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
551: MatView(Anatural,viewer);
552: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
553: MatDestroy(&Anatural);
554: return(0);
555: }
557: PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer)
558: {
559: DM da;
561: Mat Anatural,Aapp;
562: AO ao;
563: PetscInt rstart,rend,*app,i,m,n,M,N;
564: IS is;
565: MPI_Comm comm;
568: PetscObjectGetComm((PetscObject)A,&comm);
569: MatGetDM(A, &da);
570: if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA");
572: /* Load the matrix in natural ordering */
573: MatCreate(PetscObjectComm((PetscObject)A),&Anatural);
574: MatSetType(Anatural,((PetscObject)A)->type_name);
575: MatGetSize(A,&M,&N);
576: MatGetLocalSize(A,&m,&n);
577: MatSetSizes(Anatural,m,n,M,N);
578: MatLoad(Anatural,viewer);
580: /* Map natural ordering to application ordering and create IS */
581: DMDAGetAO(da,&ao);
582: MatGetOwnershipRange(Anatural,&rstart,&rend);
583: PetscMalloc1(rend-rstart,&app);
584: for (i=rstart; i<rend; i++) app[i-rstart] = i;
585: AOPetscToApplication(ao,rend-rstart,app);
586: ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);
588: /* Do permutation and replace header */
589: MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);
590: MatHeaderReplace(A,&Aapp);
591: ISDestroy(&is);
592: MatDestroy(&Anatural);
593: return(0);
594: }
596: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
597: {
599: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P;
600: Mat A;
601: MPI_Comm comm;
602: MatType Atype;
603: PetscSection section, sectionGlobal;
604: void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL,(*sell)(void)=NULL,(*is)(void)=NULL;
605: MatType mtype;
606: PetscMPIInt size;
607: DM_DA *dd = (DM_DA*)da->data;
610: MatInitializePackage();
611: mtype = da->mattype;
613: DMGetDefaultSection(da, §ion);
614: if (section) {
615: PetscInt bs = -1;
616: PetscInt localSize;
617: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
619: DMGetDefaultGlobalSection(da, §ionGlobal);
620: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
621: MatCreate(PetscObjectComm((PetscObject)da),&A);
622: MatSetSizes(A,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
623: MatSetType(A,mtype);
624: PetscStrcmp(mtype,MATSHELL,&isShell);
625: PetscStrcmp(mtype,MATBAIJ,&isBlock);
626: PetscStrcmp(mtype,MATSEQBAIJ,&isSeqBlock);
627: PetscStrcmp(mtype,MATMPIBAIJ,&isMPIBlock);
628: PetscStrcmp(mtype,MATSBAIJ,&isSymBlock);
629: PetscStrcmp(mtype,MATSEQSBAIJ,&isSymSeqBlock);
630: PetscStrcmp(mtype,MATMPISBAIJ,&isSymMPIBlock);
631: /* Check for symmetric storage */
632: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
633: if (isSymmetric) {
634: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
635: }
636: if (!isShell) {
637: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
639: if (bs < 0) {
640: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
641: PetscInt pStart, pEnd, p, dof;
643: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
644: for (p = pStart; p < pEnd; ++p) {
645: PetscSectionGetDof(sectionGlobal, p, &dof);
646: if (dof) {
647: bs = dof;
648: break;
649: }
650: }
651: } else {
652: bs = 1;
653: }
654: /* Must have same blocksize on all procs (some might have no points) */
655: bsLocal = bs;
656: MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));
657: }
658: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
659: /* DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix); */
660: PetscFree4(dnz, onz, dnzu, onzu);
661: }
662: }
663: /*
664: m
665: ------------------------------------------------------
666: | |
667: | |
668: | ---------------------- |
669: | | | |
670: n | ny | | |
671: | | | |
672: | .--------------------- |
673: | (xs,ys) nx |
674: | . |
675: | (gxs,gys) |
676: | |
677: -----------------------------------------------------
678: */
680: /*
681: nc - number of components per grid point
682: col - number of colors needed in one direction for single component problem
684: */
685: M = dd->M;
686: N = dd->N;
687: P = dd->P;
688: dim = da->dim;
689: dof = dd->w;
690: /* DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0); */
691: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
692: PetscObjectGetComm((PetscObject)da,&comm);
693: MatCreate(comm,&A);
694: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);
695: MatSetType(A,mtype);
696: MatSetDM(A,da);
697: if (da->structure_only) {
698: MatSetOption(A,MAT_STRUCTURE_ONLY,PETSC_TRUE);
699: }
700: MatGetType(A,&Atype);
701: /*
702: We do not provide a getmatrix function in the DMDA operations because
703: the basic DMDA does not know about matrices. We think of DMDA as being more
704: more low-level than matrices. This is kind of cheating but, cause sometimes
705: we think of DMDA has higher level than matrices.
707: We could switch based on Atype (or mtype), but we do not since the
708: specialized setting routines depend only the particular preallocation
709: details of the matrix, not the type itself.
710: */
711: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
712: if (!aij) {
713: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
714: }
715: if (!aij) {
716: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
717: if (!baij) {
718: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
719: }
720: if (!baij) {
721: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
722: if (!sbaij) {
723: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
724: }
725: if (!sbaij) {
726: PetscObjectQueryFunction((PetscObject)A,"MatMPISELLSetPreallocation_C",&sell);
727: if (!sell) {
728: PetscObjectQueryFunction((PetscObject)A,"MatSeqSELLSetPreallocation_C",&sell);
729: }
730: }
731: if (!sell) {
732: PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
733: }
734: }
735: }
736: if (aij) {
737: if (dim == 1) {
738: if (dd->ofill) {
739: DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);
740: } else {
741: DMCreateMatrix_DA_1d_MPIAIJ(da,A);
742: }
743: } else if (dim == 2) {
744: if (dd->ofill) {
745: DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);
746: } else {
747: DMCreateMatrix_DA_2d_MPIAIJ(da,A);
748: }
749: } else if (dim == 3) {
750: if (dd->ofill) {
751: DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);
752: } else {
753: DMCreateMatrix_DA_3d_MPIAIJ(da,A);
754: }
755: }
756: } else if (baij) {
757: if (dim == 2) {
758: DMCreateMatrix_DA_2d_MPIBAIJ(da,A);
759: } else if (dim == 3) {
760: DMCreateMatrix_DA_3d_MPIBAIJ(da,A);
761: } 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);
762: } else if (sbaij) {
763: if (dim == 2) {
764: DMCreateMatrix_DA_2d_MPISBAIJ(da,A);
765: } else if (dim == 3) {
766: DMCreateMatrix_DA_3d_MPISBAIJ(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 (sell) {
769: if (dim == 2) {
770: DMCreateMatrix_DA_2d_MPISELL(da,A);
771: } else if (dim == 3) {
772: DMCreateMatrix_DA_3d_MPISELL(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 if (is) {
775: DMCreateMatrix_DA_IS(da,A);
776: } else {
777: ISLocalToGlobalMapping ltog;
779: MatSetBlockSize(A,dof);
780: MatSetUp(A);
781: DMGetLocalToGlobalMapping(da,<og);
782: MatSetLocalToGlobalMapping(A,ltog,ltog);
783: }
784: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
785: MatSetStencil(A,dim,dims,starts,dof);
786: MatSetDM(A,da);
787: MPI_Comm_size(comm,&size);
788: if (size > 1) {
789: /* change viewer to display matrix in natural ordering */
790: MatSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);
791: MatSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);
792: }
793: MatSetFromOptions(A);
794: *J = A;
795: return(0);
796: }
798: /* ---------------------------------------------------------------------------------*/
799: PetscErrorCode DMCreateMatrix_DA_IS(DM dm,Mat J)
800: {
801: DM_DA *da = (DM_DA*)dm->data;
802: Mat lJ;
803: ISLocalToGlobalMapping ltog;
804: IS is_loc_filt, is_glob;
805: const PetscInt *e_loc,*idx;
806: PetscInt i,nel,nen,dnz,nv,dof,dim,*gidx,nb;
807: PetscErrorCode ierr;
809: /* The l2g map of DMDA has all ghosted nodes, and e_loc is a subset of all the local nodes (including the ghosted)
810: We need to filter the local indices that are represented through the DMDAGetElements decomposition
811: This is because the size of the local matrices in MATIS is the local size of the l2g map */
813: dof = da->w;
814: dim = dm->dim;
816: MatSetBlockSize(J,dof);
818: /* get local elements indices in local DMDA numbering */
819: DMDAGetElements(dm,&nel,&nen,&e_loc); /* this will throw an error if the stencil type is not DMDA_STENCIL_BOX */
820: ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nel*nen,e_loc,PETSC_COPY_VALUES,&is_loc_filt);
821: DMDARestoreElements(dm,&nel,&nen,&e_loc);
823: /* obtain a consistent local ordering for MATIS */
824: ISSortRemoveDups(is_loc_filt);
825: ISBlockGetLocalSize(is_loc_filt,&nb);
826: DMGetLocalToGlobalMapping(dm,<og);
827: ISLocalToGlobalMappingGetSize(ltog,&nv);
828: PetscMalloc1(PetscMax(nb,nv/dof),&gidx);
829: ISBlockGetIndices(is_loc_filt,&idx);
830: ISLocalToGlobalMappingApplyBlock(ltog,nb,idx,gidx);
831: ISBlockRestoreIndices(is_loc_filt,&idx);
832: ISCreateBlock(PetscObjectComm((PetscObject)dm),dof,nb,gidx,PETSC_USE_POINTER,&is_glob);
833: ISLocalToGlobalMappingCreateIS(is_glob,<og);
834: ISDestroy(&is_glob);
835: MatSetLocalToGlobalMapping(J,ltog,ltog);
836: ISLocalToGlobalMappingDestroy(<og);
838: /* We also attach a l2g map to the local matrices to have MatSetValueLocal to work */
839: MatISGetLocalMat(J,&lJ);
840: ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);
841: ISDestroy(&is_loc_filt);
842: ISCreateStride(PetscObjectComm((PetscObject)lJ),nv/dof,0,1,&is_glob);
843: ISGetIndices(is_glob,&idx);
844: ISGlobalToLocalMappingApplyBlock(ltog,IS_GTOLM_MASK,nv/dof,idx,&nb,gidx);
845: ISRestoreIndices(is_glob,&idx);
846: ISDestroy(&is_glob);
847: ISLocalToGlobalMappingDestroy(<og);
848: ISCreateBlock(PETSC_COMM_SELF,dof,nb,gidx,PETSC_USE_POINTER,&is_loc_filt);
849: ISLocalToGlobalMappingCreateIS(is_loc_filt,<og);
850: ISDestroy(&is_loc_filt);
851: MatSetLocalToGlobalMapping(lJ,ltog,ltog);
852: ISLocalToGlobalMappingDestroy(<og);
853: PetscFree(gidx);
855: /* Preallocation (not exact) */
856: switch (da->elementtype) {
857: case DMDA_ELEMENT_P1:
858: case DMDA_ELEMENT_Q1:
859: dnz = 1;
860: for (i=0; i<dim; i++) dnz *= 3;
861: dnz *= dof;
862: break;
863: default:
864: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unhandled element type %d",da->elementtype);
865: break;
866: }
867: MatSeqAIJSetPreallocation(lJ,dnz,NULL);
868: MatSeqBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);
869: MatSeqSBAIJSetPreallocation(lJ,dof,dnz/dof,NULL);
870: MatISRestoreLocalMat(J,&lJ);
871: return(0);
872: }
874: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da,Mat J)
875: {
876: PetscErrorCode ierr;
877: 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;
878: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
879: MPI_Comm comm;
880: PetscScalar *values;
881: DMBoundaryType bx,by;
882: ISLocalToGlobalMapping ltog;
883: DMDAStencilType st;
886: /*
887: nc - number of components per grid point
888: col - number of colors needed in one direction for single component problem
890: */
891: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
892: col = 2*s + 1;
893: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
894: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
895: PetscObjectGetComm((PetscObject)da,&comm);
897: PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
898: DMGetLocalToGlobalMapping(da,<og);
900: MatSetBlockSize(J,nc);
901: /* determine the matrix preallocation information */
902: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
903: for (i=xs; i<xs+nx; i++) {
905: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
906: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
908: for (j=ys; j<ys+ny; j++) {
909: slot = i - gxs + gnx*(j - gys);
911: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
912: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
914: cnt = 0;
915: for (k=0; k<nc; k++) {
916: for (l=lstart; l<lend+1; l++) {
917: for (p=pstart; p<pend+1; p++) {
918: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
919: cols[cnt++] = k + nc*(slot + gnx*l + p);
920: }
921: }
922: }
923: rows[k] = k + nc*(slot);
924: }
925: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
926: }
927: }
928: MatSetBlockSize(J,nc);
929: MatSeqSELLSetPreallocation(J,0,dnz);
930: MatMPISELLSetPreallocation(J,0,dnz,0,onz);
931: MatPreallocateFinalize(dnz,onz);
933: MatSetLocalToGlobalMapping(J,ltog,ltog);
935: /*
936: For each node in the grid: we get the neighbors in the local (on processor ordering
937: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
938: PETSc ordering.
939: */
940: if (!da->prealloc_only) {
941: PetscCalloc1(col*col*nc*nc,&values);
942: for (i=xs; i<xs+nx; i++) {
944: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
945: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
947: for (j=ys; j<ys+ny; j++) {
948: slot = i - gxs + gnx*(j - gys);
950: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
951: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
953: cnt = 0;
954: for (k=0; k<nc; k++) {
955: for (l=lstart; l<lend+1; l++) {
956: for (p=pstart; p<pend+1; p++) {
957: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
958: cols[cnt++] = k + nc*(slot + gnx*l + p);
959: }
960: }
961: }
962: rows[k] = k + nc*(slot);
963: }
964: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
965: }
966: }
967: PetscFree(values);
968: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
969: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
970: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
971: }
972: PetscFree2(rows,cols);
973: return(0);
974: }
976: PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da,Mat J)
977: {
978: PetscErrorCode ierr;
979: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
980: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
981: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
982: MPI_Comm comm;
983: PetscScalar *values;
984: DMBoundaryType bx,by,bz;
985: ISLocalToGlobalMapping ltog;
986: DMDAStencilType st;
989: /*
990: nc - number of components per grid point
991: col - number of colors needed in one direction for single component problem
993: */
994: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
995: col = 2*s + 1;
996: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
997: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
998: PetscObjectGetComm((PetscObject)da,&comm);
1000: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1001: DMGetLocalToGlobalMapping(da,<og);
1003: MatSetBlockSize(J,nc);
1004: /* determine the matrix preallocation information */
1005: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1006: for (i=xs; i<xs+nx; i++) {
1007: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1008: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1009: for (j=ys; j<ys+ny; j++) {
1010: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1011: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1012: for (k=zs; k<zs+nz; k++) {
1013: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1014: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1016: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1018: cnt = 0;
1019: for (l=0; l<nc; l++) {
1020: for (ii=istart; ii<iend+1; ii++) {
1021: for (jj=jstart; jj<jend+1; jj++) {
1022: for (kk=kstart; kk<kend+1; kk++) {
1023: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1024: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1025: }
1026: }
1027: }
1028: }
1029: rows[l] = l + nc*(slot);
1030: }
1031: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1032: }
1033: }
1034: }
1035: MatSetBlockSize(J,nc);
1036: MatSeqSELLSetPreallocation(J,0,dnz);
1037: MatMPISELLSetPreallocation(J,0,dnz,0,onz);
1038: MatPreallocateFinalize(dnz,onz);
1039: MatSetLocalToGlobalMapping(J,ltog,ltog);
1041: /*
1042: For each node in the grid: we get the neighbors in the local (on processor ordering
1043: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1044: PETSc ordering.
1045: */
1046: if (!da->prealloc_only) {
1047: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1048: for (i=xs; i<xs+nx; i++) {
1049: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1050: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1051: for (j=ys; j<ys+ny; j++) {
1052: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1053: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1054: for (k=zs; k<zs+nz; k++) {
1055: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1056: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1058: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1060: cnt = 0;
1061: for (l=0; l<nc; l++) {
1062: for (ii=istart; ii<iend+1; ii++) {
1063: for (jj=jstart; jj<jend+1; jj++) {
1064: for (kk=kstart; kk<kend+1; kk++) {
1065: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1066: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1067: }
1068: }
1069: }
1070: }
1071: rows[l] = l + nc*(slot);
1072: }
1073: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1074: }
1075: }
1076: }
1077: PetscFree(values);
1078: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1079: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1080: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1081: }
1082: PetscFree2(rows,cols);
1083: return(0);
1084: }
1086: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J)
1087: {
1088: PetscErrorCode ierr;
1089: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,M,N;
1090: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
1091: MPI_Comm comm;
1092: PetscScalar *values;
1093: DMBoundaryType bx,by;
1094: ISLocalToGlobalMapping ltog;
1095: DMDAStencilType st;
1096: PetscBool removedups = PETSC_FALSE;
1099: /*
1100: nc - number of components per grid point
1101: col - number of colors needed in one direction for single component problem
1103: */
1104: DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1105: col = 2*s + 1;
1106: /*
1107: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1108: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1109: */
1110: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1111: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1112: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1113: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1114: PetscObjectGetComm((PetscObject)da,&comm);
1116: PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);
1117: DMGetLocalToGlobalMapping(da,<og);
1119: MatSetBlockSize(J,nc);
1120: /* determine the matrix preallocation information */
1121: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1122: for (i=xs; i<xs+nx; i++) {
1124: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1125: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1127: for (j=ys; j<ys+ny; j++) {
1128: slot = i - gxs + gnx*(j - gys);
1130: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1131: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1133: cnt = 0;
1134: for (k=0; k<nc; k++) {
1135: for (l=lstart; l<lend+1; l++) {
1136: for (p=pstart; p<pend+1; p++) {
1137: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1138: cols[cnt++] = k + nc*(slot + gnx*l + p);
1139: }
1140: }
1141: }
1142: rows[k] = k + nc*(slot);
1143: }
1144: if (removedups) {
1145: MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1146: } else {
1147: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1148: }
1149: }
1150: }
1151: MatSetBlockSize(J,nc);
1152: MatSeqAIJSetPreallocation(J,0,dnz);
1153: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1154: MatPreallocateFinalize(dnz,onz);
1156: MatSetLocalToGlobalMapping(J,ltog,ltog);
1158: /*
1159: For each node in the grid: we get the neighbors in the local (on processor ordering
1160: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1161: PETSc ordering.
1162: */
1163: if (!da->prealloc_only) {
1164: PetscCalloc1(col*col*nc*nc,&values);
1165: for (i=xs; i<xs+nx; i++) {
1167: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1168: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1170: for (j=ys; j<ys+ny; j++) {
1171: slot = i - gxs + gnx*(j - gys);
1173: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1174: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1176: cnt = 0;
1177: for (k=0; k<nc; k++) {
1178: for (l=lstart; l<lend+1; l++) {
1179: for (p=pstart; p<pend+1; p++) {
1180: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
1181: cols[cnt++] = k + nc*(slot + gnx*l + p);
1182: }
1183: }
1184: }
1185: rows[k] = k + nc*(slot);
1186: }
1187: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1188: }
1189: }
1190: PetscFree(values);
1191: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1192: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1193: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1194: }
1195: PetscFree2(rows,cols);
1196: return(0);
1197: }
1199: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
1200: {
1201: PetscErrorCode ierr;
1202: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1203: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N;
1204: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
1205: DM_DA *dd = (DM_DA*)da->data;
1206: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
1207: MPI_Comm comm;
1208: PetscScalar *values;
1209: DMBoundaryType bx,by;
1210: ISLocalToGlobalMapping ltog;
1211: DMDAStencilType st;
1212: PetscBool removedups = PETSC_FALSE;
1215: /*
1216: nc - number of components per grid point
1217: col - number of colors needed in one direction for single component problem
1219: */
1220: DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);
1221: col = 2*s + 1;
1222: /*
1223: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1224: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1225: */
1226: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1227: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1228: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1229: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1230: PetscObjectGetComm((PetscObject)da,&comm);
1232: PetscMalloc1(col*col*nc,&cols);
1233: DMGetLocalToGlobalMapping(da,<og);
1235: MatSetBlockSize(J,nc);
1236: /* determine the matrix preallocation information */
1237: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
1238: for (i=xs; i<xs+nx; i++) {
1240: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1241: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1243: for (j=ys; j<ys+ny; j++) {
1244: slot = i - gxs + gnx*(j - gys);
1246: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1247: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1249: for (k=0; k<nc; k++) {
1250: cnt = 0;
1251: for (l=lstart; l<lend+1; l++) {
1252: for (p=pstart; p<pend+1; p++) {
1253: if (l || p) {
1254: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1255: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1256: }
1257: } else {
1258: if (dfill) {
1259: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1260: } else {
1261: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1262: }
1263: }
1264: }
1265: }
1266: row = k + nc*(slot);
1267: maxcnt = PetscMax(maxcnt,cnt);
1268: if (removedups) {
1269: MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1270: } else {
1271: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1272: }
1273: }
1274: }
1275: }
1276: MatSeqAIJSetPreallocation(J,0,dnz);
1277: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1278: MatPreallocateFinalize(dnz,onz);
1279: MatSetLocalToGlobalMapping(J,ltog,ltog);
1281: /*
1282: For each node in the grid: we get the neighbors in the local (on processor ordering
1283: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1284: PETSc ordering.
1285: */
1286: if (!da->prealloc_only) {
1287: PetscCalloc1(maxcnt,&values);
1288: for (i=xs; i<xs+nx; i++) {
1290: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1291: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1293: for (j=ys; j<ys+ny; j++) {
1294: slot = i - gxs + gnx*(j - gys);
1296: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1297: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1299: for (k=0; k<nc; k++) {
1300: cnt = 0;
1301: for (l=lstart; l<lend+1; l++) {
1302: for (p=pstart; p<pend+1; p++) {
1303: if (l || p) {
1304: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
1305: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
1306: }
1307: } else {
1308: if (dfill) {
1309: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
1310: } else {
1311: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
1312: }
1313: }
1314: }
1315: }
1316: row = k + nc*(slot);
1317: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1318: }
1319: }
1320: }
1321: PetscFree(values);
1322: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1323: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1324: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1325: }
1326: PetscFree(cols);
1327: return(0);
1328: }
1330: /* ---------------------------------------------------------------------------------*/
1332: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1333: {
1334: PetscErrorCode ierr;
1335: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1336: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1337: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
1338: MPI_Comm comm;
1339: PetscScalar *values;
1340: DMBoundaryType bx,by,bz;
1341: ISLocalToGlobalMapping ltog;
1342: DMDAStencilType st;
1343: PetscBool removedups = PETSC_FALSE;
1346: /*
1347: nc - number of components per grid point
1348: col - number of colors needed in one direction for single component problem
1350: */
1351: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
1352: col = 2*s + 1;
1354: /*
1355: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1356: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1357: */
1358: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
1359: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
1360: if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;
1362: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1363: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1364: PetscObjectGetComm((PetscObject)da,&comm);
1366: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1367: DMGetLocalToGlobalMapping(da,<og);
1369: MatSetBlockSize(J,nc);
1370: /* determine the matrix preallocation information */
1371: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1372: for (i=xs; i<xs+nx; i++) {
1373: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1374: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1375: for (j=ys; j<ys+ny; j++) {
1376: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1377: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1378: for (k=zs; k<zs+nz; k++) {
1379: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1380: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1382: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1384: cnt = 0;
1385: for (l=0; l<nc; l++) {
1386: for (ii=istart; ii<iend+1; ii++) {
1387: for (jj=jstart; jj<jend+1; jj++) {
1388: for (kk=kstart; kk<kend+1; kk++) {
1389: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1390: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1391: }
1392: }
1393: }
1394: }
1395: rows[l] = l + nc*(slot);
1396: }
1397: if (removedups) {
1398: MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1399: } else {
1400: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1401: }
1402: }
1403: }
1404: }
1405: MatSetBlockSize(J,nc);
1406: MatSeqAIJSetPreallocation(J,0,dnz);
1407: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1408: MatPreallocateFinalize(dnz,onz);
1409: MatSetLocalToGlobalMapping(J,ltog,ltog);
1411: /*
1412: For each node in the grid: we get the neighbors in the local (on processor ordering
1413: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1414: PETSc ordering.
1415: */
1416: if (!da->prealloc_only) {
1417: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1418: for (i=xs; i<xs+nx; i++) {
1419: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1420: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1421: for (j=ys; j<ys+ny; j++) {
1422: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1423: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1424: for (k=zs; k<zs+nz; k++) {
1425: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1426: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1428: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1430: cnt = 0;
1431: for (l=0; l<nc; l++) {
1432: for (ii=istart; ii<iend+1; ii++) {
1433: for (jj=jstart; jj<jend+1; jj++) {
1434: for (kk=kstart; kk<kend+1; kk++) {
1435: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1436: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1437: }
1438: }
1439: }
1440: }
1441: rows[l] = l + nc*(slot);
1442: }
1443: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1444: }
1445: }
1446: }
1447: PetscFree(values);
1448: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1449: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1450: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1451: }
1452: PetscFree2(rows,cols);
1453: return(0);
1454: }
1456: /* ---------------------------------------------------------------------------------*/
1458: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1459: {
1460: PetscErrorCode ierr;
1461: DM_DA *dd = (DM_DA*)da->data;
1462: PetscInt xs,nx,i,j,gxs,gnx,row,k,l;
1463: PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols;
1464: PetscInt *ofill = dd->ofill,*dfill = dd->dfill;
1465: PetscScalar *values;
1466: DMBoundaryType bx;
1467: ISLocalToGlobalMapping ltog;
1468: PetscMPIInt rank,size;
1471: if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1472: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1473: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1475: /*
1476: nc - number of components per grid point
1478: */
1479: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1480: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1481: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1483: MatSetBlockSize(J,nc);
1484: PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);
1486: /*
1487: note should be smaller for first and last process with no periodic
1488: does not handle dfill
1489: */
1490: cnt = 0;
1491: /* coupling with process to the left */
1492: for (i=0; i<s; i++) {
1493: for (j=0; j<nc; j++) {
1494: ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1495: cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1496: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1497: cnt++;
1498: }
1499: }
1500: for (i=s; i<nx-s; i++) {
1501: for (j=0; j<nc; j++) {
1502: cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1503: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1504: cnt++;
1505: }
1506: }
1507: /* coupling with process to the right */
1508: for (i=nx-s; i<nx; i++) {
1509: for (j=0; j<nc; j++) {
1510: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1511: cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1512: maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]);
1513: cnt++;
1514: }
1515: }
1517: MatSeqAIJSetPreallocation(J,0,cols);
1518: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1519: PetscFree2(cols,ocols);
1521: DMGetLocalToGlobalMapping(da,<og);
1522: MatSetLocalToGlobalMapping(J,ltog,ltog);
1524: /*
1525: For each node in the grid: we get the neighbors in the local (on processor ordering
1526: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1527: PETSc ordering.
1528: */
1529: if (!da->prealloc_only) {
1530: PetscCalloc2(maxcnt,&values,maxcnt,&cols);
1532: row = xs*nc;
1533: /* coupling with process to the left */
1534: for (i=xs; i<xs+s; i++) {
1535: for (j=0; j<nc; j++) {
1536: cnt = 0;
1537: if (rank) {
1538: for (l=0; l<s; l++) {
1539: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1540: }
1541: }
1542: if (dfill) {
1543: for (k=dfill[j]; k<dfill[j+1]; k++) {
1544: cols[cnt++] = i*nc + dfill[k];
1545: }
1546: } else {
1547: for (k=0; k<nc; k++) {
1548: cols[cnt++] = i*nc + k;
1549: }
1550: }
1551: for (l=0; l<s; l++) {
1552: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1553: }
1554: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1555: row++;
1556: }
1557: }
1558: for (i=xs+s; i<xs+nx-s; i++) {
1559: for (j=0; j<nc; j++) {
1560: cnt = 0;
1561: for (l=0; l<s; l++) {
1562: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1563: }
1564: if (dfill) {
1565: for (k=dfill[j]; k<dfill[j+1]; k++) {
1566: cols[cnt++] = i*nc + dfill[k];
1567: }
1568: } else {
1569: for (k=0; k<nc; k++) {
1570: cols[cnt++] = i*nc + k;
1571: }
1572: }
1573: for (l=0; l<s; l++) {
1574: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1575: }
1576: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1577: row++;
1578: }
1579: }
1580: /* coupling with process to the right */
1581: for (i=xs+nx-s; i<xs+nx; i++) {
1582: for (j=0; j<nc; j++) {
1583: cnt = 0;
1584: for (l=0; l<s; l++) {
1585: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1586: }
1587: if (dfill) {
1588: for (k=dfill[j]; k<dfill[j+1]; k++) {
1589: cols[cnt++] = i*nc + dfill[k];
1590: }
1591: } else {
1592: for (k=0; k<nc; k++) {
1593: cols[cnt++] = i*nc + k;
1594: }
1595: }
1596: if (rank < size-1) {
1597: for (l=0; l<s; l++) {
1598: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1599: }
1600: }
1601: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1602: row++;
1603: }
1604: }
1605: PetscFree2(values,cols);
1606: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1607: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1608: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1609: }
1610: return(0);
1611: }
1613: /* ---------------------------------------------------------------------------------*/
1615: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1616: {
1617: PetscErrorCode ierr;
1618: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1619: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1620: PetscInt istart,iend;
1621: PetscScalar *values;
1622: DMBoundaryType bx;
1623: ISLocalToGlobalMapping ltog;
1626: /*
1627: nc - number of components per grid point
1628: col - number of colors needed in one direction for single component problem
1630: */
1631: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1632: col = 2*s + 1;
1634: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1635: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1637: MatSetBlockSize(J,nc);
1638: MatSeqAIJSetPreallocation(J,col*nc,0);
1639: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1641: DMGetLocalToGlobalMapping(da,<og);
1642: MatSetLocalToGlobalMapping(J,ltog,ltog);
1644: /*
1645: For each node in the grid: we get the neighbors in the local (on processor ordering
1646: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1647: PETSc ordering.
1648: */
1649: if (!da->prealloc_only) {
1650: PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1651: PetscCalloc1(col*nc*nc,&values);
1652: for (i=xs; i<xs+nx; i++) {
1653: istart = PetscMax(-s,gxs - i);
1654: iend = PetscMin(s,gxs + gnx - i - 1);
1655: slot = i - gxs;
1657: cnt = 0;
1658: for (l=0; l<nc; l++) {
1659: for (i1=istart; i1<iend+1; i1++) {
1660: cols[cnt++] = l + nc*(slot + i1);
1661: }
1662: rows[l] = l + nc*(slot);
1663: }
1664: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1665: }
1666: PetscFree(values);
1667: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1668: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1669: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1670: PetscFree2(rows,cols);
1671: }
1672: return(0);
1673: }
1675: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1676: {
1677: PetscErrorCode ierr;
1678: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1679: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1680: PetscInt istart,iend,jstart,jend,ii,jj;
1681: MPI_Comm comm;
1682: PetscScalar *values;
1683: DMBoundaryType bx,by;
1684: DMDAStencilType st;
1685: ISLocalToGlobalMapping ltog;
1688: /*
1689: nc - number of components per grid point
1690: col - number of colors needed in one direction for single component problem
1691: */
1692: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1693: col = 2*s + 1;
1695: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1696: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1697: PetscObjectGetComm((PetscObject)da,&comm);
1699: PetscMalloc1(col*col*nc*nc,&cols);
1701: DMGetLocalToGlobalMapping(da,<og);
1703: /* determine the matrix preallocation information */
1704: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1705: for (i=xs; i<xs+nx; i++) {
1706: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1707: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1708: for (j=ys; j<ys+ny; j++) {
1709: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1710: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1711: slot = i - gxs + gnx*(j - gys);
1713: /* Find block columns in block row */
1714: cnt = 0;
1715: for (ii=istart; ii<iend+1; ii++) {
1716: for (jj=jstart; jj<jend+1; jj++) {
1717: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1718: cols[cnt++] = slot + ii + gnx*jj;
1719: }
1720: }
1721: }
1722: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1723: }
1724: }
1725: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1726: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1727: MatPreallocateFinalize(dnz,onz);
1729: MatSetLocalToGlobalMapping(J,ltog,ltog);
1731: /*
1732: For each node in the grid: we get the neighbors in the local (on processor ordering
1733: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1734: PETSc ordering.
1735: */
1736: if (!da->prealloc_only) {
1737: PetscCalloc1(col*col*nc*nc,&values);
1738: for (i=xs; i<xs+nx; i++) {
1739: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1740: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1741: for (j=ys; j<ys+ny; j++) {
1742: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1743: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1744: slot = i - gxs + gnx*(j - gys);
1745: cnt = 0;
1746: for (ii=istart; ii<iend+1; ii++) {
1747: for (jj=jstart; jj<jend+1; jj++) {
1748: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1749: cols[cnt++] = slot + ii + gnx*jj;
1750: }
1751: }
1752: }
1753: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1754: }
1755: }
1756: PetscFree(values);
1757: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1758: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1759: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1760: }
1761: PetscFree(cols);
1762: return(0);
1763: }
1765: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1766: {
1767: PetscErrorCode ierr;
1768: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1769: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1770: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1771: MPI_Comm comm;
1772: PetscScalar *values;
1773: DMBoundaryType bx,by,bz;
1774: DMDAStencilType st;
1775: ISLocalToGlobalMapping ltog;
1778: /*
1779: nc - number of components per grid point
1780: col - number of colors needed in one direction for single component problem
1782: */
1783: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1784: col = 2*s + 1;
1786: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1787: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1788: PetscObjectGetComm((PetscObject)da,&comm);
1790: PetscMalloc1(col*col*col,&cols);
1792: DMGetLocalToGlobalMapping(da,<og);
1794: /* determine the matrix preallocation information */
1795: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1796: for (i=xs; i<xs+nx; i++) {
1797: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1798: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1799: for (j=ys; j<ys+ny; j++) {
1800: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1801: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1802: for (k=zs; k<zs+nz; k++) {
1803: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1804: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1806: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1808: /* Find block columns in block row */
1809: cnt = 0;
1810: for (ii=istart; ii<iend+1; ii++) {
1811: for (jj=jstart; jj<jend+1; jj++) {
1812: for (kk=kstart; kk<kend+1; kk++) {
1813: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1814: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1815: }
1816: }
1817: }
1818: }
1819: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1820: }
1821: }
1822: }
1823: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1824: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1825: MatPreallocateFinalize(dnz,onz);
1827: MatSetLocalToGlobalMapping(J,ltog,ltog);
1829: /*
1830: For each node in the grid: we get the neighbors in the local (on processor ordering
1831: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1832: PETSc ordering.
1833: */
1834: if (!da->prealloc_only) {
1835: PetscCalloc1(col*col*col*nc*nc,&values);
1836: for (i=xs; i<xs+nx; i++) {
1837: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1838: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1839: for (j=ys; j<ys+ny; j++) {
1840: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1841: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1842: for (k=zs; k<zs+nz; k++) {
1843: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1844: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1846: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1848: cnt = 0;
1849: for (ii=istart; ii<iend+1; ii++) {
1850: for (jj=jstart; jj<jend+1; jj++) {
1851: for (kk=kstart; kk<kend+1; kk++) {
1852: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1853: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1854: }
1855: }
1856: }
1857: }
1858: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1859: }
1860: }
1861: }
1862: PetscFree(values);
1863: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1864: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1865: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1866: }
1867: PetscFree(cols);
1868: return(0);
1869: }
1871: /*
1872: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1873: identify in the local ordering with periodic domain.
1874: */
1875: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1876: {
1878: PetscInt i,n;
1881: ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1882: ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1883: for (i=0,n=0; i<*cnt; i++) {
1884: if (col[i] >= *row) col[n++] = col[i];
1885: }
1886: *cnt = n;
1887: return(0);
1888: }
1890: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1891: {
1892: PetscErrorCode ierr;
1893: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1894: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1895: PetscInt istart,iend,jstart,jend,ii,jj;
1896: MPI_Comm comm;
1897: PetscScalar *values;
1898: DMBoundaryType bx,by;
1899: DMDAStencilType st;
1900: ISLocalToGlobalMapping ltog;
1903: /*
1904: nc - number of components per grid point
1905: col - number of colors needed in one direction for single component problem
1906: */
1907: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1908: col = 2*s + 1;
1910: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1911: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1912: PetscObjectGetComm((PetscObject)da,&comm);
1914: PetscMalloc1(col*col*nc*nc,&cols);
1916: DMGetLocalToGlobalMapping(da,<og);
1918: /* determine the matrix preallocation information */
1919: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1920: for (i=xs; i<xs+nx; i++) {
1921: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1922: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1923: for (j=ys; j<ys+ny; j++) {
1924: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1925: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1926: slot = i - gxs + gnx*(j - gys);
1928: /* Find block columns in block row */
1929: cnt = 0;
1930: for (ii=istart; ii<iend+1; ii++) {
1931: for (jj=jstart; jj<jend+1; jj++) {
1932: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1933: cols[cnt++] = slot + ii + gnx*jj;
1934: }
1935: }
1936: }
1937: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1938: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1939: }
1940: }
1941: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1942: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1943: MatPreallocateFinalize(dnz,onz);
1945: MatSetLocalToGlobalMapping(J,ltog,ltog);
1947: /*
1948: For each node in the grid: we get the neighbors in the local (on processor ordering
1949: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1950: PETSc ordering.
1951: */
1952: if (!da->prealloc_only) {
1953: PetscCalloc1(col*col*nc*nc,&values);
1954: for (i=xs; i<xs+nx; i++) {
1955: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1956: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1957: for (j=ys; j<ys+ny; j++) {
1958: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1959: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1960: slot = i - gxs + gnx*(j - gys);
1962: /* Find block columns in block row */
1963: cnt = 0;
1964: for (ii=istart; ii<iend+1; ii++) {
1965: for (jj=jstart; jj<jend+1; jj++) {
1966: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1967: cols[cnt++] = slot + ii + gnx*jj;
1968: }
1969: }
1970: }
1971: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1972: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1973: }
1974: }
1975: PetscFree(values);
1976: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1977: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1978: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1979: }
1980: PetscFree(cols);
1981: return(0);
1982: }
1984: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1985: {
1986: PetscErrorCode ierr;
1987: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1988: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1989: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1990: MPI_Comm comm;
1991: PetscScalar *values;
1992: DMBoundaryType bx,by,bz;
1993: DMDAStencilType st;
1994: ISLocalToGlobalMapping ltog;
1997: /*
1998: nc - number of components per grid point
1999: col - number of colors needed in one direction for single component problem
2000: */
2001: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
2002: col = 2*s + 1;
2004: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2005: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2006: PetscObjectGetComm((PetscObject)da,&comm);
2008: /* create the matrix */
2009: PetscMalloc1(col*col*col,&cols);
2011: DMGetLocalToGlobalMapping(da,<og);
2013: /* determine the matrix preallocation information */
2014: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
2015: for (i=xs; i<xs+nx; i++) {
2016: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2017: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2018: for (j=ys; j<ys+ny; j++) {
2019: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2020: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2021: for (k=zs; k<zs+nz; k++) {
2022: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2023: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2025: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2027: /* Find block columns in block row */
2028: cnt = 0;
2029: for (ii=istart; ii<iend+1; ii++) {
2030: for (jj=jstart; jj<jend+1; jj++) {
2031: for (kk=kstart; kk<kend+1; kk++) {
2032: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2033: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2034: }
2035: }
2036: }
2037: }
2038: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2039: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
2040: }
2041: }
2042: }
2043: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
2044: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
2045: MatPreallocateFinalize(dnz,onz);
2047: MatSetLocalToGlobalMapping(J,ltog,ltog);
2049: /*
2050: For each node in the grid: we get the neighbors in the local (on processor ordering
2051: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2052: PETSc ordering.
2053: */
2054: if (!da->prealloc_only) {
2055: PetscCalloc1(col*col*col*nc*nc,&values);
2056: for (i=xs; i<xs+nx; i++) {
2057: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2058: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2059: for (j=ys; j<ys+ny; j++) {
2060: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2061: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2062: for (k=zs; k<zs+nz; k++) {
2063: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2064: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2066: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2068: cnt = 0;
2069: for (ii=istart; ii<iend+1; ii++) {
2070: for (jj=jstart; jj<jend+1; jj++) {
2071: for (kk=kstart; kk<kend+1; kk++) {
2072: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
2073: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
2074: }
2075: }
2076: }
2077: }
2078: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
2079: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
2080: }
2081: }
2082: }
2083: PetscFree(values);
2084: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2085: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2086: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2087: }
2088: PetscFree(cols);
2089: return(0);
2090: }
2092: /* ---------------------------------------------------------------------------------*/
2094: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
2095: {
2096: PetscErrorCode ierr;
2097: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
2098: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz;
2099: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P;
2100: DM_DA *dd = (DM_DA*)da->data;
2101: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
2102: MPI_Comm comm;
2103: PetscScalar *values;
2104: DMBoundaryType bx,by,bz;
2105: ISLocalToGlobalMapping ltog;
2106: DMDAStencilType st;
2107: PetscBool removedups = PETSC_FALSE;
2110: /*
2111: nc - number of components per grid point
2112: col - number of colors needed in one direction for single component problem
2114: */
2115: DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);
2116: col = 2*s + 1;
2117: 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\
2118: by 2*stencil_width + 1\n");
2119: 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\
2120: by 2*stencil_width + 1\n");
2121: 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\
2122: by 2*stencil_width + 1\n");
2124: /*
2125: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2126: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2127: */
2128: if (M == 1 && 2*s >= m) removedups = PETSC_TRUE;
2129: if (N == 1 && 2*s >= n) removedups = PETSC_TRUE;
2130: if (P == 1 && 2*s >= p) removedups = PETSC_TRUE;
2132: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
2133: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
2134: PetscObjectGetComm((PetscObject)da,&comm);
2136: PetscMalloc1(col*col*col*nc,&cols);
2137: DMGetLocalToGlobalMapping(da,<og);
2139: /* determine the matrix preallocation information */
2140: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
2142: MatSetBlockSize(J,nc);
2143: for (i=xs; i<xs+nx; i++) {
2144: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2145: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2146: for (j=ys; j<ys+ny; j++) {
2147: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2148: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2149: for (k=zs; k<zs+nz; k++) {
2150: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2151: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2153: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2155: for (l=0; l<nc; l++) {
2156: cnt = 0;
2157: for (ii=istart; ii<iend+1; ii++) {
2158: for (jj=jstart; jj<jend+1; jj++) {
2159: for (kk=kstart; kk<kend+1; kk++) {
2160: if (ii || jj || kk) {
2161: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2162: 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);
2163: }
2164: } else {
2165: if (dfill) {
2166: 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);
2167: } else {
2168: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2169: }
2170: }
2171: }
2172: }
2173: }
2174: row = l + nc*(slot);
2175: maxcnt = PetscMax(maxcnt,cnt);
2176: if (removedups) {
2177: MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2178: } else {
2179: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
2180: }
2181: }
2182: }
2183: }
2184: }
2185: MatSeqAIJSetPreallocation(J,0,dnz);
2186: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
2187: MatPreallocateFinalize(dnz,onz);
2188: MatSetLocalToGlobalMapping(J,ltog,ltog);
2190: /*
2191: For each node in the grid: we get the neighbors in the local (on processor ordering
2192: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2193: PETSc ordering.
2194: */
2195: if (!da->prealloc_only) {
2196: PetscCalloc1(maxcnt,&values);
2197: for (i=xs; i<xs+nx; i++) {
2198: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
2199: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
2200: for (j=ys; j<ys+ny; j++) {
2201: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
2202: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
2203: for (k=zs; k<zs+nz; k++) {
2204: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
2205: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
2207: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
2209: for (l=0; l<nc; l++) {
2210: cnt = 0;
2211: for (ii=istart; ii<iend+1; ii++) {
2212: for (jj=jstart; jj<jend+1; jj++) {
2213: for (kk=kstart; kk<kend+1; kk++) {
2214: if (ii || jj || kk) {
2215: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
2216: 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);
2217: }
2218: } else {
2219: if (dfill) {
2220: 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);
2221: } else {
2222: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
2223: }
2224: }
2225: }
2226: }
2227: }
2228: row = l + nc*(slot);
2229: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
2230: }
2231: }
2232: }
2233: }
2234: PetscFree(values);
2235: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2236: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2237: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2238: }
2239: PetscFree(cols);
2240: return(0);
2241: }