Actual source code: dadd.c
1: #include <petsc/private/dmdaimpl.h>
3: /*@
4: DMDACreatePatchIS - Creates an index set corresponding to a patch of the DA.
6: Collective
8: Input Parameters:
9: + da - the DMDA
10: . lower - a matstencil with i, j and k corresponding to the lower corner of the patch
11: . upper - a matstencil with i, j and k corresponding to the upper corner of the patch
12: - offproc - indicate whether the returned IS will contain off process indices
14: Output Parameters:
15: . is - the IS corresponding to the patch
17: Level: developer
19: Notes:
20: This routine always returns an IS on the DMDA's comm, if offproc is set to PETSC_TRUE,
21: the routine returns an IS with all the indices requested regardless of whether these indices
22: are present on the requesting rank or not. Thus, it is upon the caller to ensure that
23: the indices returned in this mode are appropriate. If offproc is set to PETSC_FALSE,
24: the IS only returns the subset of indices that are present on the requesting rank and there
25: is no duplication of indices.
27: .seealso: DMDACreateDomainDecomposition(), DMDACreateDomainDecompositionScatters()
28: @*/
29: PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is, PetscBool offproc)
30: {
31: PetscInt ms=0,ns=0,ps=0;
32: PetscInt mw=0,nw=0,pw=0;
33: PetscInt me=1,ne=1,pe=1;
34: PetscInt mr=0,nr=0,pr=0;
35: PetscInt ii,jj,kk;
36: PetscInt si,sj,sk;
37: PetscInt i,j,k,l,idx=0;
38: PetscInt base;
39: PetscInt xm=1,ym=1,zm=1;
40: PetscInt ox,oy,oz;
41: PetscInt m,n,p,M,N,P,dof;
42: const PetscInt *lx,*ly,*lz;
43: PetscInt nindices;
44: PetscInt *indices;
45: DM_DA *dd = (DM_DA*)da->data;
46: PetscBool skip_i=PETSC_TRUE, skip_j=PETSC_TRUE, skip_k=PETSC_TRUE;
47: PetscBool valid_j=PETSC_FALSE, valid_k=PETSC_FALSE; /* DMDA has at least 1 dimension */
51: M = dd->M; N = dd->N; P = dd->P;
52: m = dd->m; n = dd->n; p = dd->p;
53: dof = dd->w;
55: nindices = -1;
56: if (PetscLikely(upper->i - lower->i)) {
57: nindices = nindices*(upper->i - lower->i);
58: skip_i=PETSC_FALSE;
59: }
60: if (N>1) {
61: valid_j = PETSC_TRUE;
62: if (PetscLikely(upper->j - lower->j)) {
63: nindices = nindices*(upper->j - lower->j);
64: skip_j=PETSC_FALSE;
65: }
66: }
67: if (P>1) {
68: valid_k = PETSC_TRUE;
69: if (PetscLikely(upper->k - lower->k)) {
70: nindices = nindices*(upper->k - lower->k);
71: skip_k=PETSC_FALSE;
72: }
73: }
74: if (PetscLikely(nindices<0)) {
75: if (PetscUnlikely(skip_i && skip_j && skip_k)) {
76: nindices = 0;
77: } else nindices = nindices*(-1);
78: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Lower and Upper stencils are identical! Please check inputs.");
80: PetscMalloc1(nindices*dof,&indices);
81: DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);
83: if (!valid_k) {
84: k = 0;
85: upper->k=0;
86: lower->k=0;
87: }
88: if (!valid_j) {
89: j = 0;
90: upper->j=0;
91: lower->j=0;
92: }
94: if (offproc) {
95: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
96: /* start at index 0 on processor 0 */
97: mr = 0;
98: nr = 0;
99: pr = 0;
100: ms = 0;
101: ns = 0;
102: ps = 0;
103: if (lx) me = lx[0];
104: if (ly) ne = ly[0];
105: if (lz) pe = lz[0];
106: /*
107: If no indices are to be returned, create an empty is,
108: this prevents hanging in while loops
109: */
110: if (skip_i && skip_j && skip_k) goto createis;
111: /*
112: do..while loops to ensure the block gets entered once,
113: regardless of control condition being met, necessary for
114: cases when a subset of skip_i/j/k is true
115: */
116: if (skip_k) k = upper->k-oz; else k = lower->k-oz;
117: do {
118: if (skip_j) j = upper->j-oy; else j = lower->j-oy;
119: do {
120: if (skip_i) i = upper->i-ox; else i = lower->i-ox;
121: do {
122: /* "actual" indices rather than ones outside of the domain */
123: ii = i;
124: jj = j;
125: kk = k;
126: if (ii < 0) ii = M + ii;
127: if (jj < 0) jj = N + jj;
128: if (kk < 0) kk = P + kk;
129: if (ii > M-1) ii = ii - M;
130: if (jj > N-1) jj = jj - N;
131: if (kk > P-1) kk = kk - P;
132: /* gone out of processor range on x axis */
133: while (ii > me-1 || ii < ms) {
134: if (mr == m-1) {
135: ms = 0;
136: me = lx[0];
137: mr = 0;
138: } else {
139: mr++;
140: ms = me;
141: me += lx[mr];
142: }
143: }
144: /* gone out of processor range on y axis */
145: while (jj > ne-1 || jj < ns) {
146: if (nr == n-1) {
147: ns = 0;
148: ne = ly[0];
149: nr = 0;
150: } else {
151: nr++;
152: ns = ne;
153: ne += ly[nr];
154: }
155: }
156: /* gone out of processor range on z axis */
157: while (kk > pe-1 || kk < ps) {
158: if (pr == p-1) {
159: ps = 0;
160: pe = lz[0];
161: pr = 0;
162: } else {
163: pr++;
164: ps = pe;
165: pe += lz[pr];
166: }
167: }
168: /* compute the vector base on owning processor */
169: xm = me - ms;
170: ym = ne - ns;
171: zm = pe - ps;
172: base = ms*ym*zm + ns*M*zm + ps*M*N;
173: /* compute the local coordinates on owning processor */
174: si = ii - ms;
175: sj = jj - ns;
176: sk = kk - ps;
177: for (l=0;l<dof;l++) {
178: indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
179: idx++;
180: }
181: i++;
182: } while (i<upper->i-ox);
183: j++;
184: } while (j<upper->j-oy);
185: k++;
186: } while (k<upper->k-oz);
187: }
189: if (!offproc) {
190: DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw);
191: me = ms + mw;
192: if (N>1) ne = ns + nw;
193: if (P>1) pe = ps + pw;
194: /* Account for DM offsets */
195: ms = ms - ox; me = me - ox;
196: ns = ns - oy; ne = ne - oy;
197: ps = ps - oz; pe = pe - oz;
199: /* compute the vector base on owning processor */
200: xm = me - ms;
201: ym = ne - ns;
202: zm = pe - ps;
203: base = ms*ym*zm + ns*M*zm + ps*M*N;
204: /*
205: if no indices are to be returned, create an empty is,
206: this prevents hanging in while loops
207: */
208: if (skip_i && skip_j && skip_k) goto createis;
209: /*
210: do..while loops to ensure the block gets entered once,
211: regardless of control condition being met, necessary for
212: cases when a subset of skip_i/j/k is true
213: */
214: if (skip_k) k = upper->k-oz; else k = lower->k-oz;
215: do {
216: if (skip_j) j = upper->j-oy; else j = lower->j-oy;
217: do {
218: if (skip_i) i = upper->i-ox; else i = lower->i-ox;
219: do {
220: if (k>=ps && k<=pe-1) {
221: if (j>=ns && j<=ne-1) {
222: if (i>=ms && i<=me-1) {
223: /* compute the local coordinates on owning processor */
224: si = i - ms;
225: sj = j - ns;
226: sk = k - ps;
227: for (l=0; l<dof; l++) {
228: indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
229: idx++;
230: }
231: }
232: }
233: }
234: i++;
235: } while (i<upper->i-ox);
236: j++;
237: } while (j<upper->j-oy);
238: k++;
239: } while (k<upper->k-oz);
241: PetscRealloc((size_t)(idx*sizeof(PetscInt)), (void*)&indices);
242: }
244: createis:
245: ISCreateGeneral(PetscObjectComm((PetscObject)da),idx,indices,PETSC_OWN_POINTER,is);
246: return(0);
247: }
249: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
250: {
251: DM *da;
252: PetscInt dim,size,i,j,k,idx;
254: DMDALocalInfo info;
255: PetscInt xsize,ysize,zsize;
256: PetscInt xo,yo,zo;
257: PetscInt xs,ys,zs;
258: PetscInt xm=1,ym=1,zm=1;
259: PetscInt xol,yol,zol;
260: PetscInt m=1,n=1,p=1;
261: PetscInt M,N,P;
262: PetscInt pm,mtmp;
265: DMDAGetLocalInfo(dm,&info);
266: DMDAGetOverlap(dm,&xol,&yol,&zol);
267: DMDAGetNumLocalSubDomains(dm,&size);
268: PetscMalloc1(size,&da);
270: if (nlocal) *nlocal = size;
272: dim = info.dim;
274: M = info.xm;
275: N = info.ym;
276: P = info.zm;
278: if (dim == 1) {
279: m = size;
280: } else if (dim == 2) {
281: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
282: while (m > 0) {
283: n = size/m;
284: if (m*n*p == size) break;
285: m--;
286: }
287: } else if (dim == 3) {
288: n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); if (!n) n = 1;
289: while (n > 0) {
290: pm = size/n;
291: if (n*pm == size) break;
292: n--;
293: }
294: if (!n) n = 1;
295: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
296: if (!m) m = 1;
297: while (m > 0) {
298: p = size/(m*n);
299: if (m*n*p == size) break;
300: m--;
301: }
302: if (M > P && m < p) {mtmp = m; m = p; p = mtmp;}
303: }
305: zs = info.zs;
306: idx = 0;
307: for (k = 0; k < p; k++) {
308: ys = info.ys;
309: for (j = 0; j < n; j++) {
310: xs = info.xs;
311: for (i = 0; i < m; i++) {
312: if (dim == 1) {
313: xm = M/m + ((M % m) > i);
314: } else if (dim == 2) {
315: xm = M/m + ((M % m) > i);
316: ym = N/n + ((N % n) > j);
317: } else if (dim == 3) {
318: xm = M/m + ((M % m) > i);
319: ym = N/n + ((N % n) > j);
320: zm = P/p + ((P % p) > k);
321: }
323: xsize = xm;
324: ysize = ym;
325: zsize = zm;
326: xo = xs;
327: yo = ys;
328: zo = zs;
330: DMDACreate(PETSC_COMM_SELF,&(da[idx]));
331: DMSetOptionsPrefix(da[idx],"sub_");
332: DMSetDimension(da[idx], info.dim);
333: DMDASetDof(da[idx], info.dof);
335: DMDASetStencilType(da[idx],info.st);
336: DMDASetStencilWidth(da[idx],info.sw);
338: if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
339: xsize += xol;
340: xo -= xol;
341: }
342: if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
343: ysize += yol;
344: yo -= yol;
345: }
346: if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
347: zsize += zol;
348: zo -= zol;
349: }
351: if (info.bx == DM_BOUNDARY_PERIODIC || (xs+xm != info.mx)) xsize += xol;
352: if (info.by == DM_BOUNDARY_PERIODIC || (ys+ym != info.my)) ysize += yol;
353: if (info.bz == DM_BOUNDARY_PERIODIC || (zs+zm != info.mz)) zsize += zol;
355: if (info.bx != DM_BOUNDARY_PERIODIC) {
356: if (xo < 0) {
357: xsize += xo;
358: xo = 0;
359: }
360: if (xo+xsize > info.mx-1) {
361: xsize -= xo+xsize - info.mx;
362: }
363: }
364: if (info.by != DM_BOUNDARY_PERIODIC) {
365: if (yo < 0) {
366: ysize += yo;
367: yo = 0;
368: }
369: if (yo+ysize > info.my-1) {
370: ysize -= yo+ysize - info.my;
371: }
372: }
373: if (info.bz != DM_BOUNDARY_PERIODIC) {
374: if (zo < 0) {
375: zsize += zo;
376: zo = 0;
377: }
378: if (zo+zsize > info.mz-1) {
379: zsize -= zo+zsize - info.mz;
380: }
381: }
383: DMDASetSizes(da[idx], xsize, ysize, zsize);
384: DMDASetNumProcs(da[idx], 1, 1, 1);
385: DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);
387: /* set up as a block instead */
388: DMSetUp(da[idx]);
390: /* nonoverlapping region */
391: DMDASetNonOverlappingRegion(da[idx],xs,ys,zs,xm,ym,zm);
393: /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
394: DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);
395: xs += xm;
396: idx++;
397: }
398: ys += ym;
399: }
400: zs += zm;
401: }
402: *sdm = da;
403: return(0);
404: }
406: /*
407: Fills the local vector problem on the subdomain from the global problem.
409: Right now this assumes one subdomain per processor.
411: */
412: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat)
413: {
415: DMDALocalInfo info,subinfo;
416: DM subdm;
417: MatStencil upper,lower;
418: IS idis,isis,odis,osis,gdis;
419: Vec svec,dvec,slvec;
420: PetscInt xm,ym,zm,xs,ys,zs;
421: PetscInt i;
422: PetscBool patchis_offproc = PETSC_TRUE;
425: /* allocate the arrays of scatters */
426: if (iscat) {PetscMalloc1(nsubdms,iscat);}
427: if (oscat) {PetscMalloc1(nsubdms,oscat);}
428: if (lscat) {PetscMalloc1(nsubdms,lscat);}
430: DMDAGetLocalInfo(dm,&info);
431: for (i = 0; i < nsubdms; i++) {
432: subdm = subdms[i];
433: DMDAGetLocalInfo(subdm,&subinfo);
434: DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);
436: /* create the global and subdomain index sets for the inner domain */
437: lower.i = xs;
438: lower.j = ys;
439: lower.k = zs;
440: upper.i = xs+xm;
441: upper.j = ys+ym;
442: upper.k = zs+zm;
443: DMDACreatePatchIS(dm,&lower,&upper,&idis,patchis_offproc);
444: DMDACreatePatchIS(subdm,&lower,&upper,&isis,patchis_offproc);
446: /* create the global and subdomain index sets for the outer subdomain */
447: lower.i = subinfo.xs;
448: lower.j = subinfo.ys;
449: lower.k = subinfo.zs;
450: upper.i = subinfo.xs+subinfo.xm;
451: upper.j = subinfo.ys+subinfo.ym;
452: upper.k = subinfo.zs+subinfo.zm;
453: DMDACreatePatchIS(dm,&lower,&upper,&odis,patchis_offproc);
454: DMDACreatePatchIS(subdm,&lower,&upper,&osis,patchis_offproc);
456: /* global and subdomain ISes for the local indices of the subdomain */
457: /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
458: lower.i = subinfo.gxs;
459: lower.j = subinfo.gys;
460: lower.k = subinfo.gzs;
461: upper.i = subinfo.gxs+subinfo.gxm;
462: upper.j = subinfo.gys+subinfo.gym;
463: upper.k = subinfo.gzs+subinfo.gzm;
464: DMDACreatePatchIS(dm,&lower,&upper,&gdis,patchis_offproc);
466: /* form the scatter */
467: DMGetGlobalVector(dm,&dvec);
468: DMGetGlobalVector(subdm,&svec);
469: DMGetLocalVector(subdm,&slvec);
471: if (iscat) {VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[i]);}
472: if (oscat) {VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[i]);}
473: if (lscat) {VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[i]);}
475: DMRestoreGlobalVector(dm,&dvec);
476: DMRestoreGlobalVector(subdm,&svec);
477: DMRestoreLocalVector(subdm,&slvec);
479: ISDestroy(&idis);
480: ISDestroy(&isis);
482: ISDestroy(&odis);
483: ISDestroy(&osis);
485: ISDestroy(&gdis);
486: }
487: return(0);
488: }
490: PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois)
491: {
493: PetscInt i;
494: DMDALocalInfo info,subinfo;
495: MatStencil lower,upper;
496: PetscBool patchis_offproc = PETSC_TRUE;
499: DMDAGetLocalInfo(dm,&info);
500: if (iis) {PetscMalloc1(n,iis);}
501: if (ois) {PetscMalloc1(n,ois);}
503: for (i = 0;i < n; i++) {
504: DMDAGetLocalInfo(subdm[i],&subinfo);
505: if (iis) {
506: /* create the inner IS */
507: lower.i = info.xs;
508: lower.j = info.ys;
509: lower.k = info.zs;
510: upper.i = info.xs+info.xm;
511: upper.j = info.ys+info.ym;
512: upper.k = info.zs+info.zm;
513: DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i],patchis_offproc);
514: }
516: if (ois) {
517: /* create the outer IS */
518: lower.i = subinfo.xs;
519: lower.j = subinfo.ys;
520: lower.k = subinfo.zs;
521: upper.i = subinfo.xs+subinfo.xm;
522: upper.j = subinfo.ys+subinfo.ym;
523: upper.k = subinfo.zs+subinfo.zm;
524: DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i],patchis_offproc);
525: }
526: }
527: return(0);
528: }
530: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm)
531: {
533: DM *sdm;
534: PetscInt n,i;
537: DMDASubDomainDA_Private(dm,&n,&sdm);
538: if (names) {
539: PetscMalloc1(n,names);
540: for (i=0;i<n;i++) (*names)[i] = NULL;
541: }
542: DMDASubDomainIS_Private(dm,n,sdm,iis,ois);
543: if (subdm) *subdm = sdm;
544: else {
545: for (i=0;i<n;i++) {
546: DMDestroy(&sdm[i]);
547: }
548: }
549: if (len) *len = n;
550: return(0);
551: }