Actual source code: fdda.c
petsc-3.5.4 2015-05-23
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,&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,&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,&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,&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,&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,&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,&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,&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: MatView(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: MPI_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 = dd->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: /* determine the matrix preallocation information */
823: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
824: for (i=xs; i<xs+nx; i++) {
826: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
827: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
829: for (j=ys; j<ys+ny; j++) {
830: slot = i - gxs + gnx*(j - gys);
832: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
833: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
835: cnt = 0;
836: for (k=0; k<nc; k++) {
837: for (l=lstart; l<lend+1; l++) {
838: for (p=pstart; p<pend+1; p++) {
839: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
840: cols[cnt++] = k + nc*(slot + gnx*l + p);
841: }
842: }
843: }
844: rows[k] = k + nc*(slot);
845: }
846: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
847: }
848: }
849: MatSetBlockSize(J,nc);
850: MatSeqAIJSetPreallocation(J,0,dnz);
851: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
852: MatPreallocateFinalize(dnz,onz);
854: MatSetLocalToGlobalMapping(J,ltog,ltog);
856: /*
857: For each node in the grid: we get the neighbors in the local (on processor ordering
858: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
859: PETSc ordering.
860: */
861: if (!da->prealloc_only) {
862: PetscCalloc1(col*col*nc*nc,&values);
863: for (i=xs; i<xs+nx; i++) {
865: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
866: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
868: for (j=ys; j<ys+ny; j++) {
869: slot = i - gxs + gnx*(j - gys);
871: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
872: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
874: cnt = 0;
875: for (k=0; k<nc; k++) {
876: for (l=lstart; l<lend+1; l++) {
877: for (p=pstart; p<pend+1; p++) {
878: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
879: cols[cnt++] = k + nc*(slot + gnx*l + p);
880: }
881: }
882: }
883: rows[k] = k + nc*(slot);
884: }
885: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
886: }
887: }
888: PetscFree(values);
889: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
890: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
891: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
892: }
893: PetscFree2(rows,cols);
894: return(0);
895: }
899: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J)
900: {
901: PetscErrorCode ierr;
902: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
903: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p;
904: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
905: DM_DA *dd = (DM_DA*)da->data;
906: PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill;
907: MPI_Comm comm;
908: PetscScalar *values;
909: DMBoundaryType bx,by;
910: ISLocalToGlobalMapping ltog;
911: DMDAStencilType st;
914: /*
915: nc - number of components per grid point
916: col - number of colors needed in one direction for single component problem
918: */
919: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
920: col = 2*s + 1;
921: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
922: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
923: PetscObjectGetComm((PetscObject)da,&comm);
925: PetscMalloc1(col*col*nc*nc,&cols);
926: DMGetLocalToGlobalMapping(da,<og);
928: /* determine the matrix preallocation information */
929: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
930: for (i=xs; i<xs+nx; i++) {
932: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
933: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
935: for (j=ys; j<ys+ny; j++) {
936: slot = i - gxs + gnx*(j - gys);
938: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
939: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
941: for (k=0; k<nc; k++) {
942: cnt = 0;
943: for (l=lstart; l<lend+1; l++) {
944: for (p=pstart; p<pend+1; p++) {
945: if (l || p) {
946: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
947: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
948: }
949: } else {
950: if (dfill) {
951: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
952: } else {
953: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
954: }
955: }
956: }
957: }
958: row = k + nc*(slot);
959: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
960: }
961: }
962: }
963: MatSeqAIJSetPreallocation(J,0,dnz);
964: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
965: MatPreallocateFinalize(dnz,onz);
966: MatSetLocalToGlobalMapping(J,ltog,ltog);
968: /*
969: For each node in the grid: we get the neighbors in the local (on processor ordering
970: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
971: PETSc ordering.
972: */
973: if (!da->prealloc_only) {
974: PetscCalloc1(col*col*nc*nc,&values);
975: for (i=xs; i<xs+nx; i++) {
977: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
978: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
980: for (j=ys; j<ys+ny; j++) {
981: slot = i - gxs + gnx*(j - gys);
983: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
984: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
986: for (k=0; k<nc; k++) {
987: cnt = 0;
988: for (l=lstart; l<lend+1; l++) {
989: for (p=pstart; p<pend+1; p++) {
990: if (l || p) {
991: if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
992: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
993: }
994: } else {
995: if (dfill) {
996: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
997: } else {
998: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
999: }
1000: }
1001: }
1002: }
1003: row = k + nc*(slot);
1004: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1005: }
1006: }
1007: }
1008: PetscFree(values);
1009: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1010: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1011: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1012: }
1013: PetscFree(cols);
1014: return(0);
1015: }
1017: /* ---------------------------------------------------------------------------------*/
1021: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)
1022: {
1023: PetscErrorCode ierr;
1024: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1025: PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL;
1026: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1027: MPI_Comm comm;
1028: PetscScalar *values;
1029: DMBoundaryType bx,by,bz;
1030: ISLocalToGlobalMapping ltog;
1031: DMDAStencilType st;
1034: /*
1035: nc - number of components per grid point
1036: col - number of colors needed in one direction for single component problem
1038: */
1039: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1040: col = 2*s + 1;
1042: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1043: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1044: PetscObjectGetComm((PetscObject)da,&comm);
1046: PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);
1047: DMGetLocalToGlobalMapping(da,<og);
1049: /* determine the matrix preallocation information */
1050: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1051: for (i=xs; i<xs+nx; i++) {
1052: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1053: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1054: for (j=ys; j<ys+ny; j++) {
1055: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1056: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1057: for (k=zs; k<zs+nz; k++) {
1058: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1059: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1061: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1063: cnt = 0;
1064: for (l=0; l<nc; l++) {
1065: for (ii=istart; ii<iend+1; ii++) {
1066: for (jj=jstart; jj<jend+1; jj++) {
1067: for (kk=kstart; kk<kend+1; kk++) {
1068: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1069: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1070: }
1071: }
1072: }
1073: }
1074: rows[l] = l + nc*(slot);
1075: }
1076: MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);
1077: }
1078: }
1079: }
1080: MatSetBlockSize(J,nc);
1081: MatSeqAIJSetPreallocation(J,0,dnz);
1082: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1083: MatPreallocateFinalize(dnz,onz);
1084: MatSetLocalToGlobalMapping(J,ltog,ltog);
1086: /*
1087: For each node in the grid: we get the neighbors in the local (on processor ordering
1088: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1089: PETSc ordering.
1090: */
1091: if (!da->prealloc_only) {
1092: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1093: for (i=xs; i<xs+nx; i++) {
1094: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1095: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1096: for (j=ys; j<ys+ny; j++) {
1097: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1098: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1099: for (k=zs; k<zs+nz; k++) {
1100: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1101: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1103: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1105: cnt = 0;
1106: for (l=0; l<nc; l++) {
1107: for (ii=istart; ii<iend+1; ii++) {
1108: for (jj=jstart; jj<jend+1; jj++) {
1109: for (kk=kstart; kk<kend+1; kk++) {
1110: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1111: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1112: }
1113: }
1114: }
1115: }
1116: rows[l] = l + nc*(slot);
1117: }
1118: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1119: }
1120: }
1121: }
1122: PetscFree(values);
1123: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1124: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1125: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1126: }
1127: PetscFree2(rows,cols);
1128: return(0);
1129: }
1131: /* ---------------------------------------------------------------------------------*/
1135: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J)
1136: {
1137: PetscErrorCode ierr;
1138: DM_DA *dd = (DM_DA*)da->data;
1139: PetscInt xs,nx,i,j,gxs,gnx,row,k,l;
1140: PetscInt m,dim,s,*cols = NULL,nc,col,cnt,*ocols;
1141: PetscInt *ofill = dd->ofill,*dfill = dd->dfill;
1142: PetscScalar *values;
1143: DMBoundaryType bx;
1144: ISLocalToGlobalMapping ltog;
1145: PetscMPIInt rank,size;
1148: if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions");
1149: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
1150: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
1152: /*
1153: nc - number of components per grid point
1154: col - number of colors needed in one direction for single component problem
1156: */
1157: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1158: col = 2*s + 1;
1160: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1161: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1163: MatSetBlockSize(J,nc);
1164: PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);
1166: if (nx < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Need at least two grid points per process");
1167: /*
1168: note should be smaller for first and last process with no periodic
1169: does not handle dfill
1170: */
1171: cnt = 0;
1172: /* coupling with process to the left */
1173: for (i=0; i<s; i++) {
1174: for (j=0; j<nc; j++) {
1175: ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j]));
1176: cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]);
1177: cnt++;
1178: }
1179: }
1180: for (i=s; i<nx-s; i++) {
1181: for (j=0; j<nc; j++) {
1182: cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]);
1183: cnt++;
1184: }
1185: }
1186: /* coupling with process to the right */
1187: for (i=nx-s; i<nx; i++) {
1188: for (j=0; j<nc; j++) {
1189: ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j]));
1190: cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]);
1191: cnt++;
1192: }
1193: }
1195: MatSeqAIJSetPreallocation(J,0,cols);
1196: MatMPIAIJSetPreallocation(J,0,cols,0,ocols);
1197: PetscFree2(cols,ocols);
1199: DMGetLocalToGlobalMapping(da,<og);
1200: MatSetLocalToGlobalMapping(J,ltog,ltog);
1202: /*
1203: For each node in the grid: we get the neighbors in the local (on processor ordering
1204: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1205: PETSc ordering.
1206: */
1207: if (!da->prealloc_only) {
1208: PetscMalloc1(col*nc*nc,&cols);
1209: PetscCalloc1(col*nc*nc,&values);
1211: row = xs*nc;
1212: /* coupling with process to the left */
1213: for (i=xs; i<xs+s; i++) {
1214: for (j=0; j<nc; j++) {
1215: cnt = 0;
1216: if (rank) {
1217: for (l=0; l<s; l++) {
1218: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1219: }
1220: }
1221: if (dfill) {
1222: for (k=dfill[j]; k<dfill[j+1]; k++) {
1223: cols[cnt++] = i*nc + dfill[k];
1224: }
1225: } else {
1226: for (k=0; k<nc; k++) {
1227: cols[cnt++] = i*nc + k;
1228: }
1229: }
1230: for (l=0; l<s; l++) {
1231: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1232: }
1233: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1234: row++;
1235: }
1236: }
1237: for (i=xs+s; i<xs+nx-s; i++) {
1238: for (j=0; j<nc; j++) {
1239: cnt = 0;
1240: for (l=0; l<s; l++) {
1241: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1242: }
1243: if (dfill) {
1244: for (k=dfill[j]; k<dfill[j+1]; k++) {
1245: cols[cnt++] = i*nc + dfill[k];
1246: }
1247: } else {
1248: for (k=0; k<nc; k++) {
1249: cols[cnt++] = i*nc + k;
1250: }
1251: }
1252: for (l=0; l<s; l++) {
1253: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1254: }
1255: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1256: row++;
1257: }
1258: }
1259: /* coupling with process to the right */
1260: for (i=xs+nx-s; i<xs+nx; i++) {
1261: for (j=0; j<nc; j++) {
1262: cnt = 0;
1263: for (l=0; l<s; l++) {
1264: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k];
1265: }
1266: if (dfill) {
1267: for (k=dfill[j]; k<dfill[j+1]; k++) {
1268: cols[cnt++] = i*nc + dfill[k];
1269: }
1270: } else {
1271: for (k=0; k<nc; k++) {
1272: cols[cnt++] = i*nc + k;
1273: }
1274: }
1275: if (rank < size-1) {
1276: for (l=0; l<s; l++) {
1277: for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k];
1278: }
1279: }
1280: MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);
1281: row++;
1282: }
1283: }
1284: PetscFree(values);
1285: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1286: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1287: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1288: PetscFree(cols);
1289: }
1290: return(0);
1291: }
1293: /* ---------------------------------------------------------------------------------*/
1297: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J)
1298: {
1299: PetscErrorCode ierr;
1300: PetscInt xs,nx,i,i1,slot,gxs,gnx;
1301: PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l;
1302: PetscInt istart,iend;
1303: PetscScalar *values;
1304: DMBoundaryType bx;
1305: ISLocalToGlobalMapping ltog;
1308: /*
1309: nc - number of components per grid point
1310: col - number of colors needed in one direction for single component problem
1312: */
1313: DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);
1314: col = 2*s + 1;
1316: DMDAGetCorners(da,&xs,0,0,&nx,0,0);
1317: DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
1319: MatSetBlockSize(J,nc);
1320: MatSeqAIJSetPreallocation(J,col*nc,0);
1321: MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);
1323: DMGetLocalToGlobalMapping(da,<og);
1324: MatSetLocalToGlobalMapping(J,ltog,ltog);
1326: /*
1327: For each node in the grid: we get the neighbors in the local (on processor ordering
1328: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1329: PETSc ordering.
1330: */
1331: if (!da->prealloc_only) {
1332: PetscMalloc2(nc,&rows,col*nc*nc,&cols);
1333: PetscCalloc1(col*nc*nc,&values);
1334: for (i=xs; i<xs+nx; i++) {
1335: istart = PetscMax(-s,gxs - i);
1336: iend = PetscMin(s,gxs + gnx - i - 1);
1337: slot = i - gxs;
1339: cnt = 0;
1340: for (l=0; l<nc; l++) {
1341: for (i1=istart; i1<iend+1; i1++) {
1342: cols[cnt++] = l + nc*(slot + i1);
1343: }
1344: rows[l] = l + nc*(slot);
1345: }
1346: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
1347: }
1348: PetscFree(values);
1349: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1350: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1351: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1352: PetscFree2(rows,cols);
1353: }
1354: return(0);
1355: }
1359: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J)
1360: {
1361: PetscErrorCode ierr;
1362: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1363: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1364: PetscInt istart,iend,jstart,jend,ii,jj;
1365: MPI_Comm comm;
1366: PetscScalar *values;
1367: DMBoundaryType bx,by;
1368: DMDAStencilType st;
1369: ISLocalToGlobalMapping ltog;
1372: /*
1373: nc - number of components per grid point
1374: col - number of colors needed in one direction for single component problem
1375: */
1376: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1377: col = 2*s + 1;
1379: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1380: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1381: PetscObjectGetComm((PetscObject)da,&comm);
1383: PetscMalloc1(col*col*nc*nc,&cols);
1385: DMGetLocalToGlobalMapping(da,<og);
1387: /* determine the matrix preallocation information */
1388: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1389: for (i=xs; i<xs+nx; i++) {
1390: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1391: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1392: for (j=ys; j<ys+ny; j++) {
1393: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1394: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1395: slot = i - gxs + gnx*(j - gys);
1397: /* Find block columns in block row */
1398: cnt = 0;
1399: for (ii=istart; ii<iend+1; ii++) {
1400: for (jj=jstart; jj<jend+1; jj++) {
1401: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1402: cols[cnt++] = slot + ii + gnx*jj;
1403: }
1404: }
1405: }
1406: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1407: }
1408: }
1409: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1410: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1411: MatPreallocateFinalize(dnz,onz);
1413: MatSetLocalToGlobalMapping(J,ltog,ltog);
1415: /*
1416: For each node in the grid: we get the neighbors in the local (on processor ordering
1417: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1418: PETSc ordering.
1419: */
1420: if (!da->prealloc_only) {
1421: PetscCalloc1(col*col*nc*nc,&values);
1422: for (i=xs; i<xs+nx; i++) {
1423: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1424: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1425: for (j=ys; j<ys+ny; j++) {
1426: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1427: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1428: slot = i - gxs + gnx*(j - gys);
1429: cnt = 0;
1430: for (ii=istart; ii<iend+1; ii++) {
1431: for (jj=jstart; jj<jend+1; jj++) {
1432: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1433: cols[cnt++] = slot + ii + gnx*jj;
1434: }
1435: }
1436: }
1437: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1438: }
1439: }
1440: PetscFree(values);
1441: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1442: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1443: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1444: }
1445: PetscFree(cols);
1446: return(0);
1447: }
1451: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J)
1452: {
1453: PetscErrorCode ierr;
1454: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1455: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1456: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1457: MPI_Comm comm;
1458: PetscScalar *values;
1459: DMBoundaryType bx,by,bz;
1460: DMDAStencilType st;
1461: ISLocalToGlobalMapping ltog;
1464: /*
1465: nc - number of components per grid point
1466: col - number of colors needed in one direction for single component problem
1468: */
1469: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1470: col = 2*s + 1;
1472: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1473: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1474: PetscObjectGetComm((PetscObject)da,&comm);
1476: PetscMalloc1(col*col*col,&cols);
1478: DMGetLocalToGlobalMapping(da,<og);
1480: /* determine the matrix preallocation information */
1481: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1482: for (i=xs; i<xs+nx; i++) {
1483: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1484: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1485: for (j=ys; j<ys+ny; j++) {
1486: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1487: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1488: for (k=zs; k<zs+nz; k++) {
1489: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1490: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1492: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1494: /* Find block columns in block row */
1495: cnt = 0;
1496: for (ii=istart; ii<iend+1; ii++) {
1497: for (jj=jstart; jj<jend+1; jj++) {
1498: for (kk=kstart; kk<kend+1; kk++) {
1499: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1500: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1501: }
1502: }
1503: }
1504: }
1505: MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);
1506: }
1507: }
1508: }
1509: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1510: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1511: MatPreallocateFinalize(dnz,onz);
1513: MatSetLocalToGlobalMapping(J,ltog,ltog);
1515: /*
1516: For each node in the grid: we get the neighbors in the local (on processor ordering
1517: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1518: PETSc ordering.
1519: */
1520: if (!da->prealloc_only) {
1521: PetscCalloc1(col*col*col*nc*nc,&values);
1522: for (i=xs; i<xs+nx; i++) {
1523: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1524: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1525: for (j=ys; j<ys+ny; j++) {
1526: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1527: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1528: for (k=zs; k<zs+nz; k++) {
1529: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1530: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1532: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1534: cnt = 0;
1535: for (ii=istart; ii<iend+1; ii++) {
1536: for (jj=jstart; jj<jend+1; jj++) {
1537: for (kk=kstart; kk<kend+1; kk++) {
1538: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1539: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1540: }
1541: }
1542: }
1543: }
1544: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1545: }
1546: }
1547: }
1548: PetscFree(values);
1549: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1550: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1551: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1552: }
1553: PetscFree(cols);
1554: return(0);
1555: }
1559: /*
1560: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1561: identify in the local ordering with periodic domain.
1562: */
1563: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[])
1564: {
1566: PetscInt i,n;
1569: ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);
1570: ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);
1571: for (i=0,n=0; i<*cnt; i++) {
1572: if (col[i] >= *row) col[n++] = col[i];
1573: }
1574: *cnt = n;
1575: return(0);
1576: }
1580: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J)
1581: {
1582: PetscErrorCode ierr;
1583: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1584: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1585: PetscInt istart,iend,jstart,jend,ii,jj;
1586: MPI_Comm comm;
1587: PetscScalar *values;
1588: DMBoundaryType bx,by;
1589: DMDAStencilType st;
1590: ISLocalToGlobalMapping ltog;
1593: /*
1594: nc - number of components per grid point
1595: col - number of colors needed in one direction for single component problem
1596: */
1597: DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);
1598: col = 2*s + 1;
1600: DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1601: DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1602: PetscObjectGetComm((PetscObject)da,&comm);
1604: PetscMalloc1(col*col*nc*nc,&cols);
1606: DMGetLocalToGlobalMapping(da,<og);
1608: /* determine the matrix preallocation information */
1609: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1610: for (i=xs; i<xs+nx; i++) {
1611: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1612: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1613: for (j=ys; j<ys+ny; j++) {
1614: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1615: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1616: slot = i - gxs + gnx*(j - gys);
1618: /* Find block columns in block row */
1619: cnt = 0;
1620: for (ii=istart; ii<iend+1; ii++) {
1621: for (jj=jstart; jj<jend+1; jj++) {
1622: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1623: cols[cnt++] = slot + ii + gnx*jj;
1624: }
1625: }
1626: }
1627: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1628: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1629: }
1630: }
1631: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1632: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1633: MatPreallocateFinalize(dnz,onz);
1635: MatSetLocalToGlobalMapping(J,ltog,ltog);
1637: /*
1638: For each node in the grid: we get the neighbors in the local (on processor ordering
1639: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1640: PETSc ordering.
1641: */
1642: if (!da->prealloc_only) {
1643: PetscCalloc1(col*col*nc*nc,&values);
1644: for (i=xs; i<xs+nx; i++) {
1645: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1646: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1647: for (j=ys; j<ys+ny; j++) {
1648: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1649: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1650: slot = i - gxs + gnx*(j - gys);
1652: /* Find block columns in block row */
1653: cnt = 0;
1654: for (ii=istart; ii<iend+1; ii++) {
1655: for (jj=jstart; jj<jend+1; jj++) {
1656: if (st == DMDA_STENCIL_BOX || !ii || !jj) {
1657: cols[cnt++] = slot + ii + gnx*jj;
1658: }
1659: }
1660: }
1661: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1662: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1663: }
1664: }
1665: PetscFree(values);
1666: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1667: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1668: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1669: }
1670: PetscFree(cols);
1671: return(0);
1672: }
1676: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J)
1677: {
1678: PetscErrorCode ierr;
1679: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1680: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1681: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1682: MPI_Comm comm;
1683: PetscScalar *values;
1684: DMBoundaryType bx,by,bz;
1685: DMDAStencilType st;
1686: ISLocalToGlobalMapping ltog;
1689: /*
1690: nc - number of components per grid point
1691: col - number of colors needed in one direction for single component problem
1692: */
1693: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1694: col = 2*s + 1;
1696: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1697: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1698: PetscObjectGetComm((PetscObject)da,&comm);
1700: /* create the matrix */
1701: PetscMalloc1(col*col*col,&cols);
1703: DMGetLocalToGlobalMapping(da,<og);
1705: /* determine the matrix preallocation information */
1706: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1707: for (i=xs; i<xs+nx; i++) {
1708: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1709: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1710: for (j=ys; j<ys+ny; j++) {
1711: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1712: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1713: for (k=zs; k<zs+nz; k++) {
1714: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1715: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1717: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1719: /* Find block columns in block row */
1720: cnt = 0;
1721: for (ii=istart; ii<iend+1; ii++) {
1722: for (jj=jstart; jj<jend+1; jj++) {
1723: for (kk=kstart; kk<kend+1; kk++) {
1724: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1725: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1726: }
1727: }
1728: }
1729: }
1730: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1731: MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);
1732: }
1733: }
1734: }
1735: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1736: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1737: MatPreallocateFinalize(dnz,onz);
1739: MatSetLocalToGlobalMapping(J,ltog,ltog);
1741: /*
1742: For each node in the grid: we get the neighbors in the local (on processor ordering
1743: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1744: PETSc ordering.
1745: */
1746: if (!da->prealloc_only) {
1747: PetscCalloc1(col*col*col*nc*nc,&values);
1748: for (i=xs; i<xs+nx; i++) {
1749: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1750: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1751: for (j=ys; j<ys+ny; j++) {
1752: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1753: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1754: for (k=zs; k<zs+nz; k++) {
1755: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1756: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1758: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1760: cnt = 0;
1761: for (ii=istart; ii<iend+1; ii++) {
1762: for (jj=jstart; jj<jend+1; jj++) {
1763: for (kk=kstart; kk<kend+1; kk++) {
1764: if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) {
1765: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1766: }
1767: }
1768: }
1769: }
1770: L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);
1771: MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1772: }
1773: }
1774: }
1775: PetscFree(values);
1776: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1777: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1778: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1779: }
1780: PetscFree(cols);
1781: return(0);
1782: }
1784: /* ---------------------------------------------------------------------------------*/
1788: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J)
1789: {
1790: PetscErrorCode ierr;
1791: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1792: PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,l,p,*dnz,*onz;
1793: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1794: DM_DA *dd = (DM_DA*)da->data;
1795: PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill;
1796: MPI_Comm comm;
1797: PetscScalar *values;
1798: DMBoundaryType bx,by,bz;
1799: ISLocalToGlobalMapping ltog;
1800: DMDAStencilType st;
1803: /*
1804: nc - number of components per grid point
1805: col - number of colors needed in one direction for single component problem
1807: */
1808: DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);
1809: col = 2*s + 1;
1810: 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\
1811: by 2*stencil_width + 1\n");
1812: 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\
1813: by 2*stencil_width + 1\n");
1814: 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\
1815: by 2*stencil_width + 1\n");
1817: DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1818: DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1819: PetscObjectGetComm((PetscObject)da,&comm);
1821: PetscMalloc1(col*col*col*nc,&cols);
1822: DMGetLocalToGlobalMapping(da,<og);
1824: /* determine the matrix preallocation information */
1825: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1828: for (i=xs; i<xs+nx; i++) {
1829: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1830: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1831: for (j=ys; j<ys+ny; j++) {
1832: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1833: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1834: for (k=zs; k<zs+nz; k++) {
1835: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1836: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1838: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1840: for (l=0; l<nc; l++) {
1841: cnt = 0;
1842: for (ii=istart; ii<iend+1; ii++) {
1843: for (jj=jstart; jj<jend+1; jj++) {
1844: for (kk=kstart; kk<kend+1; kk++) {
1845: if (ii || jj || kk) {
1846: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1847: 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);
1848: }
1849: } else {
1850: if (dfill) {
1851: 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);
1852: } else {
1853: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1854: }
1855: }
1856: }
1857: }
1858: }
1859: row = l + nc*(slot);
1860: MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);
1861: }
1862: }
1863: }
1864: }
1865: MatSeqAIJSetPreallocation(J,0,dnz);
1866: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1867: MatPreallocateFinalize(dnz,onz);
1868: MatSetLocalToGlobalMapping(J,ltog,ltog);
1870: /*
1871: For each node in the grid: we get the neighbors in the local (on processor ordering
1872: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1873: PETSc ordering.
1874: */
1875: if (!da->prealloc_only) {
1876: PetscCalloc1(col*col*col*nc*nc*nc,&values);
1877: for (i=xs; i<xs+nx; i++) {
1878: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i));
1879: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1));
1880: for (j=ys; j<ys+ny; j++) {
1881: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j));
1882: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1));
1883: for (k=zs; k<zs+nz; k++) {
1884: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k));
1885: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1));
1887: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1889: for (l=0; l<nc; l++) {
1890: cnt = 0;
1891: for (ii=istart; ii<iend+1; ii++) {
1892: for (jj=jstart; jj<jend+1; jj++) {
1893: for (kk=kstart; kk<kend+1; kk++) {
1894: if (ii || jj || kk) {
1895: if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1896: 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);
1897: }
1898: } else {
1899: if (dfill) {
1900: 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);
1901: } else {
1902: for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1903: }
1904: }
1905: }
1906: }
1907: }
1908: row = l + nc*(slot);
1909: MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);
1910: }
1911: }
1912: }
1913: }
1914: PetscFree(values);
1915: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1916: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1917: MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1918: }
1919: PetscFree(cols);
1920: return(0);
1921: }