Actual source code: pack.c

      1: #define PETSCDM_DLL
      2:
     3:  #include petscda.h
     4:  #include petscmat.h

      6: /*
      7:    rstart is where an array/subvector starts in the global parallel vector, so arrays
      8:    rstarts are meaningless (and set to the previous one) except on processor 0
      9: */

     11: typedef enum {VECPACK_ARRAY, VECPACK_DA, VECPACK_VECSCATTER} VecPackLinkType;

     13: struct VecPackLink {
     14:   DA                 da;
     15:   PetscInt           n,rstart;      /* rstart is relative to this processor */
     16:   VecPackLinkType    type;
     17:   struct VecPackLink *next;
     18: };

     20: typedef struct _VecPackOps *VecPackOps;
     21: struct _VecPackOps {
     22:   PetscErrorCode (*view)(VecPack,PetscViewer);
     23:   PetscErrorCode (*createglobalvector)(VecPack,Vec*);
     24:   PetscErrorCode (*getcoloring)(VecPack,ISColoringType,ISColoring*);
     25:   PetscErrorCode (*getmatrix)(VecPack,MatType,Mat*);
     26:   PetscErrorCode (*getinterpolation)(VecPack,VecPack,Mat*,Vec*);
     27:   PetscErrorCode (*refine)(VecPack,MPI_Comm,VecPack*);
     28: };

     30: struct _p_VecPack {
     31:   PETSCHEADER(struct _VecPackOps);
     32:   PetscMPIInt        rank;
     33:   PetscInt           n,N,rstart;   /* rstart is relative to all processors */
     34:   Vec                globalvector;
     35:   PetscInt           nDA,nredundant;
     36:   struct VecPackLink *next;
     37: };

     41: /*@C
     42:     VecPackCreate - Creates a vector packer, used to generate "composite"
     43:       vectors made up of several subvectors.

     45:     Collective on MPI_Comm

     47:     Input Parameter:
     48: .   comm - the processors that will share the global vector

     50:     Output Parameters:
     51: .   packer - the packer object

     53:     Level: advanced

     55: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
     56:          VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()
     57:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

     59: @*/
     60: PetscErrorCode PETSCDM_DLLEXPORT VecPackCreate(MPI_Comm comm,VecPack *packer)
     61: {
     63:   VecPack        p;

     67:   *packer = PETSC_NULL;
     68: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
     69:   DMInitializePackage(PETSC_NULL);
     70: #endif

     72:   PetscHeaderCreate(p,_p_VecPack,struct _VecPackOps,DA_COOKIE,0,"VecPack",comm,VecPackDestroy,0);
     73:   p->n            = 0;
     74:   p->next         = PETSC_NULL;
     75:   p->comm         = comm;
     76:   p->globalvector = PETSC_NULL;
     77:   p->nredundant   = 0;
     78:   p->nDA          = 0;
     79:   MPI_Comm_rank(comm,&p->rank);

     81:   p->ops->createglobalvector = VecPackCreateGlobalVector;
     82:   p->ops->refine             = VecPackRefine;
     83:   p->ops->getinterpolation   = VecPackGetInterpolation;
     84:   *packer = p;
     85:   return(0);
     86: }

     90: /*@C
     91:     VecPackDestroy - Destroys a vector packer.

     93:     Collective on VecPack

     95:     Input Parameter:
     96: .   packer - the packer object

     98:     Level: advanced

    100: .seealso VecPackCreate(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
    101:          VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()

    103: @*/
    104: PetscErrorCode PETSCDM_DLLEXPORT VecPackDestroy(VecPack packer)
    105: {
    106:   PetscErrorCode     ierr;
    107:   struct VecPackLink *next = packer->next,*prev;

    110:   if (--packer->refct > 0) return(0);
    111:   while (next) {
    112:     prev = next;
    113:     next = next->next;
    114:     if (prev->type == VECPACK_DA) {
    115:       DADestroy(prev->da);
    116:     }
    117:     PetscFree(prev);
    118:   }
    119:   if (packer->globalvector) {
    120:     VecDestroy(packer->globalvector);
    121:   }
    122:   PetscHeaderDestroy(packer);
    123:   return(0);
    124: }

    126: /* --------------------------------------------------------------------------------------*/

    130: PetscErrorCode VecPackGetAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
    131: {
    133:   PetscScalar    *varray;

    136:   if (array) {
    137:     if (!packer->rank) {
    138:       VecGetArray(vec,&varray);
    139:       *array  = varray + mine->rstart;
    140:       VecRestoreArray(vec,&varray);
    141:     } else {
    142:       *array = 0;
    143:     }
    144:   }
    145:   return(0);
    146: }

    150: PetscErrorCode VecPackGetAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
    151: {
    153:   PetscScalar    *array;

    156:   if (global) {
    157:     DAGetGlobalVector(mine->da,global);
    158:     VecGetArray(vec,&array);
    159:     VecPlaceArray(*global,array+mine->rstart);
    160:     VecRestoreArray(vec,&array);
    161:   }
    162:   return(0);
    163: }

    167: PetscErrorCode VecPackRestoreAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
    168: {
    170:   return(0);
    171: }

    175: PetscErrorCode VecPackRestoreAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
    176: {

    180:   if (global) {
    181:     VecResetArray(*global);
    182:     DARestoreGlobalVector(mine->da,global);
    183:   }
    184:   return(0);
    185: }

    189: PetscErrorCode VecPackScatter_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
    190: {
    192:   PetscScalar    *varray;

    196:   if (!packer->rank) {
    197:     VecGetArray(vec,&varray);
    198:     PetscMemcpy(array,varray+mine->rstart,mine->n*sizeof(PetscScalar));
    199:     VecRestoreArray(vec,&varray);
    200:   }
    201:   MPI_Bcast(array,mine->n,MPIU_SCALAR,0,packer->comm);
    202:   return(0);
    203: }

    207: PetscErrorCode VecPackScatter_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
    208: {
    210:   PetscScalar    *array;
    211:   Vec            global;

    214:   DAGetGlobalVector(mine->da,&global);
    215:   VecGetArray(vec,&array);
    216:   VecPlaceArray(global,array+mine->rstart);
    217:   DAGlobalToLocalBegin(mine->da,global,INSERT_VALUES,local);
    218:   DAGlobalToLocalEnd(mine->da,global,INSERT_VALUES,local);
    219:   VecRestoreArray(vec,&array);
    220:   VecResetArray(global);
    221:   DARestoreGlobalVector(mine->da,&global);
    222:   return(0);
    223: }

    227: PetscErrorCode VecPackGather_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
    228: {
    230:   PetscScalar    *varray;

    233:   if (!packer->rank) {
    234:     VecGetArray(vec,&varray);
    235:     if (varray+mine->rstart == array) SETERRQ(PETSC_ERR_ARG_WRONG,"You need not VecPackGather() into objects obtained via VecPackGetAccess()");
    236:     PetscMemcpy(varray+mine->rstart,array,mine->n*sizeof(PetscScalar));
    237:     VecRestoreArray(vec,&varray);
    238:   }
    239:   return(0);
    240: }

    244: PetscErrorCode VecPackGather_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
    245: {
    247:   PetscScalar    *array;
    248:   Vec            global;

    251:   DAGetGlobalVector(mine->da,&global);
    252:   VecGetArray(vec,&array);
    253:   VecPlaceArray(global,array+mine->rstart);
    254:   DALocalToGlobal(mine->da,local,INSERT_VALUES,global);
    255:   VecRestoreArray(vec,&array);
    256:   VecResetArray(global);
    257:   DARestoreGlobalVector(mine->da,&global);
    258:   return(0);
    259: }

    261: /* ----------------------------------------------------------------------------------*/

    263: #include <stdarg.h>

    267: /*@C
    268:     VecPackGetAccess - Allows one to access the individual packed vectors in their global
    269:        representation.

    271:     Collective on VecPack

    273:     Input Parameter:
    274: +    packer - the packer object
    275: .    gvec - the global vector
    276: -    ... - the individual sequential or parallel objects (arrays or vectors)
    277:
    278:     Level: advanced

    280: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    281:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
    282:          VecPackRestoreAccess(), VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    284: @*/
    285: PetscErrorCode PETSCDM_DLLEXPORT VecPackGetAccess(VecPack packer,Vec gvec,...)
    286: {
    287:   va_list            Argp;
    288:   PetscErrorCode     ierr;
    289:   struct VecPackLink *next = packer->next;

    292:   if (!packer->globalvector) {
    293:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
    294:   }

    296:   /* loop over packed objects, handling one at at time */
    297:   va_start(Argp,gvec);
    298:   while (next) {
    299:     if (next->type == VECPACK_ARRAY) {
    300:       PetscScalar **array;
    301:       array = va_arg(Argp, PetscScalar**);
    302:       VecPackGetAccess_Array(packer,next,gvec,array);
    303:     } else if (next->type == VECPACK_DA) {
    304:       Vec *vec;
    305:       vec  = va_arg(Argp, Vec*);
    306:       VecPackGetAccess_DA(packer,next,gvec,vec);
    307:     } else {
    308:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    309:     }
    310:     next = next->next;
    311:   }
    312:   va_end(Argp);
    313:   return(0);
    314: }

    318: /*@C
    319:     VecPackRestoreAccess - Allows one to access the individual packed vectors in their global
    320:        representation.

    322:     Collective on VecPack

    324:     Input Parameter:
    325: +    packer - the packer object
    326: .    gvec - the global vector
    327: -    ... - the individual sequential or parallel objects (arrays or vectors)
    328:
    329:     Level: advanced

    331: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    332:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
    333:          VecPackRestoreAccess()

    335: @*/
    336: PetscErrorCode PETSCDM_DLLEXPORT VecPackRestoreAccess(VecPack packer,Vec gvec,...)
    337: {
    338:   va_list            Argp;
    339:   PetscErrorCode     ierr;
    340:   struct VecPackLink *next = packer->next;

    343:   if (!packer->globalvector) {
    344:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
    345:   }

    347:   /* loop over packed objects, handling one at at time */
    348:   va_start(Argp,gvec);
    349:   while (next) {
    350:     if (next->type == VECPACK_ARRAY) {
    351:       PetscScalar **array;
    352:       array = va_arg(Argp, PetscScalar**);
    353:       VecPackRestoreAccess_Array(packer,next,gvec,array);
    354:     } else if (next->type == VECPACK_DA) {
    355:       Vec *vec;
    356:       vec  = va_arg(Argp, Vec*);
    357:       VecPackRestoreAccess_DA(packer,next,gvec,vec);
    358:     } else {
    359:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    360:     }
    361:     next = next->next;
    362:   }
    363:   va_end(Argp);
    364:   return(0);
    365: }

    369: /*@C
    370:     VecPackScatter - Scatters from a global packed vector into its individual local vectors

    372:     Collective on VecPack

    374:     Input Parameter:
    375: +    packer - the packer object
    376: .    gvec - the global vector
    377: -    ... - the individual sequential objects (arrays or vectors)
    378:
    379:     Level: advanced

    381: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    382:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    383:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    385: @*/
    386: PetscErrorCode PETSCDM_DLLEXPORT VecPackScatter(VecPack packer,Vec gvec,...)
    387: {
    388:   va_list            Argp;
    389:   PetscErrorCode     ierr;
    390:   struct VecPackLink *next = packer->next;

    393:   if (!packer->globalvector) {
    394:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
    395:   }

    397:   /* loop over packed objects, handling one at at time */
    398:   va_start(Argp,gvec);
    399:   while (next) {
    400:     if (next->type == VECPACK_ARRAY) {
    401:       PetscScalar *array;
    402:       array = va_arg(Argp, PetscScalar*);
    403:       VecPackScatter_Array(packer,next,gvec,array);
    404:     } else if (next->type == VECPACK_DA) {
    405:       Vec vec;
    406:       vec = va_arg(Argp, Vec);
    408:       VecPackScatter_DA(packer,next,gvec,vec);
    409:     } else {
    410:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    411:     }
    412:     next = next->next;
    413:   }
    414:   va_end(Argp);
    415:   return(0);
    416: }

    420: /*@C
    421:     VecPackGather - Gathers into a global packed vector from its individual local vectors

    423:     Collective on VecPack

    425:     Input Parameter:
    426: +    packer - the packer object
    427: .    gvec - the global vector
    428: -    ... - the individual sequential objects (arrays or vectors)
    429:
    430:     Level: advanced

    432: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    433:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    434:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    436: @*/
    437: PetscErrorCode PETSCDM_DLLEXPORT VecPackGather(VecPack packer,Vec gvec,...)
    438: {
    439:   va_list            Argp;
    440:   PetscErrorCode     ierr;
    441:   struct VecPackLink *next = packer->next;

    444:   if (!packer->globalvector) {
    445:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
    446:   }

    448:   /* loop over packed objects, handling one at at time */
    449:   va_start(Argp,gvec);
    450:   while (next) {
    451:     if (next->type == VECPACK_ARRAY) {
    452:       PetscScalar *array;
    453:       array = va_arg(Argp, PetscScalar*);
    454:       VecPackGather_Array(packer,next,gvec,array);
    455:     } else if (next->type == VECPACK_DA) {
    456:       Vec vec;
    457:       vec = va_arg(Argp, Vec);
    459:       VecPackGather_DA(packer,next,gvec,vec);
    460:     } else {
    461:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    462:     }
    463:     next = next->next;
    464:   }
    465:   va_end(Argp);
    466:   return(0);
    467: }

    471: /*@C
    472:     VecPackAddArray - adds an "redundant" array to a VecPack. The array values will
    473:        be stored in part of the array on processor 0.

    475:     Collective on VecPack

    477:     Input Parameter:
    478: +    packer - the packer object
    479: -    n - the length of the array
    480:
    481:     Level: advanced

    483: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
    484:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    485:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    487: @*/
    488: PetscErrorCode PETSCDM_DLLEXPORT VecPackAddArray(VecPack packer,PetscInt n)
    489: {
    490:   struct VecPackLink *mine,*next = packer->next;
    491:   PetscErrorCode     ierr;

    494:   if (packer->globalvector) {
    495:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add an array once you have called VecPackCreateGlobalVector()");
    496:   }

    498:   /* create new link */
    499:   PetscNew(struct VecPackLink,&mine);
    500:   mine->n             = n;
    501:   mine->da            = PETSC_NULL;
    502:   mine->type          = VECPACK_ARRAY;
    503:   mine->next          = PETSC_NULL;
    504:   if (!packer->rank) packer->n += n;

    506:   /* add to end of list */
    507:   if (!next) {
    508:     packer->next = mine;
    509:   } else {
    510:     while (next->next) next = next->next;
    511:     next->next = mine;
    512:   }
    513:   packer->nredundant++;
    514:   return(0);
    515: }

    519: /*@C
    520:     VecPackAddDA - adds a DA vector to a VecPack

    522:     Collective on VecPack

    524:     Input Parameter:
    525: +    packer - the packer object
    526: -    da - the DA object
    527:
    528:     Level: advanced

    530: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
    531:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    532:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    534: @*/
    535: PetscErrorCode PETSCDM_DLLEXPORT VecPackAddDA(VecPack packer,DA da)
    536: {
    537:   PetscErrorCode     ierr;
    538:   PetscInt           n;
    539:   struct VecPackLink *mine,*next = packer->next;
    540:   Vec                global;

    543:   if (packer->globalvector) {
    544:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DA once you have called VecPackCreateGlobalVector()");
    545:   }

    547:   /* create new link */
    548:   PetscNew(struct VecPackLink,&mine);
    549:   PetscObjectReference((PetscObject)da);
    550:   DAGetGlobalVector(da,&global);
    551:   VecGetLocalSize(global,&n);
    552:   DARestoreGlobalVector(da,&global);
    553:   mine->n      = n;
    554:   mine->da     = da;
    555:   mine->type   = VECPACK_DA;
    556:   mine->next   = PETSC_NULL;
    557:   packer->n   += n;

    559:   /* add to end of list */
    560:   if (!next) {
    561:     packer->next = mine;
    562:   } else {
    563:     while (next->next) next = next->next;
    564:     next->next = mine;
    565:   }
    566:   packer->nDA++;
    567:   return(0);
    568: }

    572: /*@C
    573:     VecPackCreateGlobalVector - Creates a vector of the correct size to be gathered into
    574:         by the packer.

    576:     Collective on VecPack

    578:     Input Parameter:
    579: .    packer - the packer object

    581:     Output Parameters:
    582: .   gvec - the global vector

    584:     Level: advanced

    586:     Notes: Once this has been created you cannot add additional arrays or vectors to be packed.

    588: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
    589:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    590:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    592: @*/
    593: PetscErrorCode PETSCDM_DLLEXPORT VecPackCreateGlobalVector(VecPack packer,Vec *gvec)
    594: {
    595:   PetscErrorCode     ierr;
    596:   PetscInt           nprev = 0;
    597:   PetscMPIInt        rank;
    598:   struct VecPackLink *next = packer->next;

    601:   if (packer->globalvector) {
    602:     VecDuplicate(packer->globalvector,gvec);
    603:   } else {
    604:     VecCreateMPI(packer->comm,packer->n,PETSC_DETERMINE,gvec);
    605:     PetscObjectReference((PetscObject)*gvec);
    606:     packer->globalvector = *gvec;

    608:     VecGetSize(*gvec,&packer->N);
    609:     VecGetOwnershipRange(*gvec,&packer->rstart,PETSC_NULL);
    610:
    611:     /* now set the rstart for each linked array/vector */
    612:     MPI_Comm_rank(packer->comm,&rank);
    613:     while (next) {
    614:       next->rstart = nprev;
    615:       if (!rank || next->type != VECPACK_ARRAY) nprev += next->n;
    616:       next = next->next;
    617:     }
    618:   }
    619:   return(0);
    620: }

    624: /*@C
    625:     VecPackGetGlobalIndices - Gets the global indices for all the entries in the packed
    626:       vectors.

    628:     Collective on VecPack

    630:     Input Parameter:
    631: .    packer - the packer object

    633:     Output Parameters:
    634: .    idx - the individual indices for each packed vector/array
    635:
    636:     Level: advanced

    638:     Notes:
    639:        The idx parameters should be freed by the calling routine with PetscFree()

    641: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    642:          VecPackGather(), VecPackCreate(), VecPackGetAccess(),
    643:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    645: @*/
    646: PetscErrorCode PETSCDM_DLLEXPORT VecPackGetGlobalIndices(VecPack packer,...)
    647: {
    648:   va_list            Argp;
    649:   PetscErrorCode     ierr;
    650:   PetscInt           i,**idx,n;
    651:   struct VecPackLink *next = packer->next;
    652:   Vec                global,dglobal;
    653:   PF                 pf;
    654:   PetscScalar        *array;

    657:   VecPackCreateGlobalVector(packer,&global);

    659:   /* put 0 to N-1 into the global vector */
    660:   PFCreate(PETSC_COMM_WORLD,1,1,&pf);
    661:   PFSetType(pf,PFIDENTITY,PETSC_NULL);
    662:   PFApplyVec(pf,PETSC_NULL,global);
    663:   PFDestroy(pf);

    665:   /* loop over packed objects, handling one at at time */
    666:   va_start(Argp,packer);
    667:   while (next) {
    668:     idx = va_arg(Argp, PetscInt**);

    670:     if (next->type == VECPACK_ARRAY) {
    671:
    672:       PetscMalloc(next->n*sizeof(PetscInt),idx);
    673:       if (!packer->rank) {
    674:         VecGetArray(global,&array);
    675:         array += next->rstart;
    676:         for (i=0; i<next->n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
    677:         array -= next->rstart;
    678:         VecRestoreArray(global,&array);
    679:       }
    680:       MPI_Bcast(*idx,next->n,MPIU_INT,0,packer->comm);

    682:     } else if (next->type == VECPACK_DA) {
    683:       Vec local;

    685:       DACreateLocalVector(next->da,&local);
    686:       VecGetArray(global,&array);
    687:       array += next->rstart;
    688:       DAGetGlobalVector(next->da,&dglobal);
    689:       VecPlaceArray(dglobal,array);
    690:       DAGlobalToLocalBegin(next->da,dglobal,INSERT_VALUES,local);
    691:       DAGlobalToLocalEnd(next->da,dglobal,INSERT_VALUES,local);
    692:       array -= next->rstart;
    693:       VecRestoreArray(global,&array);
    694:       VecResetArray(dglobal);
    695:       DARestoreGlobalVector(next->da,&dglobal);

    697:       VecGetArray(local,&array);
    698:       VecGetSize(local,&n);
    699:       PetscMalloc(n*sizeof(PetscInt),idx);
    700:       for (i=0; i<n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
    701:       VecRestoreArray(local,&array);
    702:       VecDestroy(local);

    704:     } else {
    705:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    706:     }
    707:     next = next->next;
    708:   }
    709:   va_end(Argp);
    710:   VecDestroy(global);
    711:   return(0);
    712: }

    714: /* -------------------------------------------------------------------------------------*/
    717: PetscErrorCode VecPackGetLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
    718: {
    721:   PetscMalloc(mine->n*sizeof(PetscScalar),array);
    722:   return(0);
    723: }

    727: PetscErrorCode VecPackGetLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
    728: {
    731:   DAGetLocalVector(mine->da,local);
    732:   return(0);
    733: }

    737: PetscErrorCode VecPackRestoreLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
    738: {
    741:   PetscFree(*array);
    742:   return(0);
    743: }

    747: PetscErrorCode VecPackRestoreLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
    748: {
    751:   DARestoreLocalVector(mine->da,local);
    752:   return(0);
    753: }

    757: /*@C
    758:     VecPackGetLocalVectors - Gets local vectors and arrays for each part of a VecPack.'
    759:        Use VecPakcRestoreLocalVectors() to return them.

    761:     Collective on VecPack

    763:     Input Parameter:
    764: .    packer - the packer object
    765:
    766:     Output Parameter:
    767: .    ... - the individual sequential objects (arrays or vectors)
    768:
    769:     Level: advanced

    771: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    772:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    773:          VecPackRestoreLocalVectors()

    775: @*/
    776: PetscErrorCode PETSCDM_DLLEXPORT VecPackGetLocalVectors(VecPack packer,...)
    777: {
    778:   va_list            Argp;
    779:   PetscErrorCode     ierr;
    780:   struct VecPackLink *next = packer->next;

    784:   /* loop over packed objects, handling one at at time */
    785:   va_start(Argp,packer);
    786:   while (next) {
    787:     if (next->type == VECPACK_ARRAY) {
    788:       PetscScalar **array;
    789:       array = va_arg(Argp, PetscScalar**);
    790:       VecPackGetLocalVectors_Array(packer,next,array);
    791:     } else if (next->type == VECPACK_DA) {
    792:       Vec *vec;
    793:       vec = va_arg(Argp, Vec*);
    794:       VecPackGetLocalVectors_DA(packer,next,vec);
    795:     } else {
    796:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    797:     }
    798:     next = next->next;
    799:   }
    800:   va_end(Argp);
    801:   return(0);
    802: }

    806: /*@C
    807:     VecPackRestoreLocalVectors - Restores local vectors and arrays for each part of a VecPack.'
    808:        Use VecPakcRestoreLocalVectors() to return them.

    810:     Collective on VecPack

    812:     Input Parameter:
    813: .    packer - the packer object
    814:
    815:     Output Parameter:
    816: .    ... - the individual sequential objects (arrays or vectors)
    817:
    818:     Level: advanced

    820: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    821:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    822:          VecPackGetLocalVectors()

    824: @*/
    825: PetscErrorCode PETSCDM_DLLEXPORT VecPackRestoreLocalVectors(VecPack packer,...)
    826: {
    827:   va_list            Argp;
    828:   PetscErrorCode     ierr;
    829:   struct VecPackLink *next = packer->next;

    833:   /* loop over packed objects, handling one at at time */
    834:   va_start(Argp,packer);
    835:   while (next) {
    836:     if (next->type == VECPACK_ARRAY) {
    837:       PetscScalar **array;
    838:       array = va_arg(Argp, PetscScalar**);
    839:       VecPackRestoreLocalVectors_Array(packer,next,array);
    840:     } else if (next->type == VECPACK_DA) {
    841:       Vec *vec;
    842:       vec = va_arg(Argp, Vec*);
    843:       VecPackRestoreLocalVectors_DA(packer,next,vec);
    844:     } else {
    845:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    846:     }
    847:     next = next->next;
    848:   }
    849:   va_end(Argp);
    850:   return(0);
    851: }

    853: /* -------------------------------------------------------------------------------------*/
    856: PetscErrorCode VecPackGetEntries_Array(VecPack packer,struct VecPackLink *mine,PetscInt *n)
    857: {
    859:   if (n) *n = mine->n;
    860:   return(0);
    861: }

    865: PetscErrorCode VecPackGetEntries_DA(VecPack packer,struct VecPackLink *mine,DA *da)
    866: {
    868:   if (da) *da = mine->da;
    869:   return(0);
    870: }

    874: /*@C
    875:     VecPackGetEntries - Gets the DA, redundant size, etc for each entry in a VecPack.
    876:        Use VecPackRestoreEntries() to return them.

    878:     Collective on VecPack

    880:     Input Parameter:
    881: .    packer - the packer object
    882:
    883:     Output Parameter:
    884: .    ... - the individual entries, DAs or integer sizes)
    885:
    886:     Level: advanced

    888: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    889:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    890:          VecPackRestoreLocalVectors(), VecPackGetLocalVectors(), VecPackRestoreEntries(),
    891:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    893: @*/
    894: PetscErrorCode PETSCDM_DLLEXPORT VecPackGetEntries(VecPack packer,...)
    895: {
    896:   va_list            Argp;
    897:   PetscErrorCode     ierr;
    898:   struct VecPackLink *next = packer->next;

    902:   /* loop over packed objects, handling one at at time */
    903:   va_start(Argp,packer);
    904:   while (next) {
    905:     if (next->type == VECPACK_ARRAY) {
    906:       PetscInt *n;
    907:       n = va_arg(Argp, PetscInt*);
    908:       VecPackGetEntries_Array(packer,next,n);
    909:     } else if (next->type == VECPACK_DA) {
    910:       DA *da;
    911:       da = va_arg(Argp, DA*);
    912:       VecPackGetEntries_DA(packer,next,da);
    913:     } else {
    914:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    915:     }
    916:     next = next->next;
    917:   }
    918:   va_end(Argp);
    919:   return(0);
    920: }

    924: /*@C
    925:     VecPackRefine - Refines a VecPack by refining all of its DAs

    927:     Collective on VecPack

    929:     Input Parameters:
    930: +    packer - the packer object
    931: -    comm - communicator to contain the new DM object, usually PETSC_NULL

    933:     Output Parameter:
    934: .    fine - new packer
    935:
    936:     Level: advanced

    938: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    939:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    940:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    942: @*/
    943: PetscErrorCode PETSCDM_DLLEXPORT VecPackRefine(VecPack packer,MPI_Comm comm,VecPack *fine)
    944: {
    945:   PetscErrorCode     ierr;
    946:   struct VecPackLink *next = packer->next;
    947:   DA                 da;

    950:   VecPackCreate(comm,fine);

    952:   /* loop over packed objects, handling one at at time */
    953:   while (next) {
    954:     if (next->type == VECPACK_ARRAY) {
    955:       VecPackAddArray(*fine,next->n);
    956:     } else if (next->type == VECPACK_DA) {
    957:       DARefine(next->da,comm,&da);
    958:       VecPackAddDA(*fine,da);
    959:       PetscObjectDereference((PetscObject)da);
    960:     } else {
    961:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    962:     }
    963:     next = next->next;
    964:   }
    965:   return(0);
    966: }

    968:  #include petscmat.h

    970: struct MatPackLink {
    971:   Mat                A;
    972:   struct MatPackLink *next;
    973: };

    975: struct MatPack {
    976:   VecPack            right,left;
    977:   struct MatPackLink *next;
    978: };

    982: PetscErrorCode MatMultBoth_Shell_Pack(Mat A,Vec x,Vec y,PetscTruth add)
    983: {
    984:   struct MatPack     *mpack;
    985:   struct VecPackLink *xnext,*ynext;
    986:   struct MatPackLink *anext;
    987:   PetscScalar        *xarray,*yarray;
    988:   PetscErrorCode     ierr;
    989:   PetscInt           i;
    990:   Vec                xglobal,yglobal;

    993:   MatShellGetContext(A,(void**)&mpack);
    994:   xnext = mpack->right->next;
    995:   ynext = mpack->left->next;
    996:   anext = mpack->next;

    998:   while (xnext) {
    999:     if (xnext->type == VECPACK_ARRAY) {
    1000:       if (!mpack->right->rank) {
    1001:         VecGetArray(x,&xarray);
    1002:         VecGetArray(y,&yarray);
    1003:         if (add) {
    1004:           for (i=0; i<xnext->n; i++) {
    1005:             yarray[ynext->rstart+i] += xarray[xnext->rstart+i];
    1006:           }
    1007:         } else {
    1008:           PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
    1009:         }
    1010:         VecRestoreArray(x,&xarray);
    1011:         VecRestoreArray(y,&yarray);
    1012:       }
    1013:     } else if (xnext->type == VECPACK_DA) {
    1014:       VecGetArray(x,&xarray);
    1015:       VecGetArray(y,&yarray);
    1016:       DAGetGlobalVector(xnext->da,&xglobal);
    1017:       DAGetGlobalVector(ynext->da,&yglobal);
    1018:       VecPlaceArray(xglobal,xarray+xnext->rstart);
    1019:       VecPlaceArray(yglobal,yarray+ynext->rstart);
    1020:       if (add) {
    1021:         MatMultAdd(anext->A,xglobal,yglobal,yglobal);
    1022:       } else {
    1023:         MatMult(anext->A,xglobal,yglobal);
    1024:       }
    1025:       VecRestoreArray(x,&xarray);
    1026:       VecRestoreArray(y,&yarray);
    1027:       VecResetArray(xglobal);
    1028:       VecResetArray(yglobal);
    1029:       DARestoreGlobalVector(xnext->da,&xglobal);
    1030:       DARestoreGlobalVector(ynext->da,&yglobal);
    1031:       anext = anext->next;
    1032:     } else {
    1033:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    1034:     }
    1035:     xnext = xnext->next;
    1036:     ynext = ynext->next;
    1037:   }
    1038:   return(0);
    1039: }

    1043: PetscErrorCode MatMultAdd_Shell_Pack(Mat A,Vec x,Vec y,Vec z)
    1044: {
    1047:   if (z != y) SETERRQ(PETSC_ERR_SUP,"Handles y == z only");
    1048:   MatMultBoth_Shell_Pack(A,x,y,PETSC_TRUE);
    1049:   return(0);
    1050: }

    1054: PetscErrorCode MatMult_Shell_Pack(Mat A,Vec x,Vec y)
    1055: {
    1058:   MatMultBoth_Shell_Pack(A,x,y,PETSC_FALSE);
    1059:   return(0);
    1060: }

    1064: PetscErrorCode MatMultTranspose_Shell_Pack(Mat A,Vec x,Vec y)
    1065: {
    1066:   struct MatPack     *mpack;
    1067:   struct VecPackLink *xnext,*ynext;
    1068:   struct MatPackLink *anext;
    1069:   PetscScalar        *xarray,*yarray;
    1070:   PetscErrorCode     ierr;
    1071:   Vec                xglobal,yglobal;

    1074:   MatShellGetContext(A,(void**)&mpack);
    1075:   xnext = mpack->left->next;
    1076:   ynext = mpack->right->next;
    1077:   anext = mpack->next;

    1079:   while (xnext) {
    1080:     if (xnext->type == VECPACK_ARRAY) {
    1081:       if (!mpack->right->rank) {
    1082:         VecGetArray(x,&xarray);
    1083:         VecGetArray(y,&yarray);
    1084:         PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
    1085:         VecRestoreArray(x,&xarray);
    1086:         VecRestoreArray(y,&yarray);
    1087:       }
    1088:     } else if (xnext->type == VECPACK_DA) {
    1089:       VecGetArray(x,&xarray);
    1090:       VecGetArray(y,&yarray);
    1091:       DAGetGlobalVector(xnext->da,&xglobal);
    1092:       DAGetGlobalVector(ynext->da,&yglobal);
    1093:       VecPlaceArray(xglobal,xarray+xnext->rstart);
    1094:       VecPlaceArray(yglobal,yarray+ynext->rstart);
    1095:       MatMultTranspose(anext->A,xglobal,yglobal);
    1096:       VecRestoreArray(x,&xarray);
    1097:       VecRestoreArray(y,&yarray);
    1098:       VecResetArray(xglobal);
    1099:       VecResetArray(yglobal);
    1100:       DARestoreGlobalVector(xnext->da,&xglobal);
    1101:       DARestoreGlobalVector(ynext->da,&yglobal);
    1102:       anext = anext->next;
    1103:     } else {
    1104:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    1105:     }
    1106:     xnext = xnext->next;
    1107:     ynext = ynext->next;
    1108:   }
    1109:   return(0);
    1110: }

    1114: PetscErrorCode MatDestroy_Shell_Pack(Mat A)
    1115: {
    1116:   struct MatPack     *mpack;
    1117:   struct MatPackLink *anext,*oldanext;
    1118:   PetscErrorCode     ierr;

    1121:   MatShellGetContext(A,(void**)&mpack);
    1122:   anext = mpack->next;

    1124:   while (anext) {
    1125:     MatDestroy(anext->A);
    1126:     oldanext = anext;
    1127:     anext    = anext->next;
    1128:     PetscFree(oldanext);
    1129:   }
    1130:   PetscFree(mpack);
    1131:   return(0);
    1132: }

    1136: /*@C
    1137:     VecPackGetInterpolation - GetInterpolations a VecPack by refining all of its DAs

    1139:     Collective on VecPack

    1141:     Input Parameters:
    1142: +    coarse - coarse grid packer
    1143: -    fine - fine grid packer

    1145:     Output Parameter:
    1146: +    A - interpolation matrix
    1147: -    v - scaling vector
    1148:
    1149:     Level: advanced

    1151: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
    1152:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
    1153:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

    1155: @*/
    1156: PetscErrorCode PETSCDM_DLLEXPORT VecPackGetInterpolation(VecPack coarse,VecPack fine,Mat *A,Vec *v)
    1157: {
    1158:   PetscErrorCode     ierr;
    1159:   PetscInt           m,n,M,N;
    1160:   struct VecPackLink *nextc  = coarse->next;
    1161:   struct VecPackLink *nextf = fine->next;
    1162:   struct MatPackLink *nextmat,*pnextmat = 0;
    1163:   struct MatPack     *mpack;
    1164:   Vec                gcoarse,gfine;

    1167:   /* use global vectors only for determining matrix layout */
    1168:   VecPackCreateGlobalVector(coarse,&gcoarse);
    1169:   VecPackCreateGlobalVector(fine,&gfine);
    1170:   VecGetLocalSize(gcoarse,&n);
    1171:   VecGetLocalSize(gfine,&m);
    1172:   VecGetSize(gcoarse,&N);
    1173:   VecGetSize(gfine,&M);
    1174:   VecDestroy(gcoarse);
    1175:   VecDestroy(gfine);

    1177:   PetscNew(struct MatPack,&mpack);
    1178:   mpack->right = coarse;
    1179:   mpack->left  = fine;
    1180:   MatCreate(fine->comm,A);
    1181:   MatSetSizes(*A,m,n,M,N);
    1182:   MatSetType(*A,MATSHELL);
    1183:   MatShellSetContext(*A,mpack);
    1184:   MatShellSetOperation(*A,MATOP_MULT,(void(*)(void))MatMult_Shell_Pack);
    1185:   MatShellSetOperation(*A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell_Pack);
    1186:   MatShellSetOperation(*A,MATOP_MULT_ADD,(void(*)(void))MatMultAdd_Shell_Pack);
    1187:   MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell_Pack);

    1189:   /* loop over packed objects, handling one at at time */
    1190:   while (nextc) {
    1191:     if (nextc->type != nextf->type) SETERRQ(PETSC_ERR_ARG_INCOMP,"Two VecPack have different layout");

    1193:     if (nextc->type == VECPACK_ARRAY) {
    1194:       ;
    1195:     } else if (nextc->type == VECPACK_DA) {
    1196:       PetscNew(struct MatPackLink,&nextmat);
    1197:       nextmat->next = 0;
    1198:       if (pnextmat) {
    1199:         pnextmat->next = nextmat;
    1200:         pnextmat       = nextmat;
    1201:       } else {
    1202:         pnextmat    = nextmat;
    1203:         mpack->next = nextmat;
    1204:       }
    1205:       DAGetInterpolation(nextc->da,nextf->da,&nextmat->A,PETSC_NULL);
    1206:     } else {
    1207:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
    1208:     }
    1209:     nextc = nextc->next;
    1210:     nextf = nextf->next;
    1211:   }
    1212:   return(0);
    1213: }