Actual source code: fdda.c
1: /*$Id: fdda.c,v 1.75 2001/08/07 21:31:51 bsmith Exp $*/
2:
3: #include src/dm/da/daimpl.h
4: #include petscmat.h
7: EXTERN int DAGetColoring1d_MPIAIJ(DA,ISColoringType,ISColoring *);
8: EXTERN int DAGetColoring2d_MPIAIJ(DA,ISColoringType,ISColoring *);
9: EXTERN int DAGetColoring2d_5pt_MPIAIJ(DA,ISColoringType,ISColoring *);
10: EXTERN int DAGetColoring3d_MPIAIJ(DA,ISColoringType,ISColoring *);
12: #undef __FUNCT__
14: /*@C
15: DAGetColoring - Gets the coloring required for computing the Jacobian via
16: finite differences on a function defined using a stencil on the DA.
18: Collective on DA
20: Input Parameter:
21: + da - the distributed array
22: - ctype - IS_COLORING_LOCAL or IS_COLORING_GHOSTED
24: Output Parameters:
25: . coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)
27: Level: advanced
29: Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used
30: for efficient (parallel or thread based) triangular solves etc is NOT yet
31: available.
34: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring
36: @*/
37: int DAGetColoring(DA da,ISColoringType ctype,ISColoring *coloring)
38: {
39: int ierr,dim;
42: /*
43: m
44: ------------------------------------------------------
45: | |
46: | |
47: | ---------------------- |
48: | | | |
49: n | yn | | |
50: | | | |
51: | .--------------------- |
52: | (xs,ys) xn |
53: | . |
54: | (gxs,gys) |
55: | |
56: -----------------------------------------------------
57: */
59: /*
60: nc - number of components per grid point
61: col - number of colors needed in one direction for single component problem
62:
63: */
64: DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
66: /*
67: We do not provide a getcoloring function in the DA operations because
68: the basic DA does not know about matrices. We think of DA as being more
69: more low-level then matrices.
70: */
71: if (dim == 1) {
72: DAGetColoring1d_MPIAIJ(da,ctype,coloring);
73: } else if (dim == 2) {
74: DAGetColoring2d_MPIAIJ(da,ctype,coloring);
75: } else if (dim == 3) {
76: DAGetColoring3d_MPIAIJ(da,ctype,coloring);
77: } else {
78: SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
79: }
80: return(0);
81: }
83: /* ---------------------------------------------------------------------------------*/
85: #undef __FUNCT__
87: int DAGetColoring2d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
88: {
89: int ierr,xs,ys,nx,ny,*colors,i,j,ii,gxs,gys,gnx,gny;
90: int m,n,M,N,dim,w,s,k,nc,col,size;
91: MPI_Comm comm;
92: DAPeriodicType wrap;
93: DAStencilType st;
96: /*
97: nc - number of components per grid point
98: col - number of colors needed in one direction for single component problem
99:
100: */
101: DAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&w,&s,&wrap,&st);
102: nc = w;
103: col = 2*s + 1;
104: if (DAXPeriodic(wrap) && (m % col)){
105: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
106: by 2*stencil_width + 1n");
107: }
108: if (DAYPeriodic(wrap) && (n % col)){
109: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
110: by 2*stencil_width + 1n");
111: }
112: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
113: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
114: PetscObjectGetComm((PetscObject)da,&comm);
115: MPI_Comm_size(comm,&size);
117: /* create the coloring */
118: if (st == DA_STENCIL_STAR && s == 1 && !DAXPeriodic(wrap) && !DAYPeriodic(wrap)) {
119: DAGetColoring2d_5pt_MPIAIJ(da,ctype,coloring);
120: } else if (ctype == IS_COLORING_LOCAL) {
121: if (!da->localcoloring) {
122: PetscMalloc(nc*nx*ny*sizeof(int),&colors);
123: ii = 0;
124: for (j=ys; j<ys+ny; j++) {
125: for (i=xs; i<xs+nx; i++) {
126: for (k=0; k<nc; k++) {
127: colors[ii++] = k + nc*((i % col) + col*(j % col));
128: }
129: }
130: }
131: ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
132: }
133: *coloring = da->localcoloring;
134: } else if (ctype == IS_COLORING_GHOSTED) {
135: if (!da->ghostedcoloring) {
136: PetscMalloc(nc*gnx*gny*sizeof(int),&colors);
137: ii = 0;
138: for (j=gys; j<gys+gny; j++) {
139: for (i=gxs; i<gxs+gnx; i++) {
140: for (k=0; k<nc; k++) {
141: /* the complicated stuff is to handle periodic boundaries */
142: colors[ii++] = k + nc*( (((i < 0) ? m+i:((i >= m) ? i-m:i)) % col) +
143: col*(((j < 0) ? n+j:((j >= n) ? j-n:j)) % col));
144: }
145: }
146: }
147: ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
148: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
149: }
150: *coloring = da->ghostedcoloring;
151: } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
152: ISColoringReference(*coloring);
153: return(0);
154: }
156: /* ---------------------------------------------------------------------------------*/
158: #undef __FUNCT__
160: int DAGetColoring3d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
161: {
162: int ierr,xs,ys,nx,ny,*colors,i,j,gxs,gys,gnx,gny;
163: int m,n,p,dim,s,k,nc,col,size,zs,gzs,ii,l,nz,gnz,M,N,P;
164: MPI_Comm comm;
165: DAPeriodicType wrap;
166: DAStencilType st;
169: /*
170: nc - number of components per grid point
171: col - number of colors needed in one direction for single component problem
172:
173: */
174: DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&wrap,&st);
175: col = 2*s + 1;
176: if (DAXPeriodic(wrap) && (m % col)){
177: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
178: by 2*stencil_width + 1n");
179: }
180: if (DAYPeriodic(wrap) && (n % col)){
181: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
182: by 2*stencil_width + 1n");
183: }
184: if (DAZPeriodic(wrap) && (p % col)){
185: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisiblen
186: by 2*stencil_width + 1n");
187: }
189: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
190: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
191: PetscObjectGetComm((PetscObject)da,&comm);
192: MPI_Comm_size(comm,&size);
194: /* create the coloring */
195: if (ctype == IS_COLORING_LOCAL) {
196: if (!da->localcoloring) {
197: PetscMalloc(nc*nx*ny*nz*sizeof(int),&colors);
198: ii = 0;
199: for (k=zs; k<zs+nz; k++) {
200: for (j=ys; j<ys+ny; j++) {
201: for (i=xs; i<xs+nx; i++) {
202: for (l=0; l<nc; l++) {
203: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
204: }
205: }
206: }
207: }
208: ISColoringCreate(comm,nc*nx*ny*nz,colors,&da->localcoloring);
209: }
210: *coloring = da->localcoloring;
211: } else if (ctype == IS_COLORING_GHOSTED) {
212: if (!da->ghostedcoloring) {
213: PetscMalloc(nc*gnx*gny*gnz*sizeof(int),&colors);
214: ii = 0;
215: for (k=gzs; k<gzs+gnz; k++) {
216: for (j=gys; j<gys+gny; j++) {
217: for (i=gxs; i<gxs+gnx; i++) {
218: for (l=0; l<nc; l++) {
219: /* the complicated stuff is to handle periodic boundaries */
220: colors[ii++] = l + nc*( (((i < 0) ? m+i:((i >= m) ? i-m:i)) % col) +
221: col*(((j < 0) ? n+j:((j >= n) ? j-n:j)) % col) +
222: col*col*(((k < 0) ? p+k:((k >= p) ? k-p:k)) % col));
223: }
224: }
225: }
226: }
227: ISColoringCreate(comm,nc*gnx*gny*gnz,colors,&da->ghostedcoloring);
228: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
229: }
230: *coloring = da->ghostedcoloring;
231: } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
232: ISColoringReference(*coloring);
233: return(0);
234: }
236: /* ---------------------------------------------------------------------------------*/
238: #undef __FUNCT__
240: int DAGetColoring1d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
241: {
242: int ierr,xs,nx,*colors,i,i1,gxs,gnx,l;
243: int m,M,dim,s,nc,col,size;
244: MPI_Comm comm;
245: DAPeriodicType wrap;
248: /*
249: nc - number of components per grid point
250: col - number of colors needed in one direction for single component problem
251:
252: */
253: DAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&wrap,0);
254: col = 2*s + 1;
256: if (DAXPeriodic(wrap) && (m % col)) {
257: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisiblen
258: by 2*stencil_width + 1n");
259: }
261: DAGetCorners(da,&xs,0,0,&nx,0,0);
262: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
263: PetscObjectGetComm((PetscObject)da,&comm);
264: MPI_Comm_size(comm,&size);
266: /* create the coloring */
267: if (ctype == IS_COLORING_LOCAL) {
268: if (!da->localcoloring) {
269: PetscMalloc(nc*nx*sizeof(int),&colors);
270: i1 = 0;
271: for (i=xs; i<xs+nx; i++) {
272: for (l=0; l<nc; l++) {
273: colors[i1++] = l + nc*(i % col);
274: }
275: }
276: ISColoringCreate(comm,nc*nx,colors,&da->localcoloring);
277: }
278: *coloring = da->localcoloring;
279: } else if (ctype == IS_COLORING_GHOSTED) {
280: if (!da->ghostedcoloring) {
281: PetscMalloc(nc*gnx*sizeof(int),&colors);
282: i1 = 0;
283: for (i=gxs; i<gxs+gnx; i++) {
284: for (l=0; l<nc; l++) {
285: /* the complicated stuff is to handle periodic boundaries */
286: colors[i1++] = l + nc*((((i < 0) ? m+i:((i >= m) ? i-m:i)) % col));
287: }
288: }
289: ISColoringCreate(comm,nc*gnx,colors,&da->ghostedcoloring);
290: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
291: }
292: *coloring = da->ghostedcoloring;
293: } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
294: ISColoringReference(*coloring);
295: return(0);
296: }
298: #undef __FUNCT__
300: int DAGetColoring2d_5pt_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
301: {
302: int ierr,xs,ys,nx,ny,*colors,i,j,ii,gxs,gys,gnx,gny;
303: int m,n,dim,w,s,k,nc;
304: MPI_Comm comm;
307: /*
308: nc - number of components per grid point
309: col - number of colors needed in one direction for single component problem
310:
311: */
312: ierr = DAGetInfo(da,&dim,&m,&n,0,0,0,0,&w,&s,0,0);
313: nc = w;
314: ierr = DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
315: ierr = DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
316: ierr = PetscObjectGetComm((PetscObject)da,&comm);
318: /* create the coloring */
319: if (ctype == IS_COLORING_LOCAL) {
320: if (!da->localcoloring) {
321: PetscMalloc(nc*nx*ny*sizeof(int),&colors);
322: ii = 0;
323: for (j=ys; j<ys+ny; j++) {
324: for (i=xs; i<xs+nx; i++) {
325: for (k=0; k<nc; k++) {
326: colors[ii++] = k + nc*((3*j+i) % 5);
327: }
328: }
329: }
330: ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
331: }
332: *coloring = da->localcoloring;
333: } else if (ctype == IS_COLORING_GHOSTED) {
334: if (!da->ghostedcoloring) {
335: PetscMalloc(nc*gnx*gny*sizeof(int),&colors);
336: ii = 0;
337: for (j=gys; j<gys+gny; j++) {
338: for (i=gxs; i<gxs+gnx; i++) {
339: for (k=0; k<nc; k++) {
340: colors[ii++] = k + nc*((3*j+i) % 5);
341: }
342: }
343: }
344: ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
345: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
346: }
347: *coloring = da->ghostedcoloring;
348: } else SETERRQ1(1,"Unknown ISColoringType %d",ctype);
349: return(0);
350: }
352: /* =========================================================================== */
353: EXTERN int DAGetMatrix1d_MPIAIJ(DA,Mat *);
354: EXTERN int DAGetMatrix2d_MPIAIJ(DA,Mat *);
355: EXTERN int DAGetMatrix3d_MPIAIJ(DA,Mat *);
356: EXTERN int DAGetMatrix3d_MPIBAIJ(DA,Mat *);
357: EXTERN int DAGetMatrix3d_MPISBAIJ(DA,Mat *);
359: #undef __FUNCT__
361: /*@C
362: DAGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for
363: computing the Jacobian on a function defined using the stencil set in the DA.
365: Collective on DA
367: Input Parameter:
368: + da - the distributed array
369: - mtype - either MATMPIAIJ or MATMPIBAIJ
371: Output Parameters:
372: . J - matrix with the correct nonzero structure
373: (obviously without the correct Jacobian values)
375: Level: advanced
377: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
378: do not need to do it yourself.
380: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate()
382: @*/
383: int DAGetMatrix(DA da,MatType mtype,Mat *J)
384: {
385: int ierr,dim;
386: PetscTruth aij,baij,sbaij;
389: /*
390: m
391: ------------------------------------------------------
392: | |
393: | |
394: | ---------------------- |
395: | | | |
396: n | yn | | |
397: | | | |
398: | .--------------------- |
399: | (xs,ys) xn |
400: | . |
401: | (gxs,gys) |
402: | |
403: -----------------------------------------------------
404: */
406: /*
407: nc - number of components per grid point
408: col - number of colors needed in one direction for single component problem
409:
410: */
411: DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
413: /*
414: We do not provide a getcoloring function in the DA operations because
415: the basic DA does not know about matrices. We think of DA as being more
416: more low-level then matrices.
417: */
418: PetscStrcmp(MATMPIAIJ,mtype,&aij);
419: PetscStrcmp(MATMPIBAIJ,mtype,&baij);
420: PetscStrcmp(MATMPISBAIJ,mtype,&sbaij);
421: if (aij) {
422: if (dim == 1) {
423: DAGetMatrix1d_MPIAIJ(da,J);
424: } else if (dim == 2) {
425: DAGetMatrix2d_MPIAIJ(da,J);
426: } else if (dim == 3) {
427: DAGetMatrix3d_MPIAIJ(da,J);
428: }
429: } else if(baij) {
430: if (dim == 3) {
431: DAGetMatrix3d_MPIBAIJ(da,J);
432: } else {
433: SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
434: }
435: } else if(sbaij) {
436: if (dim == 3) {
437: DAGetMatrix3d_MPISBAIJ(da,J);
438: } else {
439: SETERRQ1(1,"Not done for %d dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
440: }
441: } else {
442: SETERRQ1(1,"Not done for %s matrix type, send us mail petsc-maint@mcs.anl.gov for code",mtype);
443: }
444: return(0);
445: }
447: /* ---------------------------------------------------------------------------------*/
449: #undef __FUNCT__
451: int DAGetMatrix2d_MPIAIJ(DA da,Mat *J)
452: {
453: int ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
454: int m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
455: int lstart,lend,pstart,pend,*dnz,*onz,size;
456: int dims[2],starts[2];
457: MPI_Comm comm;
458: PetscScalar *values;
459: DAPeriodicType wrap;
460: ISLocalToGlobalMapping ltog;
461: DAStencilType st;
464: /*
465: nc - number of components per grid point
466: col - number of colors needed in one direction for single component problem
467:
468: */
469: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
470: col = 2*s + 1;
471: if (DAXPeriodic(wrap) && (m % col)){
472: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
473: by 2*stencil_width + 1n");
474: }
475: if (DAYPeriodic(wrap) && (n % col)){
476: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
477: by 2*stencil_width + 1n");
478: }
479: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
480: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
481: PetscObjectGetComm((PetscObject)da,&comm);
482: MPI_Comm_size(comm,&size);
484: /* create empty Jacobian matrix */
485: MatCreate(comm,nc*nx*ny,nc*nx*ny,PETSC_DECIDE,PETSC_DECIDE,J);
487: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
488: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
489: PetscMalloc(nc*sizeof(int),&rows);
490: PetscMalloc(col*col*nc*nc*sizeof(int),&cols);
491: DAGetISLocalToGlobalMapping(da,<og);
492:
493: /* determine the matrix preallocation information */
494: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
495: for (i=xs; i<xs+nx; i++) {
497: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
498: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
500: for (j=ys; j<ys+ny; j++) {
501: slot = i - gxs + gnx*(j - gys);
503: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
504: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
506: cnt = 0;
507: for (k=0; k<nc; k++) {
508: for (l=lstart; l<lend+1; l++) {
509: for (p=pstart; p<pend+1; p++) {
510: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
511: cols[cnt++] = k + nc*(slot + gnx*l + p);
512: }
513: }
514: }
515: rows[k] = k + nc*(slot);
516: }
517: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
518: }
519: }
520: /* set matrix type and preallocation information */
521: if (size > 1) {
522: MatSetType(*J,MATMPIAIJ);
523: } else {
524: MatSetType(*J,MATSEQAIJ);
525: }
526: MatSeqAIJSetPreallocation(*J,0,dnz);
527: MatSeqBAIJSetPreallocation(*J,nc,0,dnz);
528: MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
529: MatMPIBAIJSetPreallocation(*J,nc,0,dnz,0,onz);
530: MatPreallocateFinalize(dnz,onz);
531: MatSetLocalToGlobalMapping(*J,ltog);
532: DAGetGhostCorners(da,&starts[0],&starts[1],PETSC_IGNORE,&dims[0],&dims[1],PETSC_IGNORE);
533: MatSetStencil(*J,2,dims,starts,nc);
535: /*
536: For each node in the grid: we get the neighbors in the local (on processor ordering
537: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
538: PETSc ordering.
539: */
540: for (i=xs; i<xs+nx; i++) {
541:
542: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
543: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
544:
545: for (j=ys; j<ys+ny; j++) {
546: slot = i - gxs + gnx*(j - gys);
547:
548: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
549: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
551: cnt = 0;
552: for (k=0; k<nc; k++) {
553: for (l=lstart; l<lend+1; l++) {
554: for (p=pstart; p<pend+1; p++) {
555: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
556: cols[cnt++] = k + nc*(slot + gnx*l + p);
557: }
558: }
559: }
560: rows[k] = k + nc*(slot);
561: }
562: MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
563: }
564: }
565: PetscFree(values);
566: PetscFree(rows);
567: PetscFree(cols);
568: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
569: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
570: return(0);
571: }
573: /* ---------------------------------------------------------------------------------*/
575: #undef __FUNCT__
577: int DAGetMatrix3d_MPIAIJ(DA da,Mat *J)
578: {
579: int ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
580: int m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
581: int istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,size;
582: int dims[3],starts[3];
583: MPI_Comm comm;
584: PetscScalar *values;
585: DAPeriodicType wrap;
586: ISLocalToGlobalMapping ltog;
587: DAStencilType st;
590: /*
591: nc - number of components per grid point
592: col - number of colors needed in one direction for single component problem
593:
594: */
595: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
596: col = 2*s + 1;
597: if (DAXPeriodic(wrap) && (m % col)){
598: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisiblen
599: by 2*stencil_width + 1n");
600: }
601: if (DAYPeriodic(wrap) && (n % col)){
602: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisiblen
603: by 2*stencil_width + 1n");
604: }
605: if (DAZPeriodic(wrap) && (p % col)){
606: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisiblen
607: by 2*stencil_width + 1n");
608: }
610: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
611: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
612: PetscObjectGetComm((PetscObject)da,&comm);
613: MPI_Comm_size(comm,&size);
616: /* create the matrix */
617: /* create empty Jacobian matrix */
618: MatCreate(comm,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,J);
619: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
620: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
621: PetscMalloc(nc*sizeof(int),&rows);
622: PetscMalloc(col*col*col*nc*sizeof(int),&cols);
623: DAGetISLocalToGlobalMapping(da,<og);
625: /* determine the matrix preallocation information */
626: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
627: for (i=xs; i<xs+nx; i++) {
628: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
629: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
630: for (j=ys; j<ys+ny; j++) {
631: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
632: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
633: for (k=zs; k<zs+nz; k++) {
634: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
635: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
636:
637: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
638:
639: cnt = 0;
640: for (l=0; l<nc; l++) {
641: for (ii=istart; ii<iend+1; ii++) {
642: for (jj=jstart; jj<jend+1; jj++) {
643: for (kk=kstart; kk<kend+1; kk++) {
644: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
645: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
646: }
647: }
648: }
649: }
650: rows[l] = l + nc*(slot);
651: }
652: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
653: }
654: }
655: }
656: /* set matrix type and preallocation */
657: if (size > 1) {
658: MatSetType(*J,MATMPIAIJ);
659: } else {
660: MatSetType(*J,MATSEQAIJ);
661: }
662: MatSeqAIJSetPreallocation(*J,0,dnz);
663: MatSeqBAIJSetPreallocation(*J,nc,0,dnz);
664: MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
665: MatMPIBAIJSetPreallocation(*J,nc,0,dnz,0,onz);
666: MatPreallocateFinalize(dnz,onz);
667: MatSetLocalToGlobalMapping(*J,ltog);
668: DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
669: MatSetStencil(*J,3,dims,starts,nc);
671: /*
672: For each node in the grid: we get the neighbors in the local (on processor ordering
673: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
674: PETSc ordering.
675: */
676: for (i=xs; i<xs+nx; i++) {
677: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
678: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
679: for (j=ys; j<ys+ny; j++) {
680: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
681: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
682: for (k=zs; k<zs+nz; k++) {
683: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
684: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
685:
686: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
687:
688: cnt = 0;
689: for (l=0; l<nc; l++) {
690: for (ii=istart; ii<iend+1; ii++) {
691: for (jj=jstart; jj<jend+1; jj++) {
692: for (kk=kstart; kk<kend+1; kk++) {
693: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
694: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
695: }
696: }
697: }
698: }
699: rows[l] = l + nc*(slot);
700: }
701: MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
702: }
703: }
704: }
705: PetscFree(values);
706: PetscFree(rows);
707: PetscFree(cols);
708: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
709: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
710: return(0);
711: }
713: /* ---------------------------------------------------------------------------------*/
715: #undef __FUNCT__
717: int DAGetMatrix1d_MPIAIJ(DA da,Mat *J)
718: {
719: int ierr,xs,nx,i,i1,slot,gxs,gnx;
720: int m,dim,s,*cols,nc,*rows,col,cnt,l;
721: int istart,iend,size;
722: int dims[1],starts[1];
723: MPI_Comm comm;
724: PetscScalar *values;
725: DAPeriodicType wrap;
726: ISLocalToGlobalMapping ltog;
729: /*
730: nc - number of components per grid point
731: col - number of colors needed in one direction for single component problem
732:
733: */
734: DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
735: col = 2*s + 1;
737: if (DAXPeriodic(wrap) && (m % col)) {
738: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisiblen
739: by 2*stencil_width + 1n");
740: }
743: DAGetCorners(da,&xs,0,0,&nx,0,0);
744: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
745: PetscObjectGetComm((PetscObject)da,&comm);
746: MPI_Comm_size(comm,&size);
748: /* create empty Jacobian matrix */
749:
750: ierr = MatCreate(comm,nc*nx,nc*nx,PETSC_DECIDE,PETSC_DECIDE,J);
751: if (size > 1) {
752: MatSetType(*J,MATMPIAIJ);
753: } else {
754: MatSetType(*J,MATSEQAIJ);
755: }
756: MatSeqAIJSetPreallocation(*J,col*nc,0);
757: MatSeqBAIJSetPreallocation(*J,nc,col,0);
758: MatMPIAIJSetPreallocation(*J,col*nc,0,0,0);
759: MatMPIBAIJSetPreallocation(*J,nc,col,0,0,0);
760: DAGetGhostCorners(da,&starts[0],PETSC_IGNORE,PETSC_IGNORE,&dims[0],PETSC_IGNORE,PETSC_IGNORE);
761: MatSetStencil(*J,1,dims,starts,nc);
762:
763: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
764: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
765: PetscMalloc(nc*sizeof(int),&rows);
766: PetscMalloc(col*nc*sizeof(int),&cols);
767:
768: DAGetISLocalToGlobalMapping(da,<og);
769: MatSetLocalToGlobalMapping(*J,ltog);
770:
771: /*
772: For each node in the grid: we get the neighbors in the local (on processor ordering
773: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
774: PETSc ordering.
775: */
776: for (i=xs; i<xs+nx; i++) {
777: istart = PetscMax(-s,gxs - i);
778: iend = PetscMin(s,gxs + gnx - i - 1);
779: slot = i - gxs;
780:
781: cnt = 0;
782: for (l=0; l<nc; l++) {
783: for (i1=istart; i1<iend+1; i1++) {
784: cols[cnt++] = l + nc*(slot + i1);
785: }
786: rows[l] = l + nc*(slot);
787: }
788: MatSetValuesLocal(*J,nc,rows,cnt,cols,values,INSERT_VALUES);
789: }
790: PetscFree(values);
791: PetscFree(rows);
792: PetscFree(cols);
793: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
794: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
795: return(0);
796: }
798: #undef __FUNCT__
800: int DAGetMatrix3d_MPIBAIJ(DA da,Mat *J)
801: {
802: int ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
803: int m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
804: int istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
805: MPI_Comm comm;
806: PetscScalar *values;
807: DAPeriodicType wrap;
808: DAStencilType st;
809: ISLocalToGlobalMapping ltog;
812: /*
813: nc - number of components per grid point
814: col - number of colors needed in one direction for single component problem
815:
816: */
817: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
818: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
819: col = 2*s + 1;
821: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
822: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
823: PetscObjectGetComm((PetscObject)da,&comm);
825: /* create the matrix */
826: ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
827: ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
828: ierr = PetscMalloc(col*col*col*sizeof(int),&cols);
830: DAGetISLocalToGlobalMappingBlck(da,<og);
832: /* determine the matrix preallocation information */
833: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
834: for (i=xs; i<xs+nx; i++) {
835: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
836: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
837: for (j=ys; j<ys+ny; j++) {
838: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
839: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
840: for (k=zs; k<zs+nz; k++) {
841: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
842: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
844: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
846: /* Find block columns in block row */
847: cnt = 0;
848: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
849: for (ii=istart; ii<iend+1; ii++) {
850: for (jj=jstart; jj<jend+1; jj++) {
851: for (kk=kstart; kk<kend+1; kk++) {
852: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
853: }
854: }
855: }
856: } else { /* Star stencil */
857: cnt = 0;
858: for (ii=istart; ii<iend+1; ii++) {
859: if (ii) {
860: /* jj and kk must be zero */
861: /* cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; */
862: cols[cnt++] = slot + ii;
863: } else {
864: for (jj=jstart; jj<jend+1; jj++) {
865: if (jj) {
866: /* ii and kk must be zero */
867: cols[cnt++] = slot + gnx*jj;
868: } else {
869: /* ii and jj must be zero */
870: for (kk=kstart; kk<kend+1; kk++) {
871: cols[cnt++] = slot + gnx*gny*kk;
872: }
873: }
874: }
875: }
876: }
877: }
878: MatPreallocateSetLocal(ltog,1,&slot,cnt,cols,dnz,onz);
879: }
880: }
881: }
883: /* create empty Jacobian matrix */
884: MatCreateMPIBAIJ(comm,nc,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,0,dnz,0,onz,J);
886: MatPreallocateFinalize(dnz,onz);
887: MatSetLocalToGlobalMappingBlock(*J,ltog);
889: /*
890: For each node in the grid: we get the neighbors in the local (on processor ordering
891: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
892: PETSc ordering.
893: */
895: for (i=xs; i<xs+nx; i++) {
896: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
897: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
898: for (j=ys; j<ys+ny; j++) {
899: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
900: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
901: for (k=zs; k<zs+nz; k++) {
902: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
903: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
904:
905: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
906:
907: cnt = 0;
908: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
909: for (ii=istart; ii<iend+1; ii++) {
910: for (jj=jstart; jj<jend+1; jj++) {
911: for (kk=kstart; kk<kend+1; kk++) {
912: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
913: }
914: }
915: }
916: } else { /* Star stencil */
917: cnt = 0;
918: for (ii=istart; ii<iend+1; ii++) {
919: if (ii) {
920: /* jj and kk must be zero */
921: cols[cnt++] = slot + ii;
922: } else {
923: for (jj=jstart; jj<jend+1; jj++) {
924: if (jj) {
925: /* ii and kk must be zero */
926: cols[cnt++] = slot + gnx*jj;
927: } else {
928: /* ii and jj must be zero */
929: for (kk=kstart; kk<kend+1; kk++) {
930: cols[cnt++] = slot + gnx*gny*kk;
931: }
932: }
933: }
934: }
935: }
936: }
937: MatSetValuesBlockedLocal(*J,1,&slot,cnt,cols,values,INSERT_VALUES);
938: }
939: }
940: }
941: PetscFree(values);
942: PetscFree(cols);
943: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
944: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
945: return(0);
946: }
948: /* BAD! Almost identical to the BAIJ one */
949: #undef __FUNCT__
951: int DAGetMatrix3d_MPISBAIJ(DA da,Mat *J)
952: {
953: int ierr,xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
954: int m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
955: int istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
956: MPI_Comm comm;
957: PetscScalar *values;
958: DAPeriodicType wrap;
959: DAStencilType st;
960: ISLocalToGlobalMapping ltog;
963: /*
964: nc - number of components per grid point
965: col - number of colors needed in one direction for single component problem
966:
967: */
968: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
969: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
970: col = 2*s + 1;
972: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
973: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
974: PetscObjectGetComm((PetscObject)da,&comm);
976: /* create the matrix */
977: ierr = PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
978: ierr = PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
979: ierr = PetscMalloc(col*col*col*sizeof(int),&cols);
981: DAGetISLocalToGlobalMappingBlck(da,<og);
983: /* determine the matrix preallocation information */
984: MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
985: for (i=xs; i<xs+nx; i++) {
986: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
987: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
988: for (j=ys; j<ys+ny; j++) {
989: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
990: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
991: for (k=zs; k<zs+nz; k++) {
992: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
993: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
995: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
997: /* Find block columns in block row */
998: cnt = 0;
999: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1000: for (ii=istart; ii<iend+1; ii++) {
1001: for (jj=jstart; jj<jend+1; jj++) {
1002: for (kk=kstart; kk<kend+1; kk++) {
1003: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1004: }
1005: }
1006: }
1007: } else { /* Star stencil */
1008: cnt = 0;
1009: for (ii=istart; ii<iend+1; ii++) {
1010: if (ii) {
1011: /* jj and kk must be zero */
1012: /* cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; */
1013: cols[cnt++] = slot + ii;
1014: } else {
1015: for (jj=jstart; jj<jend+1; jj++) {
1016: if (jj) {
1017: /* ii and kk must be zero */
1018: cols[cnt++] = slot + gnx*jj;
1019: } else {
1020: /* ii and jj must be zero */
1021: for (kk=kstart; kk<kend+1; kk++) {
1022: cols[cnt++] = slot + gnx*gny*kk;
1023: }
1024: }
1025: }
1026: }
1027: }
1028: }
1029: MatPreallocateSymmetricSetLocal(ltog,1,&slot,cnt,cols,dnz,onz);
1030: }
1031: }
1032: }
1034: /* create empty Jacobian matrix */
1035: MatCreateMPISBAIJ(comm,nc,nc*nx*ny*nz,nc*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE,0,dnz,0,onz,J);
1037: MatPreallocateFinalize(dnz,onz);
1038: MatSetLocalToGlobalMappingBlock(*J,ltog);
1040: /*
1041: For each node in the grid: we get the neighbors in the local (on processor ordering
1042: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1043: PETSc ordering.
1044: */
1046: for (i=xs; i<xs+nx; i++) {
1047: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1048: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1049: for (j=ys; j<ys+ny; j++) {
1050: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1051: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1052: for (k=zs; k<zs+nz; k++) {
1053: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1054: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1055:
1056: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1057:
1058: cnt = 0;
1059: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1060: for (ii=istart; ii<iend+1; ii++) {
1061: for (jj=jstart; jj<jend+1; jj++) {
1062: for (kk=kstart; kk<kend+1; kk++) {
1063: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1064: }
1065: }
1066: }
1067: } else { /* Star stencil */
1068: cnt = 0;
1069: for (ii=istart; ii<iend+1; ii++) {
1070: if (ii) {
1071: /* jj and kk must be zero */
1072: cols[cnt++] = slot + ii;
1073: } else {
1074: for (jj=jstart; jj<jend+1; jj++) {
1075: if (jj) {
1076: /* ii and kk must be zero */
1077: cols[cnt++] = slot + gnx*jj;
1078: } else {
1079: /* ii and jj must be zero */
1080: for (kk=kstart; kk<kend+1; kk++) {
1081: cols[cnt++] = slot + gnx*gny*kk;
1082: }
1083: }
1084: }
1085: }
1086: }
1087: }
1088: MatSetValuesBlockedLocal(*J,1,&slot,cnt,cols,values,INSERT_VALUES);
1089: }
1090: }
1091: }
1092: PetscFree(values);
1093: PetscFree(cols);
1094: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1095: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1096: return(0);
1097: }