Actual source code: dalocal.c

      1: #define PETSCDM_DLL

      3: /*
      4:   Code for manipulating distributed regular arrays in parallel.
      5: */

     7:  #include src/dm/da/daimpl.h

      9: /*
     10:    This allows the DA vectors to properly tell Matlab their dimensions
     11: */
     12: #if defined(PETSC_HAVE_MATLAB)
     13: #include "engine.h"   /* Matlab include file */
     14: #include "mex.h"      /* Matlab include file */
     18: PetscErrorCode PETSCDM_DLLEXPORT VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
     19: {
     21:   PetscInt n,m;
     22:   Vec          vec = (Vec)obj;
     23:   PetscScalar  *array;
     24:   mxArray      *mat;
     25:   DA           da;

     28:   PetscObjectQuery((PetscObject)vec,"DA",(PetscObject*)&da);
     29:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DA");
     30:   DAGetGhostCorners(da,0,0,0,&m,&n,0);

     32:   VecGetArray(vec,&array);
     33: #if !defined(PETSC_USE_COMPLEX)
     34:   mat  = mxCreateDoubleMatrix(m,n,mxREAL);
     35: #else
     36:   mat  = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
     37: #endif
     38:   PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));
     39:   PetscObjectName(obj);
     40:   engPutVariable((Engine *)mengine,obj->name,mat);
     41:
     42:   VecRestoreArray(vec,&array);
     43:   return(0);
     44: }
     46: #endif

     51: /*@C
     52:    DACreateLocalVector - Creates a Seq PETSc vector that
     53:    may be used with the DAXXX routines.

     55:    Not Collective

     57:    Input Parameter:
     58: .  da - the distributed array

     60:    Output Parameter:
     61: .  g - the local vector

     63:    Level: beginner

     65:    Note:
     66:    The output parameter, g, is a regular PETSc vector that should be destroyed
     67:    with a call to VecDestroy() when usage is finished.

     69: .keywords: distributed array, create, local, vector

     71: .seealso: DACreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
     72:           DACreate1d(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(),
     73:           DAGlobalToLocalEnd(), DALocalToGlobal(), DAGetLocalVector(), DARestoreLocalVector()
     74: @*/
     75: PetscErrorCode PETSCDM_DLLEXPORT DACreateLocalVector(DA da,Vec* g)
     76: {

     82:   VecCreateSeq(PETSC_COMM_SELF,da->nlocal,g);
     83:   VecSetBlockSize(*g,da->w);
     84:   PetscObjectCompose((PetscObject)*g,"DA",(PetscObject)da);
     85: #if defined(PETSC_HAVE_MATLAB)
     86:   if (da->w == 1  && da->dim == 2) {
     87:     PetscObjectComposeFunctionDynamic((PetscObject)*g,"PetscMatlabEnginePut_C","VecMatlabEnginePut_DA2d",VecMatlabEnginePut_DA2d);
     88:   }
     89: #endif
     90:   return(0);
     91: }

     95: /*@C
     96:    DAGetLocalVector - Gets a Seq PETSc vector that
     97:    may be used with the DAXXX routines.

     99:    Not Collective

    101:    Input Parameter:
    102: .  da - the distributed array

    104:    Output Parameter:
    105: .  g - the local vector

    107:    Level: beginner

    109:    Note:
    110:    The vector values are NOT initialized and may have garbage in them, so you may need
    111:    to zero them.

    113:    The output parameter, g, is a regular PETSc vector that should be returned with
    114:    DARestoreLocalVector() DO NOT call VecDestroy() on it.

    116:    VecStride*() operations can be useful when using DA with dof > 1

    118: .keywords: distributed array, create, local, vector

    120: .seealso: DACreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
    121:           DACreate1d(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(),
    122:           DAGlobalToLocalEnd(), DALocalToGlobal(), DACreateLocalVector(), DARestoreLocalVector(),
    123:           VecStrideMax(), VecStrideMin(), VecStrideNorm()
    124: @*/
    125: PetscErrorCode PETSCDM_DLLEXPORT DAGetLocalVector(DA da,Vec* g)
    126: {
    127:   PetscErrorCode ierr,i;

    132:   for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    133:     if (da->localin[i]) {
    134:       *g             = da->localin[i];
    135:       da->localin[i] = PETSC_NULL;
    136:       goto alldone;
    137:     }
    138:   }
    139:   DACreateLocalVector(da,g);

    141:   alldone:
    142:   for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    143:     if (!da->localout[i]) {
    144:       da->localout[i] = *g;
    145:       break;
    146:     }
    147:   }
    148:   return(0);
    149: }

    153: /*@C
    154:    DARestoreLocalVector - Returns a Seq PETSc vector that
    155:      obtained from DAGetLocalVector(). Do not use with vector obtained via
    156:      DACreateLocalVector().

    158:    Not Collective

    160:    Input Parameter:
    161: +  da - the distributed array
    162: -  g - the local vector

    164:    Level: beginner

    166: .keywords: distributed array, create, local, vector

    168: .seealso: DACreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
    169:           DACreate1d(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(),
    170:           DAGlobalToLocalEnd(), DALocalToGlobal(), DACreateLocalVector(), DAGetLocalVector()
    171: @*/
    172: PetscErrorCode PETSCDM_DLLEXPORT DARestoreLocalVector(DA da,Vec* g)
    173: {
    175:   PetscInt i,j;

    180:   for (j=0; j<DA_MAX_WORK_VECTORS; j++) {
    181:     if (*g == da->localout[j]) {
    182:       da->localout[j] = PETSC_NULL;
    183:       for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    184:         if (!da->localin[i]) {
    185:           da->localin[i] = *g;
    186:           goto alldone;
    187:         }
    188:       }
    189:     }
    190:   }
    191:   VecDestroy(*g);
    192:   alldone:
    193:   return(0);
    194: }

    198: /*@C
    199:    DAGetGlobalVector - Gets a MPI PETSc vector that
    200:    may be used with the DAXXX routines.

    202:    Collective on DA

    204:    Input Parameter:
    205: .  da - the distributed array

    207:    Output Parameter:
    208: .  g - the global vector

    210:    Level: beginner

    212:    Note:
    213:    The vector values are NOT initialized and may have garbage in them, so you may need
    214:    to zero them.

    216:    The output parameter, g, is a regular PETSc vector that should be returned with
    217:    DARestoreGlobalVector() DO NOT call VecDestroy() on it.

    219:    VecStride*() operations can be useful when using DA with dof > 1

    221: .keywords: distributed array, create, Global, vector

    223: .seealso: DACreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
    224:           DACreate1d(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(),
    225:           DAGlobalToLocalEnd(), DALocalToGlobal(), DACreateLocalVector(), DARestoreLocalVector()
    226:           VecStrideMax(), VecStrideMin(), VecStrideNorm()

    228: @*/
    229: PetscErrorCode PETSCDM_DLLEXPORT DAGetGlobalVector(DA da,Vec* g)
    230: {
    232:   PetscInt i;

    237:   for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    238:     if (da->globalin[i]) {
    239:       *g             = da->globalin[i];
    240:       da->globalin[i] = PETSC_NULL;
    241:       goto alldone;
    242:     }
    243:   }
    244:   DACreateGlobalVector(da,g);

    246:   alldone:
    247:   for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    248:     if (!da->globalout[i]) {
    249:       da->globalout[i] = *g;
    250:       break;
    251:     }
    252:   }
    253:   return(0);
    254: }

    258: /*@C
    259:    DARestoreGlobalVector - Returns a Seq PETSc vector that
    260:      obtained from DAGetGlobalVector(). Do not use with vector obtained via
    261:      DACreateGlobalVector().

    263:    Not Collective

    265:    Input Parameter:
    266: +  da - the distributed array
    267: -  g - the global vector

    269:    Level: beginner

    271: .keywords: distributed array, create, global, vector

    273: .seealso: DACreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
    274:           DACreate1d(), DACreate2d(), DACreate3d(), DAGlobalToGlobalBegin(),
    275:           DAGlobalToGlobalEnd(), DAGlobalToGlobal(), DACreateLocalVector(), DAGetGlobalVector()
    276: @*/
    277: PetscErrorCode PETSCDM_DLLEXPORT DARestoreGlobalVector(DA da,Vec* g)
    278: {
    280:   PetscInt i,j;

    285:   for (j=0; j<DA_MAX_WORK_VECTORS; j++) {
    286:     if (*g == da->globalout[j]) {
    287:       da->globalout[j] = PETSC_NULL;
    288:       for (i=0; i<DA_MAX_WORK_VECTORS; i++) {
    289:         if (!da->globalin[i]) {
    290:           da->globalin[i] = *g;
    291:           goto alldone;
    292:         }
    293:       }
    294:     }
    295:   }
    296:   VecDestroy(*g);
    297:   alldone:
    298:   return(0);
    299: }

    301: /* ------------------------------------------------------------------- */
    302: #if defined(PETSC_HAVE_ADIC)

    305: #include "adic/ad_utils.h"

    310: /*@C
    311:      DAGetAdicArray - Gets an array of derivative types for a DA
    312:
    313:     Input Parameter:
    314: +    da - information about my local patch
    315: -    ghosted - do you want arrays for the ghosted or nonghosted patch

    317:     Output Parameters:
    318: +    ptr - array data structured to be passed to ad_FormFunctionLocal()
    319: .    array_start - actual start of 1d array of all values that adiC can access directly (may be null)
    320: -    tdof - total number of degrees of freedom represented in array_start (may be null)

    322:      Notes:
    323:        The vector values are NOT initialized and may have garbage in them, so you may need
    324:        to zero them.

    326:        Returns the same type of object as the DAVecGetArray() except its elements are
    327:            derivative types instead of PetscScalars

    329:      Level: advanced

    331: .seealso: DARestoreAdicArray()

    333: @*/
    334: PetscErrorCode PETSCDM_DLLEXPORT DAGetAdicArray(DA da,PetscTruth ghosted,void **iptr,void **array_start,PetscInt *tdof)
    335: {
    337:   PetscInt j,i,deriv_type_size,xs,ys,xm,ym,zs,zm,itdof;
    338:   char *iarray_start;

    342:   if (ghosted) {
    343:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    344:       if (da->adarrayghostedin[i]) {
    345:         *iptr                   = da->adarrayghostedin[i];
    346:         iarray_start            = (char*)da->adstartghostedin[i];
    347:         itdof                   = da->ghostedtdof;
    348:         da->adarrayghostedin[i] = PETSC_NULL;
    349:         da->adstartghostedin[i] = PETSC_NULL;
    350:
    351:         goto done;
    352:       }
    353:     }
    354:     xs = da->Xs;
    355:     ys = da->Ys;
    356:     zs = da->Zs;
    357:     xm = da->Xe-da->Xs;
    358:     ym = da->Ye-da->Ys;
    359:     zm = da->Ze-da->Zs;
    360:   } else {
    361:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    362:       if (da->adarrayin[i]) {
    363:         *iptr            = da->adarrayin[i];
    364:         iarray_start     = (char*)da->adstartin[i];
    365:         itdof            = da->tdof;
    366:         da->adarrayin[i] = PETSC_NULL;
    367:         da->adstartin[i] = PETSC_NULL;
    368:
    369:         goto done;
    370:       }
    371:     }
    372:     xs = da->xs;
    373:     ys = da->ys;
    374:     zs = da->zs;
    375:     xm = da->xe-da->xs;
    376:     ym = da->ye-da->ys;
    377:     zm = da->ze-da->zs;
    378:   }
    379:   deriv_type_size = PetscADGetDerivTypeSize();

    381:   switch (da->dim) {
    382:     case 1: {
    383:       void *ptr;
    384:       itdof = xm;

    386:       PetscMalloc(xm*deriv_type_size,&iarray_start);

    388:       ptr   = (void*)(iarray_start - xs*deriv_type_size);
    389:       *iptr = (void*)ptr;
    390:       break;}
    391:     case 2: {
    392:       void **ptr;
    393:       itdof = xm*ym;

    395:       PetscMalloc((ym+1)*sizeof(void*)+xm*ym*deriv_type_size,&iarray_start);

    397:       ptr  = (void**)(iarray_start + xm*ym*deriv_type_size - ys*sizeof(void*));
    398:       for(j=ys;j<ys+ym;j++) {
    399:         ptr[j] = iarray_start + deriv_type_size*(xm*(j-ys) - xs);
    400:       }
    401:       *iptr = (void*)ptr;
    402:       break;}
    403:     case 3: {
    404:       void ***ptr,**bptr;
    405:       itdof = xm*ym*zm;

    407:       PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*deriv_type_size,&iarray_start);

    409:       ptr  = (void***)(iarray_start + xm*ym*zm*deriv_type_size - zs*sizeof(void*));
    410:       bptr = (void**)(iarray_start + xm*ym*zm*deriv_type_size + zm*sizeof(void**));

    412:       for(i=zs;i<zs+zm;i++) {
    413:         ptr[i] = bptr + ((i-zs)*ym - ys);
    414:       }
    415:       for (i=zs; i<zs+zm; i++) {
    416:         for (j=ys; j<ys+ym; j++) {
    417:           ptr[i][j] = iarray_start + deriv_type_size*(xm*ym*(i-zs) + xm*(j-ys) - xs);
    418:         }
    419:       }

    421:       *iptr = (void*)ptr;
    422:       break;}
    423:     default:
    424:       SETERRQ1(PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
    425:   }

    427:   done:
    428:   /* add arrays to the checked out list */
    429:   if (ghosted) {
    430:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    431:       if (!da->adarrayghostedout[i]) {
    432:         da->adarrayghostedout[i] = *iptr ;
    433:         da->adstartghostedout[i] = iarray_start;
    434:         da->ghostedtdof          = itdof;
    435:         break;
    436:       }
    437:     }
    438:   } else {
    439:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    440:       if (!da->adarrayout[i]) {
    441:         da->adarrayout[i] = *iptr ;
    442:         da->adstartout[i] = iarray_start;
    443:         da->tdof          = itdof;
    444:         break;
    445:       }
    446:     }
    447:   }
    448:   if (i == DA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_ERR_SUP,"Too many DA ADIC arrays obtained");
    449:   if (tdof)        *tdof = itdof;
    450:   if (array_start) *array_start = iarray_start;
    451:   return(0);
    452: }

    456: /*@C
    457:      DARestoreAdicArray - Restores an array of derivative types for a DA
    458:
    459:     Input Parameter:
    460: +    da - information about my local patch
    461: -    ghosted - do you want arrays for the ghosted or nonghosted patch

    463:     Output Parameters:
    464: +    ptr - array data structured to be passed to ad_FormFunctionLocal()
    465: .    array_start - actual start of 1d array of all values that adiC can access directly
    466: -    tdof - total number of degrees of freedom represented in array_start

    468:      Level: advanced

    470: .seealso: DAGetAdicArray()

    472: @*/
    473: PetscErrorCode PETSCDM_DLLEXPORT DARestoreAdicArray(DA da,PetscTruth ghosted,void **iptr,void **array_start,PetscInt *tdof)
    474: {
    475:   PetscInt  i;
    476:   void *iarray_start = 0;
    477:
    480:   if (ghosted) {
    481:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    482:       if (da->adarrayghostedout[i] == *iptr) {
    483:         iarray_start             = da->adstartghostedout[i];
    484:         da->adarrayghostedout[i] = PETSC_NULL;
    485:         da->adstartghostedout[i] = PETSC_NULL;
    486:         break;
    487:       }
    488:     }
    489:     if (!iarray_start) SETERRQ(PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
    490:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    491:       if (!da->adarrayghostedin[i]){
    492:         da->adarrayghostedin[i] = *iptr;
    493:         da->adstartghostedin[i] = iarray_start;
    494:         break;
    495:       }
    496:     }
    497:   } else {
    498:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    499:       if (da->adarrayout[i] == *iptr) {
    500:         iarray_start      = da->adstartout[i];
    501:         da->adarrayout[i] = PETSC_NULL;
    502:         da->adstartout[i] = PETSC_NULL;
    503:         break;
    504:       }
    505:     }
    506:     if (!iarray_start) SETERRQ(PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
    507:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    508:       if (!da->adarrayin[i]){
    509:         da->adarrayin[i]   = *iptr;
    510:         da->adstartin[i]   = iarray_start;
    511:         break;
    512:       }
    513:     }
    514:   }
    515:   return(0);
    516: }

    520: PetscErrorCode PETSCDM_DLLEXPORT ad_DAGetArray(DA da,PetscTruth ghosted,void **iptr)
    521: {
    524:   DAGetAdicArray(da,ghosted,iptr,0,0);
    525:   return(0);
    526: }

    530: PetscErrorCode PETSCDM_DLLEXPORT ad_DARestoreArray(DA da,PetscTruth ghosted,void **iptr)
    531: {
    534:   DARestoreAdicArray(da,ghosted,iptr,0,0);
    535:   return(0);
    536: }

    538: #endif

    542: /*@C
    543:      DAGetArray - Gets a work array for a DA
    544:
    545:     Input Parameter:
    546: +    da - information about my local patch
    547: -    ghosted - do you want arrays for the ghosted or nonghosted patch

    549:     Output Parameters:
    550: .    ptr - array data structured

    552:     Note:  The vector values are NOT initialized and may have garbage in them, so you may need
    553:            to zero them.

    555:   Level: advanced

    557: .seealso: DARestoreArray(), DAGetAdicArray()

    559: @*/
    560: PetscErrorCode PETSCDM_DLLEXPORT DAGetArray(DA da,PetscTruth ghosted,void **iptr)
    561: {
    563:   PetscInt j,i,xs,ys,xm,ym,zs,zm;
    564:   char *iarray_start;

    568:   if (ghosted) {
    569:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    570:       if (da->arrayghostedin[i]) {
    571:         *iptr                 = da->arrayghostedin[i];
    572:         iarray_start          = (char*)da->startghostedin[i];
    573:         da->arrayghostedin[i] = PETSC_NULL;
    574:         da->startghostedin[i] = PETSC_NULL;
    575:
    576:         goto done;
    577:       }
    578:     }
    579:     xs = da->Xs;
    580:     ys = da->Ys;
    581:     zs = da->Zs;
    582:     xm = da->Xe-da->Xs;
    583:     ym = da->Ye-da->Ys;
    584:     zm = da->Ze-da->Zs;
    585:   } else {
    586:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    587:       if (da->arrayin[i]) {
    588:         *iptr          = da->arrayin[i];
    589:         iarray_start   = (char*)da->startin[i];
    590:         da->arrayin[i] = PETSC_NULL;
    591:         da->startin[i] = PETSC_NULL;
    592:
    593:         goto done;
    594:       }
    595:     }
    596:     xs = da->xs;
    597:     ys = da->ys;
    598:     zs = da->zs;
    599:     xm = da->xe-da->xs;
    600:     ym = da->ye-da->ys;
    601:     zm = da->ze-da->zs;
    602:   }

    604:   switch (da->dim) {
    605:     case 1: {
    606:       void *ptr;

    608:       PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);

    610:       ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
    611:       *iptr = (void*)ptr;
    612:       break;}
    613:     case 2: {
    614:       void **ptr;

    616:       PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);

    618:       ptr  = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
    619:       for(j=ys;j<ys+ym;j++) {
    620:         ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
    621:       }
    622:       *iptr = (void*)ptr;
    623:       break;}
    624:     case 3: {
    625:       void ***ptr,**bptr;

    627:       PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);

    629:       ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
    630:       bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
    631:       for(i=zs;i<zs+zm;i++) {
    632:         ptr[i] = bptr + ((i-zs)*ym - ys);
    633:       }
    634:       for (i=zs; i<zs+zm; i++) {
    635:         for (j=ys; j<ys+ym; j++) {
    636:           ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
    637:         }
    638:       }

    640:       *iptr = (void*)ptr;
    641:       break;}
    642:     default:
    643:       SETERRQ1(PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
    644:   }

    646:   done:
    647:   /* add arrays to the checked out list */
    648:   if (ghosted) {
    649:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    650:       if (!da->arrayghostedout[i]) {
    651:         da->arrayghostedout[i] = *iptr ;
    652:         da->startghostedout[i] = iarray_start;
    653:         break;
    654:       }
    655:     }
    656:   } else {
    657:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    658:       if (!da->arrayout[i]) {
    659:         da->arrayout[i] = *iptr ;
    660:         da->startout[i] = iarray_start;
    661:         break;
    662:       }
    663:     }
    664:   }
    665:   return(0);
    666: }

    670: /*@C
    671:      DARestoreArray - Restores an array of derivative types for a DA
    672:
    673:     Input Parameter:
    674: +    da - information about my local patch
    675: .    ghosted - do you want arrays for the ghosted or nonghosted patch
    676: -    ptr - array data structured to be passed to ad_FormFunctionLocal()

    678:      Level: advanced

    680: .seealso: DAGetArray(), DAGetAdicArray()

    682: @*/
    683: PetscErrorCode PETSCDM_DLLEXPORT DARestoreArray(DA da,PetscTruth ghosted,void **iptr)
    684: {
    685:   PetscInt  i;
    686:   void *iarray_start = 0;
    687:
    690:   if (ghosted) {
    691:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    692:       if (da->arrayghostedout[i] == *iptr) {
    693:         iarray_start           = da->startghostedout[i];
    694:         da->arrayghostedout[i] = PETSC_NULL;
    695:         da->startghostedout[i] = PETSC_NULL;
    696:         break;
    697:       }
    698:     }
    699:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    700:       if (!da->arrayghostedin[i]){
    701:         da->arrayghostedin[i] = *iptr;
    702:         da->startghostedin[i] = iarray_start;
    703:         break;
    704:       }
    705:     }
    706:   } else {
    707:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    708:       if (da->arrayout[i] == *iptr) {
    709:         iarray_start    = da->startout[i];
    710:         da->arrayout[i] = PETSC_NULL;
    711:         da->startout[i] = PETSC_NULL;
    712:         break;
    713:       }
    714:     }
    715:     for (i=0; i<DA_MAX_WORK_ARRAYS; i++) {
    716:       if (!da->arrayin[i]){
    717:         da->arrayin[i]  = *iptr;
    718:         da->startin[i]  = iarray_start;
    719:         break;
    720:       }
    721:     }
    722:   }
    723:   return(0);
    724: }

    728: /*@C
    729:      DAGetAdicMFArray - Gets an array of derivative types for a DA for matrix-free ADIC.
    730:
    731:      Input Parameter:
    732: +    da - information about my local patch
    733: -    ghosted - do you want arrays for the ghosted or nonghosted patch?

    735:      Output Parameters:
    736: +    iptr - array data structured to be passed to ad_FormFunctionLocal()
    737: .    array_start - actual start of 1d array of all values that adiC can access directly (may be null)
    738: -    tdof - total number of degrees of freedom represented in array_start (may be null)

    740:      Notes:
    741:      The vector values are NOT initialized and may have garbage in them, so you may need
    742:      to zero them.

    744:      This routine returns the same type of object as the DAVecGetArray(), except its
    745:      elements are derivative types instead of PetscScalars.

    747:      Level: advanced

    749: .seealso: DARestoreAdicMFArray(), DAGetArray(), DAGetAdicArray()

    751: @*/
    752: PetscErrorCode PETSCDM_DLLEXPORT DAGetAdicMFArray(DA da,PetscTruth ghosted,void **iptr,void **array_start,PetscInt *tdof)
    753: {
    755:   PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0;
    756:   char *iarray_start;

    760:   if (ghosted) {
    761:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    762:       if (da->admfarrayghostedin[i]) {
    763:         *iptr                     = da->admfarrayghostedin[i];
    764:         iarray_start              = (char*)da->admfstartghostedin[i];
    765:         itdof                     = da->ghostedtdof;
    766:         da->admfarrayghostedin[i] = PETSC_NULL;
    767:         da->admfstartghostedin[i] = PETSC_NULL;
    768:
    769:         goto done;
    770:       }
    771:     }
    772:     xs = da->Xs;
    773:     ys = da->Ys;
    774:     zs = da->Zs;
    775:     xm = da->Xe-da->Xs;
    776:     ym = da->Ye-da->Ys;
    777:     zm = da->Ze-da->Zs;
    778:   } else {
    779:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    780:       if (da->admfarrayin[i]) {
    781:         *iptr              = da->admfarrayin[i];
    782:         iarray_start       = (char*)da->admfstartin[i];
    783:         itdof              = da->tdof;
    784:         da->admfarrayin[i] = PETSC_NULL;
    785:         da->admfstartin[i] = PETSC_NULL;
    786:
    787:         goto done;
    788:       }
    789:     }
    790:     xs = da->xs;
    791:     ys = da->ys;
    792:     zs = da->zs;
    793:     xm = da->xe-da->xs;
    794:     ym = da->ye-da->ys;
    795:     zm = da->ze-da->zs;
    796:   }

    798:   switch (da->dim) {
    799:     case 1: {
    800:       void *ptr;
    801:       itdof = xm;

    803:       PetscMalloc(xm*2*sizeof(PetscScalar),&iarray_start);

    805:       ptr   = (void*)(iarray_start - xs*2*sizeof(PetscScalar));
    806:       *iptr = (void*)ptr;
    807:       break;}
    808:     case 2: {
    809:       void **ptr;
    810:       itdof = xm*ym;

    812:       PetscMalloc((ym+1)*sizeof(void*)+xm*ym*2*sizeof(PetscScalar),&iarray_start);

    814:       ptr  = (void**)(iarray_start + xm*ym*2*sizeof(PetscScalar) - ys*sizeof(void*));
    815:       for(j=ys;j<ys+ym;j++) {
    816:         ptr[j] = iarray_start + 2*sizeof(PetscScalar)*(xm*(j-ys) - xs);
    817:       }
    818:       *iptr = (void*)ptr;
    819:       break;}
    820:     case 3: {
    821:       void ***ptr,**bptr;
    822:       itdof = xm*ym*zm;

    824:       PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*2*sizeof(PetscScalar),&iarray_start);

    826:       ptr  = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*));
    827:       bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**));
    828:       for(i=zs;i<zs+zm;i++) {
    829:         ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*);
    830:       }
    831:       for (i=zs; i<zs+zm; i++) {
    832:         for (j=ys; j<ys+ym; j++) {
    833:           ptr[i][j] = iarray_start + 2*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
    834:         }
    835:       }

    837:       *iptr = (void*)ptr;
    838:       break;}
    839:     default:
    840:       SETERRQ1(PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
    841:   }

    843:   done:
    844:   /* add arrays to the checked out list */
    845:   if (ghosted) {
    846:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    847:       if (!da->admfarrayghostedout[i]) {
    848:         da->admfarrayghostedout[i] = *iptr ;
    849:         da->admfstartghostedout[i] = iarray_start;
    850:         da->ghostedtdof            = itdof;
    851:         break;
    852:       }
    853:     }
    854:   } else {
    855:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    856:       if (!da->admfarrayout[i]) {
    857:         da->admfarrayout[i] = *iptr ;
    858:         da->admfstartout[i] = iarray_start;
    859:         da->tdof            = itdof;
    860:         break;
    861:       }
    862:     }
    863:   }
    864:   if (i == DA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_ERR_ARG_WRONG,"Too many DA ADIC arrays obtained");
    865:   if (tdof)        *tdof = itdof;
    866:   if (array_start) *array_start = iarray_start;
    867:   return(0);
    868: }

    872: /*@C
    873:      DARestoreAdicMFArray - Restores an array of derivative types for a DA.
    874:
    875:      Input Parameter:
    876: +    da - information about my local patch
    877: -    ghosted - do you want arrays for the ghosted or nonghosted patch?

    879:      Output Parameters:
    880: +    ptr - array data structure to be passed to ad_FormFunctionLocal()
    881: .    array_start - actual start of 1d array of all values that adiC can access directly
    882: -    tdof - total number of degrees of freedom represented in array_start

    884:      Level: advanced

    886: .seealso: DAGetAdicArray()

    888: @*/
    889: PetscErrorCode PETSCDM_DLLEXPORT DARestoreAdicMFArray(DA da,PetscTruth ghosted,void **iptr,void **array_start,PetscInt *tdof)
    890: {
    891:   PetscInt  i;
    892:   void *iarray_start = 0;
    893:
    896:   if (ghosted) {
    897:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    898:       if (da->admfarrayghostedout[i] == *iptr) {
    899:         iarray_start               = da->admfstartghostedout[i];
    900:         da->admfarrayghostedout[i] = PETSC_NULL;
    901:         da->admfstartghostedout[i] = PETSC_NULL;
    902:         break;
    903:       }
    904:     }
    905:     if (!iarray_start) SETERRQ(PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
    906:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    907:       if (!da->admfarrayghostedin[i]){
    908:         da->admfarrayghostedin[i] = *iptr;
    909:         da->admfstartghostedin[i] = iarray_start;
    910:         break;
    911:       }
    912:     }
    913:   } else {
    914:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    915:       if (da->admfarrayout[i] == *iptr) {
    916:         iarray_start        = da->admfstartout[i];
    917:         da->admfarrayout[i] = PETSC_NULL;
    918:         da->admfstartout[i] = PETSC_NULL;
    919:         break;
    920:       }
    921:     }
    922:     if (!iarray_start) SETERRQ(PETSC_ERR_ARG_WRONG,"Could not find array in checkout list");
    923:     for (i=0; i<DA_MAX_AD_ARRAYS; i++) {
    924:       if (!da->admfarrayin[i]){
    925:         da->admfarrayin[i] = *iptr;
    926:         da->admfstartin[i] = iarray_start;
    927:         break;
    928:       }
    929:     }
    930:   }
    931:   return(0);
    932: }

    936: PetscErrorCode PETSCDM_DLLEXPORT admf_DAGetArray(DA da,PetscTruth ghosted,void **iptr)
    937: {
    940:   DAGetAdicMFArray(da,ghosted,iptr,0,0);
    941:   return(0);
    942: }

    946: PetscErrorCode PETSCDM_DLLEXPORT admf_DARestoreArray(DA da,PetscTruth ghosted,void **iptr)
    947: {
    950:   DARestoreAdicMFArray(da,ghosted,iptr,0,0);
    951:   return(0);
    952: }