xref: /petsc/src/dm/impls/da/dalocal.c (revision 193ac0bc5128976c4ec2c1a67f3f9cb026b77f22)
1 
2 /*
3   Code for manipulating distributed regular arrays in parallel.
4 */
5 
6 #include <private/daimpl.h>    /*I   "petscdmda.h"   I*/
7 
8 /*
9    This allows the DMDA vectors to properly tell MATLAB their dimensions
10 */
11 #if defined(PETSC_HAVE_MATLAB_ENGINE)
12 #include <engine.h>   /* MATLAB include file */
13 #include <mex.h>      /* MATLAB include file */
14 EXTERN_C_BEGIN
15 #undef __FUNCT__
16 #define __FUNCT__ "VecMatlabEnginePut_DA2d"
17 PetscErrorCode  VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
18 {
19   PetscErrorCode ierr;
20   PetscInt       n,m;
21   Vec            vec = (Vec)obj;
22   PetscScalar    *array;
23   mxArray        *mat;
24   DM             da;
25 
26   PetscFunctionBegin;
27   ierr = PetscObjectQuery((PetscObject)vec,"DM",(PetscObject*)&da);CHKERRQ(ierr);
28   if (!da) SETERRQ(((PetscObject)vec)->comm,PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
29   ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr);
30 
31   ierr = VecGetArray(vec,&array);CHKERRQ(ierr);
32 #if !defined(PETSC_USE_COMPLEX)
33   mat  = mxCreateDoubleMatrix(m,n,mxREAL);
34 #else
35   mat  = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
36 #endif
37   ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr);
38   ierr = PetscObjectName(obj);CHKERRQ(ierr);
39   engPutVariable((Engine *)mengine,obj->name,mat);
40 
41   ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr);
42   PetscFunctionReturn(0);
43 }
44 EXTERN_C_END
45 #endif
46 
47 
48 #undef __FUNCT__
49 #define __FUNCT__ "DMCreateLocalVector_DA"
50 PetscErrorCode  DMCreateLocalVector_DA(DM da,Vec* g)
51 {
52   PetscErrorCode ierr;
53   DM_DA          *dd = (DM_DA*)da->data;
54 
55   PetscFunctionBegin;
56   PetscValidHeaderSpecific(da,DM_CLASSID,1);
57   PetscValidPointer(g,2);
58   ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr);
59   ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr);
60   ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr);
61   ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr);
62   ierr = PetscObjectCompose((PetscObject)*g,"DM",(PetscObject)da);CHKERRQ(ierr);
63 #if defined(PETSC_HAVE_MATLAB_ENGINE)
64   if (dd->w == 1  && dd->dim == 2) {
65     ierr = PetscObjectComposeFunctionDynamic((PetscObject)*g,"PetscMatlabEnginePut_C","VecMatlabEnginePut_DA2d",VecMatlabEnginePut_DA2d);CHKERRQ(ierr);
66   }
67 #endif
68   PetscFunctionReturn(0);
69 }
70 
71 #undef __FUNCT__
72 #define __FUNCT__ "DMGetLocalVector"
73 /*@
74    DMGetLocalVector - Gets a Seq PETSc vector that
75    may be used with the DMXXX routines. This vector has spaces for the ghost values.
76 
77    Not Collective
78 
79    Input Parameter:
80 .  dm - the distributed array
81 
82    Output Parameter:
83 .  g - the local vector
84 
85    Level: beginner
86 
87    Note:
88    The vector values are NOT initialized and may have garbage in them, so you may need
89    to zero them.
90 
91    The output parameter, g, is a regular PETSc vector that should be returned with
92    DMRestoreLocalVector() DO NOT call VecDestroy() on it.
93 
94    VecStride*() operations can be useful when using DM with dof > 1
95 
96 .keywords: distributed array, create, local, vector
97 
98 .seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
99           DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(),
100           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMCreateLocalVector(), DMRestoreLocalVector(),
101           VecStrideMax(), VecStrideMin(), VecStrideNorm()
102 @*/
103 PetscErrorCode  DMGetLocalVector(DM dm,Vec* g)
104 {
105   PetscErrorCode ierr,i;
106 
107   PetscFunctionBegin;
108   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
109   PetscValidPointer(g,2);
110   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
111     if (dm->localin[i]) {
112       *g             = dm->localin[i];
113       dm->localin[i] = PETSC_NULL;
114       goto alldone;
115     }
116   }
117   ierr = DMCreateLocalVector(dm,g);CHKERRQ(ierr);
118 
119   alldone:
120   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
121     if (!dm->localout[i]) {
122       dm->localout[i] = *g;
123       break;
124     }
125   }
126   PetscFunctionReturn(0);
127 }
128 
129 #undef __FUNCT__
130 #define __FUNCT__ "DMRestoreLocalVector"
131 /*@
132    DMRestoreLocalVector - Returns a Seq PETSc vector that
133      obtained from DMGetLocalVector(). Do not use with vector obtained via
134      DMCreateLocalVector().
135 
136    Not Collective
137 
138    Input Parameter:
139 +  dm - the distributed array
140 -  g - the local vector
141 
142    Level: beginner
143 
144 .keywords: distributed array, create, local, vector
145 
146 .seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
147           DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(),
148           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMCreateLocalVector(), DMGetLocalVector()
149 @*/
150 PetscErrorCode  DMRestoreLocalVector(DM dm,Vec* g)
151 {
152   PetscErrorCode ierr;
153   PetscInt       i,j;
154 
155   PetscFunctionBegin;
156   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
157   PetscValidPointer(g,2);
158   for (j=0; j<DM_MAX_WORK_VECTORS; j++) {
159     if (*g == dm->localout[j]) {
160       dm->localout[j] = PETSC_NULL;
161       for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
162         if (!dm->localin[i]) {
163           dm->localin[i] = *g;
164           goto alldone;
165         }
166       }
167     }
168   }
169   ierr = VecDestroy(g);CHKERRQ(ierr);
170   alldone:
171   PetscFunctionReturn(0);
172 }
173 
174 #undef __FUNCT__
175 #define __FUNCT__ "DMGetGlobalVector"
176 /*@
177    DMGetGlobalVector - Gets a MPI PETSc vector that
178    may be used with the DMXXX routines.
179 
180    Collective on DM
181 
182    Input Parameter:
183 .  dm - the distributed array
184 
185    Output Parameter:
186 .  g - the global vector
187 
188    Level: beginner
189 
190    Note:
191    The vector values are NOT initialized and may have garbage in them, so you may need
192    to zero them.
193 
194    The output parameter, g, is a regular PETSc vector that should be returned with
195    DMRestoreGlobalVector() DO NOT call VecDestroy() on it.
196 
197    VecStride*() operations can be useful when using DM with dof > 1
198 
199 .keywords: distributed array, create, Global, vector
200 
201 .seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
202           DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(),
203           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMCreateLocalVector(), DMRestoreLocalVector()
204           VecStrideMax(), VecStrideMin(), VecStrideNorm()
205 
206 @*/
207 PetscErrorCode  DMGetGlobalVector(DM dm,Vec* g)
208 {
209   PetscErrorCode ierr;
210   PetscInt       i;
211 
212   PetscFunctionBegin;
213   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
214   PetscValidPointer(g,2);
215   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
216     if (dm->globalin[i]) {
217       *g             = dm->globalin[i];
218       dm->globalin[i] = PETSC_NULL;
219       goto alldone;
220     }
221   }
222   ierr = DMCreateGlobalVector(dm,g);CHKERRQ(ierr);
223 
224   alldone:
225   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
226     if (!dm->globalout[i]) {
227       dm->globalout[i] = *g;
228       break;
229     }
230   }
231   PetscFunctionReturn(0);
232 }
233 
234 #undef __FUNCT__
235 #define __FUNCT__ "DMClearGlobalVectors"
236 /*@
237    DMClearGlobalVectors - Destroys all the global vectors that have been stashed in this DM
238 
239    Collective on DM
240 
241    Input Parameter:
242 .  dm - the distributed array
243 
244    Level: developer
245 
246 .keywords: distributed array, create, Global, vector
247 
248 .seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
249           DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(),
250           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMCreateLocalVector(), DMRestoreLocalVector()
251           VecStrideMax(), VecStrideMin(), VecStrideNorm()
252 
253 @*/
254 PetscErrorCode  DMClearGlobalVectors(DM dm)
255 {
256   PetscErrorCode ierr;
257   PetscInt       i;
258 
259   PetscFunctionBegin;
260   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
261   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
262     if (dm->globalout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Clearing DM of global vectors that has a global vector obtained with DMGetGlobalVector()");
263     ierr = VecDestroy(&dm->globalin[i]);CHKERRQ(ierr);
264   }
265   PetscFunctionReturn(0);
266 }
267 
268 #undef __FUNCT__
269 #define __FUNCT__ "DMRestoreGlobalVector"
270 /*@
271    DMRestoreGlobalVector - Returns a Seq PETSc vector that
272      obtained from DMGetGlobalVector(). Do not use with vector obtained via
273      DMCreateGlobalVector().
274 
275    Not Collective
276 
277    Input Parameter:
278 +  dm - the distributed array
279 -  g - the global vector
280 
281    Level: beginner
282 
283 .keywords: distributed array, create, global, vector
284 
285 .seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
286           DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToGlobalBegin(),
287           DMGlobalToGlobalEnd(), DMGlobalToGlobal(), DMCreateLocalVector(), DMGetGlobalVector()
288 @*/
289 PetscErrorCode  DMRestoreGlobalVector(DM dm,Vec* g)
290 {
291   PetscErrorCode ierr;
292   PetscInt       i,j;
293 
294   PetscFunctionBegin;
295   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
296   PetscValidPointer(g,2);
297   for (j=0; j<DM_MAX_WORK_VECTORS; j++) {
298     if (*g == dm->globalout[j]) {
299       dm->globalout[j] = PETSC_NULL;
300       for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
301         if (!dm->globalin[i]) {
302           dm->globalin[i] = *g;
303           goto alldone;
304         }
305       }
306     }
307   }
308   ierr = VecDestroy(g);CHKERRQ(ierr);
309   alldone:
310   PetscFunctionReturn(0);
311 }
312 
313 /* ------------------------------------------------------------------- */
314 #if defined(PETSC_HAVE_ADIC)
315 
316 EXTERN_C_BEGIN
317 #include <adic/ad_utils.h>
318 EXTERN_C_END
319 
320 #undef __FUNCT__
321 #define __FUNCT__ "DMDAGetAdicArray"
322 /*@C
323      DMDAGetAdicArray - Gets an array of derivative types for a DMDA
324 
325     Input Parameter:
326 +    da - information about my local patch
327 -    ghosted - do you want arrays for the ghosted or nonghosted patch
328 
329     Output Parameters:
330 +    vptr - array data structured to be passed to ad_FormFunctionLocal()
331 .    array_start - actual start of 1d array of all values that adiC can access directly (may be null)
332 -    tdof - total number of degrees of freedom represented in array_start (may be null)
333 
334      Notes:
335        The vector values are NOT initialized and may have garbage in them, so you may need
336        to zero them.
337 
338        Returns the same type of object as the DMDAVecGetArray() except its elements are
339            derivative types instead of PetscScalars
340 
341      Level: advanced
342 
343 .seealso: DMDARestoreAdicArray()
344 
345 @*/
346 PetscErrorCode  DMDAGetAdicArray(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
347 {
348   PetscErrorCode ierr;
349   PetscInt       j,i,deriv_type_size,xs,ys,xm,ym,zs,zm,itdof;
350   char           *iarray_start;
351   void           **iptr = (void**)vptr;
352   DM_DA          *dd = (DM_DA*)da->data;
353 
354   PetscFunctionBegin;
355   PetscValidHeaderSpecific(da,DM_CLASSID,1);
356   if (ghosted) {
357     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
358       if (dd->adarrayghostedin[i]) {
359         *iptr                   = dd->adarrayghostedin[i];
360         iarray_start            = (char*)dd->adstartghostedin[i];
361         itdof                   = dd->ghostedtdof;
362         dd->adarrayghostedin[i] = PETSC_NULL;
363         dd->adstartghostedin[i] = PETSC_NULL;
364 
365         goto done;
366       }
367     }
368     xs = dd->Xs;
369     ys = dd->Ys;
370     zs = dd->Zs;
371     xm = dd->Xe-dd->Xs;
372     ym = dd->Ye-dd->Ys;
373     zm = dd->Ze-dd->Zs;
374   } else {
375     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
376       if (dd->adarrayin[i]) {
377         *iptr            = dd->adarrayin[i];
378         iarray_start     = (char*)dd->adstartin[i];
379         itdof            = dd->tdof;
380         dd->adarrayin[i] = PETSC_NULL;
381         dd->adstartin[i] = PETSC_NULL;
382 
383         goto done;
384       }
385     }
386     xs = dd->xs;
387     ys = dd->ys;
388     zs = dd->zs;
389     xm = dd->xe-dd->xs;
390     ym = dd->ye-dd->ys;
391     zm = dd->ze-dd->zs;
392   }
393   deriv_type_size = PetscADGetDerivTypeSize();
394 
395   switch (dd->dim) {
396     case 1: {
397       void *ptr;
398       itdof = xm;
399 
400       ierr  = PetscMalloc(xm*deriv_type_size,&iarray_start);CHKERRQ(ierr);
401 
402       ptr   = (void*)(iarray_start - xs*deriv_type_size);
403       *iptr = (void*)ptr;
404       break;}
405     case 2: {
406       void **ptr;
407       itdof = xm*ym;
408 
409       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*deriv_type_size,&iarray_start);CHKERRQ(ierr);
410 
411       ptr  = (void**)(iarray_start + xm*ym*deriv_type_size - ys*sizeof(void*));
412       for(j=ys;j<ys+ym;j++) {
413         ptr[j] = iarray_start + deriv_type_size*(xm*(j-ys) - xs);
414       }
415       *iptr = (void*)ptr;
416       break;}
417     case 3: {
418       void ***ptr,**bptr;
419       itdof = xm*ym*zm;
420 
421       ierr  = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*deriv_type_size,&iarray_start);CHKERRQ(ierr);
422 
423       ptr  = (void***)(iarray_start + xm*ym*zm*deriv_type_size - zs*sizeof(void*));
424       bptr = (void**)(iarray_start + xm*ym*zm*deriv_type_size + zm*sizeof(void**));
425 
426       for(i=zs;i<zs+zm;i++) {
427         ptr[i] = bptr + ((i-zs)*ym - ys);
428       }
429       for (i=zs; i<zs+zm; i++) {
430         for (j=ys; j<ys+ym; j++) {
431           ptr[i][j] = iarray_start + deriv_type_size*(xm*ym*(i-zs) + xm*(j-ys) - xs);
432         }
433       }
434 
435       *iptr = (void*)ptr;
436       break;}
437     default:
438       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
439   }
440 
441   done:
442   /* add arrays to the checked out list */
443   if (ghosted) {
444     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
445       if (!dd->adarrayghostedout[i]) {
446         dd->adarrayghostedout[i] = *iptr ;
447         dd->adstartghostedout[i] = iarray_start;
448         dd->ghostedtdof          = itdof;
449         break;
450       }
451     }
452   } else {
453     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
454       if (!dd->adarrayout[i]) {
455         dd->adarrayout[i] = *iptr ;
456         dd->adstartout[i] = iarray_start;
457         dd->tdof          = itdof;
458         break;
459       }
460     }
461   }
462   if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Too many DMDA ADIC arrays obtained");
463   if (tdof)        *tdof = itdof;
464   if (array_start) *(void**)array_start = iarray_start;
465   PetscFunctionReturn(0);
466 }
467 
468 #undef __FUNCT__
469 #define __FUNCT__ "DMDARestoreAdicArray"
470 /*@C
471      DMDARestoreAdicArray - Restores an array of derivative types for a DMDA
472 
473     Input Parameter:
474 +    da - information about my local patch
475 -    ghosted - do you want arrays for the ghosted or nonghosted patch
476 
477     Output Parameters:
478 +    ptr - array data structured to be passed to ad_FormFunctionLocal()
479 .    array_start - actual start of 1d array of all values that adiC can access directly
480 -    tdof - total number of degrees of freedom represented in array_start
481 
482      Level: advanced
483 
484 .seealso: DMDAGetAdicArray()
485 
486 @*/
487 PetscErrorCode  DMDARestoreAdicArray(DM da,PetscBool  ghosted,void *ptr,void *array_start,PetscInt *tdof)
488 {
489   PetscInt  i;
490   void      **iptr = (void**)ptr,iarray_start = 0;
491 
492   PetscFunctionBegin;
493   PetscValidHeaderSpecific(da,DM_CLASSID,1);
494   if (ghosted) {
495     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
496       if (dd->adarrayghostedout[i] == *iptr) {
497         iarray_start             = dd->adstartghostedout[i];
498         dd->adarrayghostedout[i] = PETSC_NULL;
499         dd->adstartghostedout[i] = PETSC_NULL;
500         break;
501       }
502     }
503     if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
504     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
505       if (!dd->adarrayghostedin[i]){
506         dd->adarrayghostedin[i] = *iptr;
507         dd->adstartghostedin[i] = iarray_start;
508         break;
509       }
510     }
511   } else {
512     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
513       if (dd->adarrayout[i] == *iptr) {
514         iarray_start      = dd->adstartout[i];
515         dd->adarrayout[i] = PETSC_NULL;
516         dd->adstartout[i] = PETSC_NULL;
517         break;
518       }
519     }
520     if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
521     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
522       if (!dd->adarrayin[i]){
523         dd->adarrayin[i]   = *iptr;
524         dd->adstartin[i]   = iarray_start;
525         break;
526       }
527     }
528   }
529   PetscFunctionReturn(0);
530 }
531 
532 #undef __FUNCT__
533 #define __FUNCT__ "ad_DAGetArray"
534 PetscErrorCode  ad_DAGetArray(DM da,PetscBool  ghosted,void *iptr)
535 {
536   PetscErrorCode ierr;
537   PetscFunctionBegin;
538   ierr = DMDAGetAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr);
539   PetscFunctionReturn(0);
540 }
541 
542 #undef __FUNCT__
543 #define __FUNCT__ "ad_DARestoreArray"
544 PetscErrorCode  ad_DARestoreArray(DM da,PetscBool  ghosted,void *iptr)
545 {
546   PetscErrorCode ierr;
547   PetscFunctionBegin;
548   ierr = DMDARestoreAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr);
549   PetscFunctionReturn(0);
550 }
551 
552 #endif
553 
554 #undef __FUNCT__
555 #define __FUNCT__ "DMDAGetArray"
556 /*@C
557      DMDAGetArray - Gets a work array for a DMDA
558 
559     Input Parameter:
560 +    da - information about my local patch
561 -    ghosted - do you want arrays for the ghosted or nonghosted patch
562 
563     Output Parameters:
564 .    vptr - array data structured
565 
566     Note:  The vector values are NOT initialized and may have garbage in them, so you may need
567            to zero them.
568 
569   Level: advanced
570 
571 .seealso: DMDARestoreArray(), DMDAGetAdicArray()
572 
573 @*/
574 PetscErrorCode  DMDAGetArray(DM da,PetscBool  ghosted,void *vptr)
575 {
576   PetscErrorCode ierr;
577   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
578   char           *iarray_start;
579   void           **iptr = (void**)vptr;
580   DM_DA          *dd = (DM_DA*)da->data;
581 
582   PetscFunctionBegin;
583   PetscValidHeaderSpecific(da,DM_CLASSID,1);
584   if (ghosted) {
585     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
586       if (dd->arrayghostedin[i]) {
587         *iptr                 = dd->arrayghostedin[i];
588         iarray_start          = (char*)dd->startghostedin[i];
589         dd->arrayghostedin[i] = PETSC_NULL;
590         dd->startghostedin[i] = PETSC_NULL;
591 
592         goto done;
593       }
594     }
595     xs = dd->Xs;
596     ys = dd->Ys;
597     zs = dd->Zs;
598     xm = dd->Xe-dd->Xs;
599     ym = dd->Ye-dd->Ys;
600     zm = dd->Ze-dd->Zs;
601   } else {
602     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
603       if (dd->arrayin[i]) {
604         *iptr          = dd->arrayin[i];
605         iarray_start   = (char*)dd->startin[i];
606         dd->arrayin[i] = PETSC_NULL;
607         dd->startin[i] = PETSC_NULL;
608 
609         goto done;
610       }
611     }
612     xs = dd->xs;
613     ys = dd->ys;
614     zs = dd->zs;
615     xm = dd->xe-dd->xs;
616     ym = dd->ye-dd->ys;
617     zm = dd->ze-dd->zs;
618   }
619 
620   switch (dd->dim) {
621     case 1: {
622       void *ptr;
623 
624       ierr  = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
625 
626       ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
627       *iptr = (void*)ptr;
628       break;}
629     case 2: {
630       void **ptr;
631 
632       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
633 
634       ptr  = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
635       for(j=ys;j<ys+ym;j++) {
636         ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
637       }
638       *iptr = (void*)ptr;
639       break;}
640     case 3: {
641       void ***ptr,**bptr;
642 
643       ierr  = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
644 
645       ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
646       bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
647       for(i=zs;i<zs+zm;i++) {
648         ptr[i] = bptr + ((i-zs)*ym - ys);
649       }
650       for (i=zs; i<zs+zm; i++) {
651         for (j=ys; j<ys+ym; j++) {
652           ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
653         }
654       }
655 
656       *iptr = (void*)ptr;
657       break;}
658     default:
659       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
660   }
661 
662   done:
663   /* add arrays to the checked out list */
664   if (ghosted) {
665     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
666       if (!dd->arrayghostedout[i]) {
667         dd->arrayghostedout[i] = *iptr ;
668         dd->startghostedout[i] = iarray_start;
669         break;
670       }
671     }
672   } else {
673     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
674       if (!dd->arrayout[i]) {
675         dd->arrayout[i] = *iptr ;
676         dd->startout[i] = iarray_start;
677         break;
678       }
679     }
680   }
681   PetscFunctionReturn(0);
682 }
683 
684 #undef __FUNCT__
685 #define __FUNCT__ "DMDARestoreArray"
686 /*@C
687      DMDARestoreArray - Restores an array of derivative types for a DMDA
688 
689     Input Parameter:
690 +    da - information about my local patch
691 .    ghosted - do you want arrays for the ghosted or nonghosted patch
692 -    vptr - array data structured to be passed to ad_FormFunctionLocal()
693 
694      Level: advanced
695 
696 .seealso: DMDAGetArray(), DMDAGetAdicArray()
697 
698 @*/
699 PetscErrorCode  DMDARestoreArray(DM da,PetscBool  ghosted,void *vptr)
700 {
701   PetscInt  i;
702   void      **iptr = (void**)vptr,*iarray_start = 0;
703   DM_DA     *dd = (DM_DA*)da->data;
704 
705   PetscFunctionBegin;
706   PetscValidHeaderSpecific(da,DM_CLASSID,1);
707   if (ghosted) {
708     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
709       if (dd->arrayghostedout[i] == *iptr) {
710         iarray_start           = dd->startghostedout[i];
711         dd->arrayghostedout[i] = PETSC_NULL;
712         dd->startghostedout[i] = PETSC_NULL;
713         break;
714       }
715     }
716     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
717       if (!dd->arrayghostedin[i]){
718         dd->arrayghostedin[i] = *iptr;
719         dd->startghostedin[i] = iarray_start;
720         break;
721       }
722     }
723   } else {
724     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
725       if (dd->arrayout[i] == *iptr) {
726         iarray_start    = dd->startout[i];
727         dd->arrayout[i] = PETSC_NULL;
728         dd->startout[i] = PETSC_NULL;
729         break;
730       }
731     }
732     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
733       if (!dd->arrayin[i]){
734         dd->arrayin[i]  = *iptr;
735         dd->startin[i]  = iarray_start;
736         break;
737       }
738     }
739   }
740   PetscFunctionReturn(0);
741 }
742 
743 #undef __FUNCT__
744 #define __FUNCT__ "DMDAGetAdicMFArray"
745 /*@C
746      DMDAGetAdicMFArray - Gets an array of derivative types for a DMDA for matrix-free ADIC.
747 
748      Input Parameter:
749 +    da - information about my local patch
750 -    ghosted - do you want arrays for the ghosted or nonghosted patch?
751 
752      Output Parameters:
753 +    vptr - array data structured to be passed to ad_FormFunctionLocal()
754 .    array_start - actual start of 1d array of all values that adiC can access directly (may be null)
755 -    tdof - total number of degrees of freedom represented in array_start (may be null)
756 
757      Notes:
758      The vector values are NOT initialized and may have garbage in them, so you may need
759      to zero them.
760 
761      This routine returns the same type of object as the DMDAVecGetArray(), except its
762      elements are derivative types instead of PetscScalars.
763 
764      Level: advanced
765 
766 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray()
767 
768 @*/
769 PetscErrorCode  DMDAGetAdicMFArray(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
770 {
771   PetscErrorCode ierr;
772   PetscInt       j,i,xs,ys,xm,ym,zs,zm,itdof = 0;
773   char           *iarray_start;
774   void           **iptr = (void**)vptr;
775   DM_DA          *dd = (DM_DA*)da->data;
776 
777   PetscFunctionBegin;
778   PetscValidHeaderSpecific(da,DM_CLASSID,1);
779   if (ghosted) {
780     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
781       if (dd->admfarrayghostedin[i]) {
782         *iptr                     = dd->admfarrayghostedin[i];
783         iarray_start              = (char*)dd->admfstartghostedin[i];
784         itdof                     = dd->ghostedtdof;
785         dd->admfarrayghostedin[i] = PETSC_NULL;
786         dd->admfstartghostedin[i] = PETSC_NULL;
787 
788         goto done;
789       }
790     }
791     xs = dd->Xs;
792     ys = dd->Ys;
793     zs = dd->Zs;
794     xm = dd->Xe-dd->Xs;
795     ym = dd->Ye-dd->Ys;
796     zm = dd->Ze-dd->Zs;
797   } else {
798     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
799       if (dd->admfarrayin[i]) {
800         *iptr              = dd->admfarrayin[i];
801         iarray_start       = (char*)dd->admfstartin[i];
802         itdof              = dd->tdof;
803         dd->admfarrayin[i] = PETSC_NULL;
804         dd->admfstartin[i] = PETSC_NULL;
805 
806         goto done;
807       }
808     }
809     xs = dd->xs;
810     ys = dd->ys;
811     zs = dd->zs;
812     xm = dd->xe-dd->xs;
813     ym = dd->ye-dd->ys;
814     zm = dd->ze-dd->zs;
815   }
816 
817   switch (dd->dim) {
818     case 1: {
819       void *ptr;
820       itdof = xm;
821 
822       ierr  = PetscMalloc(xm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
823 
824       ptr   = (void*)(iarray_start - xs*2*sizeof(PetscScalar));
825       *iptr = (void*)ptr;
826       break;}
827     case 2: {
828       void **ptr;
829       itdof = xm*ym;
830 
831       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
832 
833       ptr  = (void**)(iarray_start + xm*ym*2*sizeof(PetscScalar) - ys*sizeof(void*));
834       for(j=ys;j<ys+ym;j++) {
835         ptr[j] = iarray_start + 2*sizeof(PetscScalar)*(xm*(j-ys) - xs);
836       }
837       *iptr = (void*)ptr;
838       break;}
839     case 3: {
840       void ***ptr,**bptr;
841       itdof = xm*ym*zm;
842 
843       ierr  = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
844 
845       ptr  = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*));
846       bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**));
847       for(i=zs;i<zs+zm;i++) {
848         ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*);
849       }
850       for (i=zs; i<zs+zm; i++) {
851         for (j=ys; j<ys+ym; j++) {
852           ptr[i][j] = iarray_start + 2*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
853         }
854       }
855 
856       *iptr = (void*)ptr;
857       break;}
858     default:
859       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
860   }
861 
862   done:
863   /* add arrays to the checked out list */
864   if (ghosted) {
865     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
866       if (!dd->admfarrayghostedout[i]) {
867         dd->admfarrayghostedout[i] = *iptr ;
868         dd->admfstartghostedout[i] = iarray_start;
869         dd->ghostedtdof            = itdof;
870         break;
871       }
872     }
873   } else {
874     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
875       if (!dd->admfarrayout[i]) {
876         dd->admfarrayout[i] = *iptr ;
877         dd->admfstartout[i] = iarray_start;
878         dd->tdof            = itdof;
879         break;
880       }
881     }
882   }
883   if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained");
884   if (tdof)        *tdof = itdof;
885   if (array_start) *(void**)array_start = iarray_start;
886   PetscFunctionReturn(0);
887 }
888 
889 #undef __FUNCT__
890 #define __FUNCT__ "DMDAGetAdicMFArray4"
891 PetscErrorCode  DMDAGetAdicMFArray4(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
892 {
893   PetscErrorCode ierr;
894   PetscInt       j,i,xs,ys,xm,ym,itdof = 0;
895   char           *iarray_start;
896   void           **iptr = (void**)vptr;
897   DM_DA          *dd = (DM_DA*)da->data;
898 
899   PetscFunctionBegin;
900   PetscValidHeaderSpecific(da,DM_CLASSID,1);
901   if (ghosted) {
902     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
903       if (dd->admfarrayghostedin[i]) {
904         *iptr                     = dd->admfarrayghostedin[i];
905         iarray_start              = (char*)dd->admfstartghostedin[i];
906         itdof                     = dd->ghostedtdof;
907         dd->admfarrayghostedin[i] = PETSC_NULL;
908         dd->admfstartghostedin[i] = PETSC_NULL;
909 
910         goto done;
911       }
912     }
913     xs = dd->Xs;
914     ys = dd->Ys;
915     xm = dd->Xe-dd->Xs;
916     ym = dd->Ye-dd->Ys;
917   } else {
918     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
919       if (dd->admfarrayin[i]) {
920         *iptr              = dd->admfarrayin[i];
921         iarray_start       = (char*)dd->admfstartin[i];
922         itdof              = dd->tdof;
923         dd->admfarrayin[i] = PETSC_NULL;
924         dd->admfstartin[i] = PETSC_NULL;
925 
926         goto done;
927       }
928     }
929     xs = dd->xs;
930     ys = dd->ys;
931     xm = dd->xe-dd->xs;
932     ym = dd->ye-dd->ys;
933   }
934 
935   switch (dd->dim) {
936     case 2: {
937       void **ptr;
938       itdof = xm*ym;
939 
940       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*5*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
941 
942       ptr  = (void**)(iarray_start + xm*ym*5*sizeof(PetscScalar) - ys*sizeof(void*));
943       for(j=ys;j<ys+ym;j++) {
944         ptr[j] = iarray_start + 5*sizeof(PetscScalar)*(xm*(j-ys) - xs);
945       }
946       *iptr = (void*)ptr;
947       break;}
948     default:
949       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
950   }
951 
952   done:
953   /* add arrays to the checked out list */
954   if (ghosted) {
955     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
956       if (!dd->admfarrayghostedout[i]) {
957         dd->admfarrayghostedout[i] = *iptr ;
958         dd->admfstartghostedout[i] = iarray_start;
959         dd->ghostedtdof            = itdof;
960         break;
961       }
962     }
963   } else {
964     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
965       if (!dd->admfarrayout[i]) {
966         dd->admfarrayout[i] = *iptr ;
967         dd->admfstartout[i] = iarray_start;
968         dd->tdof            = itdof;
969         break;
970       }
971     }
972   }
973   if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained");
974   if (tdof)        *tdof = itdof;
975   if (array_start) *(void**)array_start = iarray_start;
976   PetscFunctionReturn(0);
977 }
978 
979 #undef __FUNCT__
980 #define __FUNCT__ "DMDAGetAdicMFArray9"
981 PetscErrorCode  DMDAGetAdicMFArray9(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
982 {
983   PetscErrorCode ierr;
984   PetscInt       j,i,xs,ys,xm,ym,itdof = 0;
985   char           *iarray_start;
986   void           **iptr = (void**)vptr;
987   DM_DA          *dd = (DM_DA*)da->data;
988 
989   PetscFunctionBegin;
990   PetscValidHeaderSpecific(da,DM_CLASSID,1);
991   if (ghosted) {
992     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
993       if (dd->admfarrayghostedin[i]) {
994         *iptr                     = dd->admfarrayghostedin[i];
995         iarray_start              = (char*)dd->admfstartghostedin[i];
996         itdof                     = dd->ghostedtdof;
997         dd->admfarrayghostedin[i] = PETSC_NULL;
998         dd->admfstartghostedin[i] = PETSC_NULL;
999 
1000         goto done;
1001       }
1002     }
1003     xs = dd->Xs;
1004     ys = dd->Ys;
1005     xm = dd->Xe-dd->Xs;
1006     ym = dd->Ye-dd->Ys;
1007   } else {
1008     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1009       if (dd->admfarrayin[i]) {
1010         *iptr              = dd->admfarrayin[i];
1011         iarray_start       = (char*)dd->admfstartin[i];
1012         itdof              = dd->tdof;
1013         dd->admfarrayin[i] = PETSC_NULL;
1014         dd->admfstartin[i] = PETSC_NULL;
1015 
1016         goto done;
1017       }
1018     }
1019     xs = dd->xs;
1020     ys = dd->ys;
1021     xm = dd->xe-dd->xs;
1022     ym = dd->ye-dd->ys;
1023   }
1024 
1025   switch (dd->dim) {
1026     case 2: {
1027       void **ptr;
1028       itdof = xm*ym;
1029 
1030       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*10*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1031 
1032       ptr  = (void**)(iarray_start + xm*ym*10*sizeof(PetscScalar) - ys*sizeof(void*));
1033       for(j=ys;j<ys+ym;j++) {
1034         ptr[j] = iarray_start + 10*sizeof(PetscScalar)*(xm*(j-ys) - xs);
1035       }
1036       *iptr = (void*)ptr;
1037       break;}
1038     default:
1039       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1040   }
1041 
1042   done:
1043   /* add arrays to the checked out list */
1044   if (ghosted) {
1045     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1046       if (!dd->admfarrayghostedout[i]) {
1047         dd->admfarrayghostedout[i] = *iptr ;
1048         dd->admfstartghostedout[i] = iarray_start;
1049         dd->ghostedtdof            = itdof;
1050         break;
1051       }
1052     }
1053   } else {
1054     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1055       if (!dd->admfarrayout[i]) {
1056         dd->admfarrayout[i] = *iptr ;
1057         dd->admfstartout[i] = iarray_start;
1058         dd->tdof            = itdof;
1059         break;
1060       }
1061     }
1062   }
1063   if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained");
1064   if (tdof)        *tdof = itdof;
1065   if (array_start) *(void**)array_start = iarray_start;
1066   PetscFunctionReturn(0);
1067 }
1068 
1069 #undef __FUNCT__
1070 #define __FUNCT__ "DMDAGetAdicMFArrayb"
1071 /*@C
1072      DMDAGetAdicMFArrayb - Gets an array of derivative types for a DMDA for matrix-free ADIC.
1073 
1074      Input Parameter:
1075 +    da - information about my local patch
1076 -    ghosted - do you want arrays for the ghosted or nonghosted patch?
1077 
1078      Output Parameters:
1079 +    vptr - array data structured to be passed to ad_FormFunctionLocal()
1080 .    array_start - actual start of 1d array of all values that adiC can access directly (may be null)
1081 -    tdof - total number of degrees of freedom represented in array_start (may be null)
1082 
1083      Notes:
1084      The vector values are NOT initialized and may have garbage in them, so you may need
1085      to zero them.
1086 
1087      This routine returns the same type of object as the DMDAVecGetArray(), except its
1088      elements are derivative types instead of PetscScalars.
1089 
1090      Level: advanced
1091 
1092 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray()
1093 
1094 @*/
1095 PetscErrorCode  DMDAGetAdicMFArrayb(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
1096 {
1097   PetscErrorCode ierr;
1098   PetscInt       j,i,xs,ys,xm,ym,zs,zm,itdof = 0;
1099   char           *iarray_start;
1100   void           **iptr = (void**)vptr;
1101   DM_DA          *dd = (DM_DA*)da->data;
1102   PetscInt       bs = dd->w,bs1 = bs+1;
1103 
1104   PetscFunctionBegin;
1105   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1106   if (ghosted) {
1107     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1108       if (dd->admfarrayghostedin[i]) {
1109         *iptr                     = dd->admfarrayghostedin[i];
1110         iarray_start              = (char*)dd->admfstartghostedin[i];
1111         itdof                     = dd->ghostedtdof;
1112         dd->admfarrayghostedin[i] = PETSC_NULL;
1113         dd->admfstartghostedin[i] = PETSC_NULL;
1114 
1115         goto done;
1116       }
1117     }
1118     xs = dd->Xs;
1119     ys = dd->Ys;
1120     zs = dd->Zs;
1121     xm = dd->Xe-dd->Xs;
1122     ym = dd->Ye-dd->Ys;
1123     zm = dd->Ze-dd->Zs;
1124   } else {
1125     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1126       if (dd->admfarrayin[i]) {
1127         *iptr              = dd->admfarrayin[i];
1128         iarray_start       = (char*)dd->admfstartin[i];
1129         itdof              = dd->tdof;
1130         dd->admfarrayin[i] = PETSC_NULL;
1131         dd->admfstartin[i] = PETSC_NULL;
1132 
1133         goto done;
1134       }
1135     }
1136     xs = dd->xs;
1137     ys = dd->ys;
1138     zs = dd->zs;
1139     xm = dd->xe-dd->xs;
1140     ym = dd->ye-dd->ys;
1141     zm = dd->ze-dd->zs;
1142   }
1143 
1144   switch (dd->dim) {
1145     case 1: {
1146       void *ptr;
1147       itdof = xm;
1148 
1149       ierr  = PetscMalloc(xm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1150 
1151       ptr   = (void*)(iarray_start - xs*bs1*sizeof(PetscScalar));
1152       *iptr = (void*)ptr;
1153       break;}
1154     case 2: {
1155       void **ptr;
1156       itdof = xm*ym;
1157 
1158       ierr  = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1159 
1160       ptr  = (void**)(iarray_start + xm*ym*bs1*sizeof(PetscScalar) - ys*sizeof(void*));
1161       for(j=ys;j<ys+ym;j++) {
1162         ptr[j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*(j-ys) - xs);
1163       }
1164       *iptr = (void*)ptr;
1165       break;}
1166     case 3: {
1167       void ***ptr,**bptr;
1168       itdof = xm*ym*zm;
1169 
1170       ierr  = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1171 
1172       ptr  = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*));
1173       bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**));
1174       for(i=zs;i<zs+zm;i++) {
1175         ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*);
1176       }
1177       for (i=zs; i<zs+zm; i++) {
1178         for (j=ys; j<ys+ym; j++) {
1179           ptr[i][j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1180         }
1181       }
1182 
1183       *iptr = (void*)ptr;
1184       break;}
1185     default:
1186       SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1187   }
1188 
1189   done:
1190   /* add arrays to the checked out list */
1191   if (ghosted) {
1192     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1193       if (!dd->admfarrayghostedout[i]) {
1194         dd->admfarrayghostedout[i] = *iptr ;
1195         dd->admfstartghostedout[i] = iarray_start;
1196         dd->ghostedtdof            = itdof;
1197         break;
1198       }
1199     }
1200   } else {
1201     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1202       if (!dd->admfarrayout[i]) {
1203         dd->admfarrayout[i] = *iptr ;
1204         dd->admfstartout[i] = iarray_start;
1205         dd->tdof            = itdof;
1206         break;
1207       }
1208     }
1209   }
1210   if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained");
1211   if (tdof)        *tdof = itdof;
1212   if (array_start) *(void**)array_start = iarray_start;
1213   PetscFunctionReturn(0);
1214 }
1215 
1216 #undef __FUNCT__
1217 #define __FUNCT__ "DMDARestoreAdicMFArray"
1218 /*@C
1219      DMDARestoreAdicMFArray - Restores an array of derivative types for a DMDA.
1220 
1221      Input Parameter:
1222 +    da - information about my local patch
1223 -    ghosted - do you want arrays for the ghosted or nonghosted patch?
1224 
1225      Output Parameters:
1226 +    ptr - array data structure to be passed to ad_FormFunctionLocal()
1227 .    array_start - actual start of 1d array of all values that adiC can access directly
1228 -    tdof - total number of degrees of freedom represented in array_start
1229 
1230      Level: advanced
1231 
1232 .seealso: DMDAGetAdicArray()
1233 
1234 @*/
1235 PetscErrorCode  DMDARestoreAdicMFArray(DM da,PetscBool  ghosted,void *vptr,void *array_start,PetscInt *tdof)
1236 {
1237   PetscInt  i;
1238   void      **iptr = (void**)vptr,*iarray_start = 0;
1239   DM_DA     *dd = (DM_DA*)da->data;
1240 
1241   PetscFunctionBegin;
1242   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1243   if (ghosted) {
1244     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1245       if (dd->admfarrayghostedout[i] == *iptr) {
1246         iarray_start               = dd->admfstartghostedout[i];
1247         dd->admfarrayghostedout[i] = PETSC_NULL;
1248         dd->admfstartghostedout[i] = PETSC_NULL;
1249         break;
1250       }
1251     }
1252     if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
1253     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1254       if (!dd->admfarrayghostedin[i]){
1255         dd->admfarrayghostedin[i] = *iptr;
1256         dd->admfstartghostedin[i] = iarray_start;
1257         break;
1258       }
1259     }
1260   } else {
1261     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1262       if (dd->admfarrayout[i] == *iptr) {
1263         iarray_start        = dd->admfstartout[i];
1264         dd->admfarrayout[i] = PETSC_NULL;
1265         dd->admfstartout[i] = PETSC_NULL;
1266         break;
1267       }
1268     }
1269     if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
1270     for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) {
1271       if (!dd->admfarrayin[i]){
1272         dd->admfarrayin[i] = *iptr;
1273         dd->admfstartin[i] = iarray_start;
1274         break;
1275       }
1276     }
1277   }
1278   PetscFunctionReturn(0);
1279 }
1280 
1281 #undef __FUNCT__
1282 #define __FUNCT__ "admf_DAGetArray"
1283 PetscErrorCode  admf_DAGetArray(DM da,PetscBool  ghosted,void *iptr)
1284 {
1285   PetscErrorCode ierr;
1286   PetscFunctionBegin;
1287   ierr = DMDAGetAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr);
1288   PetscFunctionReturn(0);
1289 }
1290 
1291 #undef __FUNCT__
1292 #define __FUNCT__ "admf_DARestoreArray"
1293 PetscErrorCode  admf_DARestoreArray(DM da,PetscBool  ghosted,void *iptr)
1294 {
1295   PetscErrorCode ierr;
1296   PetscFunctionBegin;
1297   ierr = DMDARestoreAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr);
1298   PetscFunctionReturn(0);
1299 }
1300 
1301