xref: /petsc/src/dm/impls/da/dadd.c (revision ff9846d91d616d1fb8b6a5c35db72d596e376598)
1 #include <petsc-private/daimpl.h>
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "DMDACreatePatchIS"
5 
6 PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is)
7 {
8   PetscErrorCode ierr;
9   PetscInt       i,j,k,idx;
10   PetscInt       ii,jj,kk;
11   Vec            v;
12   PetscInt       n,pn,bs;
13   PetscMPIInt    rank;
14   PetscSF        sf,psf;
15   PetscLayout    map;
16   MPI_Comm       comm;
17   PetscInt       *natidx,*globidx,*leafidx;
18   PetscInt       *pnatidx,*pleafidx;
19   PetscInt       base;
20   PetscInt       ox,oy,oz;
21   DM_DA          *dd;
22   PetscFunctionBegin;
23 
24   comm = ((PetscObject)da)->comm;
25   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
26   dd = (DM_DA *)da->data;
27 
28   /* construct the natural mapping */
29   ierr = DMGetGlobalVector(da,&v);CHKERRQ(ierr);
30   ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
31   ierr = VecGetOwnershipRange(v,&base,PETSC_NULL);CHKERRQ(ierr);
32   ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
33   ierr = DMRestoreGlobalVector(da,&v);CHKERRQ(ierr);
34 
35   /* construct the layout */
36   ierr = PetscLayoutCreate(comm,&map);CHKERRQ(ierr);
37   ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr);
38   ierr = PetscLayoutSetLocalSize(map,n);CHKERRQ(ierr);
39   ierr = PetscLayoutSetUp(map);CHKERRQ(ierr);
40 
41   /* construct the list of natural indices on this process when PETSc ordering is considered */
42   ierr = DMDAGetOffset(da,&ox,&oy,&oz);CHKERRQ(ierr);
43   ierr = PetscMalloc(sizeof(PetscInt)*n,&natidx);CHKERRQ(ierr);
44   ierr = PetscMalloc(sizeof(PetscInt)*n,&globidx);CHKERRQ(ierr);
45   ierr = PetscMalloc(sizeof(PetscInt)*n,&leafidx);CHKERRQ(ierr);
46   idx = 0;
47   for (k=dd->zs;k<dd->ze;k++) {
48     for (j=dd->ys;j<dd->ye;j++) {
49       for (i=dd->xs;i<dd->xe;i++) {
50         natidx[idx] = i + dd->w*(j*dd->M + k*dd->M*dd->N);
51         globidx[idx] = base + idx;
52         leafidx[idx] = 0;
53         idx++;
54       }
55     }
56   }
57 
58   if (idx != n) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE, "for some reason the count is wrong.");
59 
60   /* construct the SF going from the natural indices to the local set of PETSc indices */
61   ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr);
62   ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr);
63   ierr = PetscSFSetGraphLayout(sf,map,n,PETSC_NULL,PETSC_OWN_POINTER,natidx);CHKERRQ(ierr);
64 
65   /* broadcast the global indices over to the corresponding natural indices */
66   ierr = PetscSFGatherBegin(sf,MPIU_INT,globidx,leafidx);CHKERRQ(ierr);
67   ierr = PetscSFGatherEnd(sf,MPIU_INT,globidx,leafidx);CHKERRQ(ierr);
68   ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
69 
70 
71   pn = dd->w*(upper->k - lower->k)*(upper->j - lower->j)*(upper->i - lower->i);
72   ierr = PetscMalloc(sizeof(PetscInt)*pn,&pnatidx);CHKERRQ(ierr);
73   ierr = PetscMalloc(sizeof(PetscInt)*pn,&pleafidx);CHKERRQ(ierr);
74   idx = 0;
75   for (k=lower->k-oz;k<upper->k-oz;k++) {
76     for (j=lower->j-oy;j<upper->j-oy;j++) {
77       for (i=dd->w*(lower->i-ox);i<dd->w*(upper->i-ox);i++) {
78         ii = i % (dd->w*dd->M);
79         jj = j % dd->N;
80         kk = k % dd->P;
81         if (ii < 0) ii = dd->w*dd->M + ii;
82         if (jj < 0) jj = dd->N + jj;
83         if (kk < 0) kk = dd->P + kk;
84         pnatidx[idx] = ii + dd->w*(jj*dd->M + kk*dd->M*dd->N);
85         idx++;
86       }
87     }
88   }
89 
90   if (idx != pn) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE, "for some reason the count is wrong");
91 
92   ierr = PetscSFCreate(comm,&psf);CHKERRQ(ierr);
93   ierr = PetscSFSetFromOptions(psf);CHKERRQ(ierr);
94   ierr = PetscSFSetGraphLayout(psf,map,pn,PETSC_NULL,PETSC_OWN_POINTER,pnatidx);CHKERRQ(ierr);
95 
96   /* broadcast the global indices through to the patch */
97   ierr = PetscSFBcastBegin(psf,MPIU_INT,leafidx,pleafidx);CHKERRQ(ierr);
98   ierr = PetscSFBcastEnd(psf,MPIU_INT,leafidx,pleafidx);CHKERRQ(ierr);
99 
100   /* create the IS */
101   ierr = ISCreateGeneral(comm,pn,pleafidx,PETSC_OWN_POINTER,is);CHKERRQ(ierr);
102 
103   ierr = PetscSFDestroy(&psf);CHKERRQ(ierr);
104 
105   ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr);
106 
107   ierr = PetscFree(globidx);CHKERRQ(ierr);
108   ierr = PetscFree(leafidx);CHKERRQ(ierr);
109   ierr = PetscFree(natidx);CHKERRQ(ierr);
110   ierr = PetscFree(pnatidx);CHKERRQ(ierr);
111 
112   PetscFunctionReturn(0);
113 }
114 
115 #undef __FUNCT__
116 #define __FUNCT__ "DMDASubDomainDA_Private"
117 PetscErrorCode DMDASubDomainDA_Private(DM dm, DM *dddm) {
118   DM               da;
119   DM_DA            *dd;
120   PetscErrorCode   ierr;
121   DMDALocalInfo    info;
122   PetscReal        lmin[3],lmax[3];
123   PetscInt         xsize,ysize,zsize;
124   PetscInt         xo,yo,zo;
125 
126   PetscFunctionBegin;
127   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
128   ierr = DMDACreate(PETSC_COMM_SELF,&da);CHKERRQ(ierr);
129   ierr = DMSetOptionsPrefix(da,"sub_");CHKERRQ(ierr);
130 
131   ierr = DMDASetDim(da, info.dim);CHKERRQ(ierr);
132   ierr = DMDASetDof(da, info.dof);CHKERRQ(ierr);
133 
134   ierr = DMDASetStencilType(da,info.st);CHKERRQ(ierr);
135   ierr = DMDASetStencilWidth(da,info.sw);CHKERRQ(ierr);
136 
137   dd = (DM_DA *)dm->data;
138 
139   /* local with overlap */
140   xsize = info.xm;
141   ysize = info.ym;
142   zsize = info.zm;
143   xo    = info.xs;
144   yo    = info.ys;
145   zo    = info.zs;
146   if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs != 0)) {
147     xsize += dd->overlap;
148     xo -= dd->overlap;
149   }
150   if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys != 0)) {
151     ysize += dd->overlap;
152     yo    -= dd->overlap;
153   }
154   if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs != 0)) {
155     zsize += dd->overlap;
156     zo    -= dd->overlap;
157   }
158 
159   if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs+info.xm != info.mx)) {
160     xsize += dd->overlap;
161   }
162   if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys+info.ym != info.my)) {
163     ysize += dd->overlap;
164   }
165   if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs+info.zm != info.mz)) {
166     zsize += dd->overlap;
167   }
168 
169   ierr = DMDASetSizes(da, xsize, ysize, zsize);CHKERRQ(ierr);
170   ierr = DMDASetNumProcs(da, 1, 1, 1);CHKERRQ(ierr);
171   ierr = DMDASetBoundaryType(da, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED);CHKERRQ(ierr);
172 
173   /* set up as a block instead */
174   ierr = DMSetUp(da);CHKERRQ(ierr);
175   ierr = DMDASetOffset(da,xo,yo,zo);CHKERRQ(ierr);
176 
177 
178   /* todo - nonuniform coordinates */
179   ierr = DMDAGetLocalBoundingBox(dm,lmin,lmax);CHKERRQ(ierr);
180   ierr = DMDASetUniformCoordinates(da,lmin[0],lmax[0],lmin[1],lmax[1],lmin[2],lmax[2]);CHKERRQ(ierr);
181 
182   dd = (DM_DA *)da->data;
183   dd->Mo = info.mx;
184   dd->No = info.my;
185   dd->Po = info.mz;
186 
187   *dddm = da;
188 
189   PetscFunctionReturn(0);
190 }
191 
192 #undef __FUNCT__
193 #define __FUNCT__ "DMCreateDomainDecompositionScatters_DA"
194 /*
195  Fills the local vector problem on the subdomain from the global problem.
196 
197  */
198 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) {
199   PetscErrorCode   ierr;
200   DMDALocalInfo    dinfo,sinfo;
201   IS               isis,idis,osis,odis,gsis,gdis;
202   PetscInt         *ididx,*isidx,*odidx,*osidx,*gdidx,*gsidx,*idx_global,n_global,*idx_sub,n_sub;
203   PetscInt         l,i,j,k,d,n_i,n_o,n_g,sl,dl,di,dj,dk,si,sj,sk;
204   Vec              dvec,svec,slvec;
205   DM               subdm;
206 
207   PetscFunctionBegin;
208 
209   /* allocate the arrays of scatters */
210   if (iscat) {ierr = PetscMalloc(sizeof(VecScatter *),iscat);CHKERRQ(ierr);}
211   if (oscat) {ierr = PetscMalloc(sizeof(VecScatter *),oscat);CHKERRQ(ierr);}
212   if (lscat) {ierr = PetscMalloc(sizeof(VecScatter *),lscat);CHKERRQ(ierr);}
213 
214   ierr = DMDAGetLocalInfo(dm,&dinfo);CHKERRQ(ierr);
215   ierr = DMDAGetGlobalIndices(dm,&n_global,&idx_global);CHKERRQ(ierr);
216   for (l = 0;l < nsubdms;l++) {
217     n_i = 0;
218     n_o = 0;
219     n_g = 0;
220     subdm = subdms[l];
221     ierr = DMDAGetLocalInfo(subdm,&sinfo);CHKERRQ(ierr);
222     ierr = DMDAGetGlobalIndices(subdm,&n_sub,&idx_sub);CHKERRQ(ierr);
223     /* count the three region sizes */
224     for (k=sinfo.gzs;k<sinfo.gzs+sinfo.gzm;k++) {
225       for (j=sinfo.gys;j<sinfo.gys+sinfo.gym;j++) {
226         for (i=sinfo.gxs;i<sinfo.gxs+sinfo.gxm;i++) {
227           for (d=0;d<sinfo.dof;d++) {
228             if (k >= sinfo.zs           && j >= sinfo.ys         && i >= sinfo.xs &&
229                 k <  sinfo.zs+sinfo.zm  && j < sinfo.ys+sinfo.ym && i < sinfo.xs+sinfo.xm) {
230 
231               /* interior - subinterior overlap */
232               if (k >= dinfo.zs            && j >= dinfo.ys          && i >= dinfo.xs &&
233                   k <  dinfo.zs+dinfo.zm  && j < dinfo.ys+dinfo.ym && i < dinfo.xs+dinfo.xm) {
234                 n_i++;
235               }
236               /* ghost - subinterior overlap */
237               if (k >= dinfo.gzs            && j >= dinfo.gys          && i >= dinfo.gxs &&
238                   k <  dinfo.gzs+dinfo.gzm  && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) {
239                 n_o++;
240               }
241             }
242 
243             /* ghost - subghost overlap */
244             if (k >= dinfo.gzs            && j >= dinfo.gys          && i >= dinfo.gxs &&
245                 k <  dinfo.gzs+dinfo.gzm  && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) {
246               n_g++;
247             }
248           }
249         }
250       }
251     }
252 
253     if (n_g == 0) SETERRQ(((PetscObject)subdm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Processor-local domain and subdomain do not intersect!");
254 
255     /* local and subdomain local index set indices */
256     ierr = PetscMalloc(sizeof(PetscInt)*n_i,&ididx);CHKERRQ(ierr);
257     ierr = PetscMalloc(sizeof(PetscInt)*n_i,&isidx);CHKERRQ(ierr);
258 
259     ierr = PetscMalloc(sizeof(PetscInt)*n_o,&odidx);CHKERRQ(ierr);
260     ierr = PetscMalloc(sizeof(PetscInt)*n_o,&osidx);CHKERRQ(ierr);
261 
262     ierr = PetscMalloc(sizeof(PetscInt)*n_g,&gdidx);CHKERRQ(ierr);
263     ierr = PetscMalloc(sizeof(PetscInt)*n_g,&gsidx);CHKERRQ(ierr);
264 
265     n_i = 0; n_o = 0;n_g = 0;
266     for (k=sinfo.gzs;k<sinfo.gzs+sinfo.gzm;k++) {
267       for (j=sinfo.gys;j<sinfo.gys+sinfo.gym;j++) {
268         for (i=sinfo.gxs;i<sinfo.gxs+sinfo.gxm;i++) {
269           for (d=0;d<sinfo.dof;d++) {
270             si = i - sinfo.gxs;
271             sj = j - sinfo.gys;
272             sk = k - sinfo.gzs;
273             sl = d + sinfo.dof*(si + sinfo.gxm*(sj + sinfo.gym*sk));
274             di = i - dinfo.gxs;
275             dj = j - dinfo.gys;
276             dk = k - dinfo.gzs;
277             dl = d + dinfo.dof*(di + dinfo.gxm*(dj + dinfo.gym*dk));
278 
279             if (k >= sinfo.zs           && j >= sinfo.ys         && i >= sinfo.xs &&
280                 k <  sinfo.zs+sinfo.zm  && j < sinfo.ys+sinfo.ym && i < sinfo.xs+sinfo.xm) {
281 
282               /* interior - subinterior overlap */
283               if (k >= dinfo.zs            && j >= dinfo.ys          && i >= dinfo.xs &&
284                   k <  dinfo.zs+dinfo.zm  && j < dinfo.ys+dinfo.ym && i < dinfo.xs+dinfo.xm) {
285                 ididx[n_i] = idx_global[dl];
286                 isidx[n_i] = idx_sub[sl];
287                 n_i++;
288               }
289               /* ghost - subinterior overlap */
290               if (k >= dinfo.gzs            && j >= dinfo.gys          && i >= dinfo.gxs &&
291                   k <  dinfo.gzs+dinfo.gzm  && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) {
292                 odidx[n_o] = idx_global[dl];
293                 osidx[n_o] = idx_sub[sl];
294                 n_o++;
295               }
296             }
297 
298             /* ghost - subghost overlap */
299             if (k >= dinfo.gzs            && j >= dinfo.gys          && i >= dinfo.gxs &&
300                 k <  dinfo.gzs+dinfo.gzm  && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) {
301               gdidx[n_g] = idx_global[dl];
302               gsidx[n_g] = sl;
303               n_g++;
304             }
305           }
306         }
307       }
308     }
309 
310     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_i,ididx,PETSC_OWN_POINTER,&idis);CHKERRQ(ierr);
311     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_i,isidx,PETSC_OWN_POINTER,&isis);CHKERRQ(ierr);
312 
313     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_o,odidx,PETSC_OWN_POINTER,&odis);CHKERRQ(ierr);
314     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_o,osidx,PETSC_OWN_POINTER,&osis);CHKERRQ(ierr);
315 
316     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_g,gdidx,PETSC_OWN_POINTER,&gdis);CHKERRQ(ierr);
317     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_g,gsidx,PETSC_OWN_POINTER,&gsis);CHKERRQ(ierr);
318 
319     /* form the scatter */
320     ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr);
321     ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr);
322     ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr);
323 
324     if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[l]);CHKERRQ(ierr);}
325     if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[l]);CHKERRQ(ierr);}
326     if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,gsis,&(*lscat)[l]);CHKERRQ(ierr);}
327 
328     ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr);
329     ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr);
330     ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr);
331 
332     ierr = ISDestroy(&idis);CHKERRQ(ierr);
333     ierr = ISDestroy(&isis);CHKERRQ(ierr);
334 
335     ierr = ISDestroy(&odis);CHKERRQ(ierr);
336     ierr = ISDestroy(&osis);CHKERRQ(ierr);
337 
338     ierr = ISDestroy(&gdis);CHKERRQ(ierr);
339     ierr = ISDestroy(&gsis);CHKERRQ(ierr);
340   }
341   PetscFunctionReturn(0);
342 
343 }
344 
345 #undef __FUNCT__
346 #define __FUNCT__ "DMDASubDomainIS_Private"
347 /* We have that the interior regions are going to be the same, but the ghost regions might not match up
348 
349 ----------
350 ----------
351 --++++++o=
352 --++++++o=
353 --++++++o=
354 --++++++o=
355 --ooooooo=
356 --========
357 
358 Therefore, for each point in the overall, we must check if it's:
359 
360 1. +: In the interior of the global dm; it lines up
361 2. o: In the overlap region -- for now the same as 1; no overlap
362 3. =: In the shared ghost region -- handled by DMCreateDomainDecompositionLocalScatter()
363 4. -: In the nonshared ghost region
364  */
365 
366 PetscErrorCode DMDASubDomainIS_Private(DM dm,DM subdm,IS *iis,IS *ois) {
367   PetscErrorCode   ierr;
368   DMDALocalInfo    info,subinfo;
369   PetscInt         *iiidx,*oiidx,*gidx,gindx;
370   PetscInt         i,j,k,d,n,nsub,nover,llindx,lindx,li,lj,lk,gi,gj,gk;
371 
372   PetscFunctionBegin;
373   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
374   ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr);
375   ierr = DMDAGetGlobalIndices(dm,&n,&gidx);CHKERRQ(ierr);
376   /* todo -- overlap */
377   nsub = info.xm*info.ym*info.zm*info.dof;
378   nover = subinfo.xm*subinfo.ym*subinfo.zm*subinfo.dof;
379   /* iis is going to have size of the local problem's global part but have a lot of fill-in */
380   ierr = PetscMalloc(sizeof(PetscInt)*nsub,&iiidx);CHKERRQ(ierr);
381   /* ois is going to have size of the local problem's global part */
382   ierr = PetscMalloc(sizeof(PetscInt)*nover,&oiidx);CHKERRQ(ierr);
383   /* loop over the ghost region of the subdm and fill in the indices */
384   for (k=subinfo.gzs;k<subinfo.gzs+subinfo.gzm;k++) {
385     for (j=subinfo.gys;j<subinfo.gys+subinfo.gym;j++) {
386       for (i=subinfo.gxs;i<subinfo.gxs+subinfo.gxm;i++) {
387         for (d=0;d<subinfo.dof;d++) {
388           li = i - subinfo.xs;
389           lj = j - subinfo.ys;
390           lk = k - subinfo.zs;
391           lindx = d + subinfo.dof*(li + subinfo.xm*(lj + subinfo.ym*lk));
392           li = i - info.xs;
393           lj = j - info.ys;
394           lk = k - info.zs;
395           llindx = d + info.dof*(li + info.xm*(lj + info.ym*lk));
396           gi = i - info.gxs;
397           gj = j - info.gys;
398           gk = k - info.gzs;
399           gindx = d + info.dof*(gi + info.gxm*(gj + info.gym*gk));
400 
401           /* check if the current point is inside the interior region */
402           if (k >= info.zs          && j >= info.ys          && i >= info.xs &&
403               k <  info.zs+info.zm  && j < info.ys+info.ym   && i < info.xs+info.xm) {
404             iiidx[llindx] = gidx[gindx];
405             oiidx[lindx] = gidx[gindx];
406             /* overlap region */
407           } else if (k >= subinfo.zs             && j >= subinfo.ys                && i >= subinfo.xs &&
408                      k <  subinfo.zs+subinfo.zm  && j < subinfo.ys+subinfo.ym   && i < subinfo.xs+subinfo.xm) {
409             oiidx[lindx] = gidx[gindx];
410           }
411         }
412       }
413     }
414   }
415 
416   /* create the index sets */
417   ierr = ISCreateGeneral(PETSC_COMM_SELF,nsub,iiidx,PETSC_OWN_POINTER,iis);CHKERRQ(ierr);
418   ierr = ISCreateGeneral(PETSC_COMM_SELF,nover,oiidx,PETSC_OWN_POINTER,ois);CHKERRQ(ierr);
419   PetscFunctionReturn(0);
420 }
421 
422 #undef __FUNCT__
423 #define __FUNCT__ "DMCreateDomainDecomposition_DA"
424 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) {
425   PetscErrorCode ierr;
426   IS             iis0,ois0;
427   DM             subdm0;
428   PetscFunctionBegin;
429   if (len)*len = 1;
430 
431   if (iis) {ierr = PetscMalloc(sizeof(IS *),iis);CHKERRQ(ierr);}
432   if (ois) {ierr = PetscMalloc(sizeof(IS *),ois);CHKERRQ(ierr);}
433   if (subdm) {ierr = PetscMalloc(sizeof(DM *),subdm);CHKERRQ(ierr);}
434   if (names) {ierr = PetscMalloc(sizeof(char *),names);CHKERRQ(ierr);}
435   ierr = DMDASubDomainDA_Private(dm,&subdm0);CHKERRQ(ierr);
436   ierr = DMDASubDomainIS_Private(dm,subdm0,&iis0,&ois0);CHKERRQ(ierr);
437   if (iis) {
438     (*iis)[0] = iis0;
439   } else {
440     ierr = ISDestroy(&iis0);CHKERRQ(ierr);
441   }
442   if (ois) {
443     (*ois)[0] = ois0;
444   } else {
445     ierr = ISDestroy(&ois0);CHKERRQ(ierr);
446   }
447   if (subdm) (*subdm)[0] = subdm0;
448   if (names) (*names)[0] = 0;
449   PetscFunctionReturn(0);
450 }
451