Actual source code: pack.c

  1: #define PETSCDM_DLL
  2: 
 3:  #include petscda.h
 4:  #include src/dm/da/daimpl.h
 5:  #include petscmat.h

  7: typedef struct _VecPackOps *VecPackOps;
  8: struct _VecPackOps {
  9:   PetscErrorCode (*view)(VecPack,PetscViewer);
 10:   PetscErrorCode (*createglobalvector)(VecPack,Vec*);
 11:   PetscErrorCode (*getcoloring)(VecPack,ISColoringType,ISColoring*);
 12:   PetscErrorCode (*getmatrix)(VecPack, MatType,Mat*);
 13:   PetscErrorCode (*getinterpolation)(VecPack,VecPack,Mat*,Vec*);
 14:   PetscErrorCode (*refine)(VecPack,MPI_Comm,VecPack*);
 15:   PetscErrorCode (*getinjection)(VecPack,VecPack,VecScatter*);

 17:   PetscErrorCode (*forminitialguess)(VecPack,PetscErrorCode (*)(void),Vec,void*);
 18:   PetscErrorCode (*formfunction)(VecPack,PetscErrorCode (*)(void),Vec,void*);
 19: };

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

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

 28: struct VecPackLink {
 29:   VecPackLinkType    type;
 30:   struct VecPackLink *next;
 31:   PetscInt           n,rstart;      /* rstart is relative to this processor */

 33:   /* only used for VECPACK_DA */
 34:   PetscInt           *grstarts;     /* global row for first unknown of this DA on each process */
 35:   DA                 da;

 37:   /* only used for VECPACK_ARRAY */
 38:   PetscInt           grstart;        /* global row for first array unknown */
 39:   PetscMPIInt        rank;          /* process where array unknowns live */
 40: };

 42: struct _p_VecPack {
 43:   PETSCHEADER(struct _VecPackOps);
 44:   PetscInt           n,N,rstart;     /* rstart is relative to all processors, n unknowns owned by this process, N is total unknowns */
 45:   PetscInt           nghost;         /* number of all local entries include DA ghost points and any shared redundant arrays */
 46:   Vec                globalvector;
 47:   PetscInt           nDA,nredundant; /* how many DA's and seperate redundant arrays used to build VecPack */
 48:   struct VecPackLink *next;
 49: };

 53: /*
 54:     Maps from 
 55: */
 56: PetscErrorCode  VecPackFormInitialGuess_DADADA(VecPack pack,PetscErrorCode (*fun)(void),Vec X,void *ctx)
 57: {
 59:   PetscErrorCode (*f)(DALocalInfo*,void*,DALocalInfo*,void*,DALocalInfo*,void*) =
 60:                            (PetscErrorCode (*)(DALocalInfo*,void*,DALocalInfo*,void*,DALocalInfo*,void*)) fun;
 61:   DALocalInfo da1,da2,da3;
 62:   DA          DA1,DA2,DA3;
 63:   void        *x1,*x2,*x3;
 64:   Vec         X1,X2,X3;

 67:   VecPackGetEntries(pack,&DA1,&DA2,&DA3);
 68:   DAGetLocalInfo(DA1,&da1);
 69:   DAGetLocalInfo(DA2,&da2);
 70:   DAGetLocalInfo(DA3,&da3);
 71:   VecPackGetAccess(pack,X,&X1,&X2,&X3);
 72:   DAVecGetArray(DA1,X1,&x1);
 73:   DAVecGetArray(DA2,X2,&x2);
 74:   DAVecGetArray(DA3,X3,&x3);


 77:   (*f)(&da1,x1,&da2,x2,&da3,x3);

 79:   return(0);
 80: }

 84: PetscErrorCode VecView_DADADA(VecPack pack,Vec X,PetscViewer viewer)
 85: {
 87:   DA             DA1,DA2,DA3;
 88:   Vec            X1,X2,X3;

 91:   VecPackGetEntries(pack,&DA1,&DA2,&DA3);
 92:   VecPackGetAccess(pack,X,&X1,&X2,&X3);
 93:   VecView(X1,viewer);
 94:   VecView(X2,viewer);
 95:   VecView(X3,viewer);
 96:   return(0);
 97: }

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

105:     Collective on MPI_Comm

107:     Input Parameter:
108: .   comm - the processors that will share the global vector

110:     Output Parameters:
111: .   packer - the packer object

113:     Level: advanced

115: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
116:          VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()
117:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

119: @*/
120: PetscErrorCode  VecPackCreate(MPI_Comm comm,VecPack *packer)
121: {
123:   VecPack        p;

127:   *packer = PETSC_NULL;
128: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
129:   DMInitializePackage(PETSC_NULL);
130: #endif

132:   PetscHeaderCreate(p,_p_VecPack,struct _VecPackOps,DA_COOKIE,0,"VecPack",comm,VecPackDestroy,0);
133:   p->n            = 0;
134:   p->next         = PETSC_NULL;
135:   p->comm         = comm;
136:   p->globalvector = PETSC_NULL;
137:   p->nredundant   = 0;
138:   p->nDA          = 0;

140:   p->ops->createglobalvector = VecPackCreateGlobalVector;
141:   p->ops->refine             = VecPackRefine;
142:   p->ops->getinterpolation   = VecPackGetInterpolation;
143:   p->ops->getmatrix          = VecPackGetMatrix;
144:   p->ops->getcoloring        = VecPackGetColoring;

146:   p->ops->forminitialguess   = VecPackFormInitialGuess_DADADA;
147:   *packer = p;
148:   return(0);
149: }

153: /*@C
154:     VecPackDestroy - Destroys a vector packer.

156:     Collective on VecPack

158:     Input Parameter:
159: .   packer - the packer object

161:     Level: advanced

163: .seealso VecPackCreate(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
164:          VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()

166: @*/
167: PetscErrorCode  VecPackDestroy(VecPack packer)
168: {
169:   PetscErrorCode     ierr;
170:   struct VecPackLink *next = packer->next,*prev;

173:   if (--packer->refct > 0) return(0);
174:   while (next) {
175:     prev = next;
176:     next = next->next;
177:     if (prev->type == VECPACK_DA) {
178:       DADestroy(prev->da);
179:     }
180:     if (prev->grstarts) {
181:       PetscFree(prev->grstarts);
182:     }
183:     PetscFree(prev);
184:   }
185:   if (packer->globalvector) {
186:     VecDestroy(packer->globalvector);
187:   }
188:   PetscHeaderDestroy(packer);
189:   return(0);
190: }

192: /* --------------------------------------------------------------------------------------*/

196: PetscErrorCode VecPackGetAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
197: {
199:   PetscScalar    *varray;
200:   PetscMPIInt    rank;

203:   MPI_Comm_rank(packer->comm,&rank);
204:   if (array) {
205:     if (rank == mine->rank) {
206:       VecGetArray(vec,&varray);
207:       *array  = varray + mine->rstart;
208:       VecRestoreArray(vec,&varray);
209:     } else {
210:       *array = 0;
211:     }
212:   }
213:   return(0);
214: }

218: PetscErrorCode VecPackGetAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
219: {
221:   PetscScalar    *array;

224:   if (global) {
225:     DAGetGlobalVector(mine->da,global);
226:     VecGetArray(vec,&array);
227:     VecPlaceArray(*global,array+mine->rstart);
228:     VecRestoreArray(vec,&array);
229:   }
230:   return(0);
231: }

235: PetscErrorCode VecPackRestoreAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
236: {
238:   return(0);
239: }

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

248:   if (global) {
249:     VecResetArray(*global);
250:     DARestoreGlobalVector(mine->da,global);
251:   }
252:   return(0);
253: }

257: PetscErrorCode VecPackScatter_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
258: {
260:   PetscScalar    *varray;
261:   PetscMPIInt    rank;

264:   MPI_Comm_rank(packer->comm,&rank);
265:   if (rank == mine->rank) {
266:     VecGetArray(vec,&varray);
267:     PetscMemcpy(array,varray+mine->rstart,mine->n*sizeof(PetscScalar));
268:     VecRestoreArray(vec,&varray);
269:   }
270:   MPI_Bcast(array,mine->n,MPIU_SCALAR,mine->rank,packer->comm);
271:   return(0);
272: }

276: PetscErrorCode VecPackScatter_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
277: {
279:   PetscScalar    *array;
280:   Vec            global;

283:   DAGetGlobalVector(mine->da,&global);
284:   VecGetArray(vec,&array);
285:   VecPlaceArray(global,array+mine->rstart);
286:   DAGlobalToLocalBegin(mine->da,global,INSERT_VALUES,local);
287:   DAGlobalToLocalEnd(mine->da,global,INSERT_VALUES,local);
288:   VecRestoreArray(vec,&array);
289:   VecResetArray(global);
290:   DARestoreGlobalVector(mine->da,&global);
291:   return(0);
292: }

296: PetscErrorCode VecPackGather_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
297: {
299:   PetscScalar    *varray;
300:   PetscMPIInt    rank;

303:   MPI_Comm_rank(packer->comm,&rank);
304:   if (rank == mine->rank) {
305:     VecGetArray(vec,&varray);
306:     if (varray+mine->rstart == array) SETERRQ(PETSC_ERR_ARG_WRONG,"You need not VecPackGather() into objects obtained via VecPackGetAccess()");
307:     PetscMemcpy(varray+mine->rstart,array,mine->n*sizeof(PetscScalar));
308:     VecRestoreArray(vec,&varray);
309:   }
310:   return(0);
311: }

315: PetscErrorCode VecPackGather_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
316: {
318:   PetscScalar    *array;
319:   Vec            global;

322:   DAGetGlobalVector(mine->da,&global);
323:   VecGetArray(vec,&array);
324:   VecPlaceArray(global,array+mine->rstart);
325:   DALocalToGlobal(mine->da,local,INSERT_VALUES,global);
326:   VecRestoreArray(vec,&array);
327:   VecResetArray(global);
328:   DARestoreGlobalVector(mine->da,&global);
329:   return(0);
330: }

332: /* ----------------------------------------------------------------------------------*/

334: #include <stdarg.h>

338: /*@C
339:     VecPackGetAccess - Allows one to access the individual packed vectors in their global
340:        representation.

342:     Collective on VecPack

344:     Input Parameter:
345: +    packer - the packer object
346: .    gvec - the global vector
347: -    ... - the individual sequential or parallel objects (arrays or vectors)
348:  
349:     Level: advanced

351: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
352:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
353:          VecPackRestoreAccess(), VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

355: @*/
356: PetscErrorCode  VecPackGetAccess(VecPack packer,Vec gvec,...)
357: {
358:   va_list            Argp;
359:   PetscErrorCode     ierr;
360:   struct VecPackLink *next = packer->next;

363:   if (!packer->globalvector) {
364:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
365:   }

367:   /* loop over packed objects, handling one at at time */
368:   va_start(Argp,gvec);
369:   while (next) {
370:     if (next->type == VECPACK_ARRAY) {
371:       PetscScalar **array;
372:       array = va_arg(Argp, PetscScalar**);
373:       VecPackGetAccess_Array(packer,next,gvec,array);
374:     } else if (next->type == VECPACK_DA) {
375:       Vec *vec;
376:       vec  = va_arg(Argp, Vec*);
377:       VecPackGetAccess_DA(packer,next,gvec,vec);
378:     } else {
379:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
380:     }
381:     next = next->next;
382:   }
383:   va_end(Argp);
384:   return(0);
385: }

389: /*@C
390:     VecPackRestoreAccess - Allows one to access the individual packed vectors in their global
391:        representation.

393:     Collective on VecPack

395:     Input Parameter:
396: +    packer - the packer object
397: .    gvec - the global vector
398: -    ... - the individual sequential or parallel objects (arrays or vectors)
399:  
400:     Level: advanced

402: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
403:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
404:          VecPackRestoreAccess()

406: @*/
407: PetscErrorCode  VecPackRestoreAccess(VecPack packer,Vec gvec,...)
408: {
409:   va_list            Argp;
410:   PetscErrorCode     ierr;
411:   struct VecPackLink *next = packer->next;

414:   if (!packer->globalvector) {
415:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
416:   }

418:   /* loop over packed objects, handling one at at time */
419:   va_start(Argp,gvec);
420:   while (next) {
421:     if (next->type == VECPACK_ARRAY) {
422:       PetscScalar **array;
423:       array = va_arg(Argp, PetscScalar**);
424:       VecPackRestoreAccess_Array(packer,next,gvec,array);
425:     } else if (next->type == VECPACK_DA) {
426:       Vec *vec;
427:       vec  = va_arg(Argp, Vec*);
428:       VecPackRestoreAccess_DA(packer,next,gvec,vec);
429:     } else {
430:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
431:     }
432:     next = next->next;
433:   }
434:   va_end(Argp);
435:   return(0);
436: }

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

443:     Collective on VecPack

445:     Input Parameter:
446: +    packer - the packer object
447: .    gvec - the global vector
448: -    ... - the individual sequential objects (arrays or vectors)
449:  
450:     Level: advanced

452: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
453:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
454:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

456: @*/
457: PetscErrorCode  VecPackScatter(VecPack packer,Vec gvec,...)
458: {
459:   va_list            Argp;
460:   PetscErrorCode     ierr;
461:   struct VecPackLink *next = packer->next;

464:   if (!packer->globalvector) {
465:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
466:   }

468:   /* loop over packed objects, handling one at at time */
469:   va_start(Argp,gvec);
470:   while (next) {
471:     if (next->type == VECPACK_ARRAY) {
472:       PetscScalar *array;
473:       array = va_arg(Argp, PetscScalar*);
474:       VecPackScatter_Array(packer,next,gvec,array);
475:     } else if (next->type == VECPACK_DA) {
476:       Vec vec;
477:       vec = va_arg(Argp, Vec);
479:       VecPackScatter_DA(packer,next,gvec,vec);
480:     } else {
481:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
482:     }
483:     next = next->next;
484:   }
485:   va_end(Argp);
486:   return(0);
487: }

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

494:     Collective on VecPack

496:     Input Parameter:
497: +    packer - the packer object
498: .    gvec - the global vector
499: -    ... - the individual sequential objects (arrays or vectors)
500:  
501:     Level: advanced

503: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
504:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
505:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

507: @*/
508: PetscErrorCode  VecPackGather(VecPack packer,Vec gvec,...)
509: {
510:   va_list            Argp;
511:   PetscErrorCode     ierr;
512:   struct VecPackLink *next = packer->next;

515:   if (!packer->globalvector) {
516:     SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
517:   }

519:   /* loop over packed objects, handling one at at time */
520:   va_start(Argp,gvec);
521:   while (next) {
522:     if (next->type == VECPACK_ARRAY) {
523:       PetscScalar *array;
524:       array = va_arg(Argp, PetscScalar*);
525:       VecPackGather_Array(packer,next,gvec,array);
526:     } else if (next->type == VECPACK_DA) {
527:       Vec vec;
528:       vec = va_arg(Argp, Vec);
530:       VecPackGather_DA(packer,next,gvec,vec);
531:     } else {
532:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
533:     }
534:     next = next->next;
535:   }
536:   va_end(Argp);
537:   return(0);
538: }

542: /*@C
543:     VecPackAddArray - adds an "redundant" array to a VecPack. The array values will 
544:        be stored in part of the array on process orank.

546:     Collective on VecPack

548:     Input Parameter:
549: +    packer - the packer object
550: .    orank - the process on which the array entries officially live, this number must be
551:              the same on all processes.
552: -    n - the length of the array
553:  
554:     Level: advanced

556: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
557:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
558:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

560: @*/
561: PetscErrorCode  VecPackAddArray(VecPack packer,PetscMPIInt orank,PetscInt n)
562: {
563:   struct VecPackLink *mine,*next = packer->next;
564:   PetscErrorCode     ierr;
565:   PetscMPIInt        rank;

568:   if (packer->globalvector) {
569:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add an array once you have called VecPackCreateGlobalVector()");
570:   }
571: #if defined(PETSC_USE_DEBUG)
572:   {
573:     PetscMPIInt        orankmax;
574:     MPI_Allreduce(&orank,&orankmax,1,MPI_INT,MPI_MAX,packer->comm);
575:     if (orank != orankmax) SETERRQ2(PETSC_ERR_ARG_INCOMP,"orank %d must be equal on all processes, another process has value %d",orank,orankmax);
576:   }
577: #endif

579:   MPI_Comm_rank(packer->comm,&rank);
580:   /* create new link */
581:   PetscNew(struct VecPackLink,&mine);
582:   mine->n             = n;
583:   mine->rank          = orank;
584:   mine->da            = PETSC_NULL;
585:   mine->type          = VECPACK_ARRAY;
586:   mine->next          = PETSC_NULL;
587:   if (rank == mine->rank) packer->n += n;

589:   /* add to end of list */
590:   if (!next) {
591:     packer->next = mine;
592:   } else {
593:     while (next->next) next = next->next;
594:     next->next = mine;
595:   }
596:   packer->nredundant++;
597:   return(0);
598: }

602: /*@C
603:     VecPackAddDA - adds a DA vector to a VecPack

605:     Collective on VecPack

607:     Input Parameter:
608: +    packer - the packer object
609: -    da - the DA object
610:  
611:     Level: advanced

613: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
614:          VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
615:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

617: @*/
618: PetscErrorCode  VecPackAddDA(VecPack packer,DA da)
619: {
620:   PetscErrorCode     ierr;
621:   PetscInt           n;
622:   struct VecPackLink *mine,*next = packer->next;
623:   Vec                global;

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

630:   /* create new link */
631:   PetscNew(struct VecPackLink,&mine);
632:   PetscObjectReference((PetscObject)da);
633:   DAGetGlobalVector(da,&global);
634:   VecGetLocalSize(global,&n);
635:   DARestoreGlobalVector(da,&global);
636:   mine->n      = n;
637:   mine->da     = da;
638:   mine->type   = VECPACK_DA;
639:   mine->next   = PETSC_NULL;
640:   packer->n   += n;

642:   /* add to end of list */
643:   if (!next) {
644:     packer->next = mine;
645:   } else {
646:     while (next->next) next = next->next;
647:     next->next = mine;
648:   }
649:   packer->nDA++;
650:   return(0);
651: }

655: /*@C
656:     VecPackCreateGlobalVector - Creates a vector of the correct size to be gathered into 
657:         by the packer.

659:     Collective on VecPack

661:     Input Parameter:
662: .    packer - the packer object

664:     Output Parameters:
665: .   gvec - the global vector

667:     Level: advanced

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

671:            This also serves as a VecPackFinalize(), perhaps that should be made a seperate 
672:            function?

674: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
675:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
676:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

678: @*/
679: PetscErrorCode  VecPackCreateGlobalVector(VecPack packer,Vec *gvec)
680: {
681:   PetscErrorCode     ierr;
682:   PetscInt           nprev = 0;
683:   PetscMPIInt        rank,size;
684:   struct VecPackLink *next = packer->next;

687:   if (packer->globalvector) {
688:     VecDuplicate(packer->globalvector,gvec);
689:   } else {
690:     VecCreateMPI(packer->comm,packer->n,PETSC_DETERMINE,gvec);
691:     PetscObjectReference((PetscObject)*gvec);
692:     packer->globalvector = *gvec;

694:     VecGetSize(*gvec,&packer->N);
695:     VecGetOwnershipRange(*gvec,&packer->rstart,PETSC_NULL);
696: 
697:     /* now set the rstart for each linked array/vector */
698:     MPI_Comm_rank(packer->comm,&rank);
699:     MPI_Comm_size(packer->comm,&size);
700:     while (next) {
701:       next->rstart = nprev;
702:       if ((rank == next->rank) || next->type != VECPACK_ARRAY) nprev += next->n;
703:       next->grstart = packer->rstart + next->rstart;
704:       if (next->type == VECPACK_ARRAY) {
705:         MPI_Bcast(&next->grstart,1,MPIU_INT,next->rank,packer->comm);
706:       } else {
707:         PetscMalloc(size*sizeof(PetscInt),&next->grstarts);
708:         MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,packer->comm);
709:       }
710:       next = next->next;
711:     }
712:     VecSetOperation(*gvec,VECOP_VIEW,(void(*)(void))VecView_DADADA);
713:   }
714:   return(0);
715: }

719: /*@C
720:     VecPackGetGlobalIndices - Gets the global indices for all the entries in the packed
721:       vectors.

723:     Collective on VecPack

725:     Input Parameter:
726: .    packer - the packer object

728:     Output Parameters:
729: .    idx - the individual indices for each packed vector/array. Note that this includes
730:            all the ghost points that individual ghosted DA's may have.
731:  
732:     Level: advanced

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

737: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
738:          VecPackGather(), VecPackCreate(), VecPackGetAccess(),
739:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

741: @*/
742: PetscErrorCode  VecPackGetGlobalIndices(VecPack packer,...)
743: {
744:   va_list            Argp;
745:   PetscErrorCode     ierr;
746:   PetscInt           i,**idx,n;
747:   struct VecPackLink *next = packer->next;
748:   Vec                global,dglobal;
749:   PF                 pf;
750:   PetscScalar        *array;
751:   PetscMPIInt        rank;

754:   VecPackCreateGlobalVector(packer,&global);
755:   MPI_Comm_rank(packer->comm,&rank);

757:   /* put 0 to N-1 into the global vector */
758:   PFCreate(PETSC_COMM_WORLD,1,1,&pf);
759:   PFSetType(pf,PFIDENTITY,PETSC_NULL);
760:   PFApplyVec(pf,PETSC_NULL,global);
761:   PFDestroy(pf);

763:   /* loop over packed objects, handling one at at time */
764:   va_start(Argp,packer);
765:   while (next) {
766:     idx = va_arg(Argp, PetscInt**);

768:     if (next->type == VECPACK_ARRAY) {
769: 
770:       PetscMalloc(next->n*sizeof(PetscInt),idx);
771:       if (rank == next->rank) {
772:         VecGetArray(global,&array);
773:         array += next->rstart;
774:         for (i=0; i<next->n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
775:         array -= next->rstart;
776:         VecRestoreArray(global,&array);
777:       }
778:       MPI_Bcast(*idx,next->n,MPIU_INT,next->rank,packer->comm);

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

783:       DACreateLocalVector(next->da,&local);
784:       VecGetArray(global,&array);
785:       array += next->rstart;
786:       DAGetGlobalVector(next->da,&dglobal);
787:       VecPlaceArray(dglobal,array);
788:       DAGlobalToLocalBegin(next->da,dglobal,INSERT_VALUES,local);
789:       DAGlobalToLocalEnd(next->da,dglobal,INSERT_VALUES,local);
790:       array -= next->rstart;
791:       VecRestoreArray(global,&array);
792:       VecResetArray(dglobal);
793:       DARestoreGlobalVector(next->da,&dglobal);

795:       VecGetArray(local,&array);
796:       VecGetSize(local,&n);
797:       PetscMalloc(n*sizeof(PetscInt),idx);
798:       for (i=0; i<n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
799:       VecRestoreArray(local,&array);
800:       VecDestroy(local);

802:     } else {
803:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
804:     }
805:     next = next->next;
806:   }
807:   va_end(Argp);
808:   VecDestroy(global);
809:   return(0);
810: }

812: /* -------------------------------------------------------------------------------------*/
815: PetscErrorCode VecPackGetLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
816: {
819:   PetscMalloc(mine->n*sizeof(PetscScalar),array);
820:   return(0);
821: }

825: PetscErrorCode VecPackGetLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
826: {
829:   DAGetLocalVector(mine->da,local);
830:   return(0);
831: }

835: PetscErrorCode VecPackRestoreLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
836: {
839:   PetscFree(*array);
840:   return(0);
841: }

845: PetscErrorCode VecPackRestoreLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
846: {
849:   DARestoreLocalVector(mine->da,local);
850:   return(0);
851: }

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

859:     Collective on VecPack

861:     Input Parameter:
862: .    packer - the packer object
863:  
864:     Output Parameter:
865: .    ... - the individual sequential objects (arrays or vectors)
866:  
867:     Level: advanced

869: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
870:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(), 
871:          VecPackRestoreLocalVectors()

873: @*/
874: PetscErrorCode  VecPackGetLocalVectors(VecPack packer,...)
875: {
876:   va_list            Argp;
877:   PetscErrorCode     ierr;
878:   struct VecPackLink *next = packer->next;


882:   /* loop over packed objects, handling one at at time */
883:   va_start(Argp,packer);
884:   while (next) {
885:     if (next->type == VECPACK_ARRAY) {
886:       PetscScalar **array;
887:       array = va_arg(Argp, PetscScalar**);
888:       VecPackGetLocalVectors_Array(packer,next,array);
889:     } else if (next->type == VECPACK_DA) {
890:       Vec *vec;
891:       vec = va_arg(Argp, Vec*);
892:       VecPackGetLocalVectors_DA(packer,next,vec);
893:     } else {
894:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
895:     }
896:     next = next->next;
897:   }
898:   va_end(Argp);
899:   return(0);
900: }

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

908:     Collective on VecPack

910:     Input Parameter:
911: .    packer - the packer object
912:  
913:     Output Parameter:
914: .    ... - the individual sequential objects (arrays or vectors)
915:  
916:     Level: advanced

918: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
919:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(), 
920:          VecPackGetLocalVectors()

922: @*/
923: PetscErrorCode  VecPackRestoreLocalVectors(VecPack packer,...)
924: {
925:   va_list            Argp;
926:   PetscErrorCode     ierr;
927:   struct VecPackLink *next = packer->next;


931:   /* loop over packed objects, handling one at at time */
932:   va_start(Argp,packer);
933:   while (next) {
934:     if (next->type == VECPACK_ARRAY) {
935:       PetscScalar **array;
936:       array = va_arg(Argp, PetscScalar**);
937:       VecPackRestoreLocalVectors_Array(packer,next,array);
938:     } else if (next->type == VECPACK_DA) {
939:       Vec *vec;
940:       vec = va_arg(Argp, Vec*);
941:       VecPackRestoreLocalVectors_DA(packer,next,vec);
942:     } else {
943:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
944:     }
945:     next = next->next;
946:   }
947:   va_end(Argp);
948:   return(0);
949: }

951: /* -------------------------------------------------------------------------------------*/
954: PetscErrorCode VecPackGetEntries_Array(VecPack packer,struct VecPackLink *mine,PetscInt *n)
955: {
957:   if (n) *n = mine->n;
958:   return(0);
959: }

963: PetscErrorCode VecPackGetEntries_DA(VecPack packer,struct VecPackLink *mine,DA *da)
964: {
966:   if (da) *da = mine->da;
967:   return(0);
968: }

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

976:     Collective on VecPack

978:     Input Parameter:
979: .    packer - the packer object
980:  
981:     Output Parameter:
982: .    ... - the individual entries, DAs or integer sizes)
983:  
984:     Level: advanced

986: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
987:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(), 
988:          VecPackRestoreLocalVectors(), VecPackGetLocalVectors(), VecPackRestoreEntries(),
989:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

991: @*/
992: PetscErrorCode  VecPackGetEntries(VecPack packer,...)
993: {
994:   va_list            Argp;
995:   PetscErrorCode     ierr;
996:   struct VecPackLink *next = packer->next;


1000:   /* loop over packed objects, handling one at at time */
1001:   va_start(Argp,packer);
1002:   while (next) {
1003:     if (next->type == VECPACK_ARRAY) {
1004:       PetscInt *n;
1005:       n = va_arg(Argp, PetscInt*);
1006:       VecPackGetEntries_Array(packer,next,n);
1007:     } else if (next->type == VECPACK_DA) {
1008:       DA *da;
1009:       da = va_arg(Argp, DA*);
1010:       VecPackGetEntries_DA(packer,next,da);
1011:     } else {
1012:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1013:     }
1014:     next = next->next;
1015:   }
1016:   va_end(Argp);
1017:   return(0);
1018: }

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

1025:     Collective on VecPack

1027:     Input Parameters:
1028: +    packer - the packer object
1029: -    comm - communicator to contain the new DM object, usually PETSC_NULL

1031:     Output Parameter:
1032: .    fine - new packer
1033:  
1034:     Level: advanced

1036: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
1037:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
1038:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

1040: @*/
1041: PetscErrorCode  VecPackRefine(VecPack packer,MPI_Comm comm,VecPack *fine)
1042: {
1043:   PetscErrorCode     ierr;
1044:   struct VecPackLink *next = packer->next;
1045:   DA                 da;

1048:   VecPackCreate(comm,fine);

1050:   /* loop over packed objects, handling one at at time */
1051:   while (next) {
1052:     if (next->type == VECPACK_ARRAY) {
1053:       VecPackAddArray(*fine,next->rank,next->n);
1054:     } else if (next->type == VECPACK_DA) {
1055:       DARefine(next->da,comm,&da);
1056:       VecPackAddDA(*fine,da);
1057:       PetscObjectDereference((PetscObject)da);
1058:     } else {
1059:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1060:     }
1061:     next = next->next;
1062:   }
1063:   return(0);
1064: }

1066:  #include petscmat.h

1068: struct MatPackLink {
1069:   Mat                A;
1070:   struct MatPackLink *next;
1071: };

1073: struct MatPack {
1074:   VecPack            right,left;
1075:   struct MatPackLink *next;
1076: };

1080: PetscErrorCode MatMultBoth_Shell_Pack(Mat A,Vec x,Vec y,PetscTruth add)
1081: {
1082:   struct MatPack     *mpack;
1083:   struct VecPackLink *xnext,*ynext;
1084:   struct MatPackLink *anext;
1085:   PetscScalar        *xarray,*yarray;
1086:   PetscErrorCode     ierr;
1087:   PetscInt           i;
1088:   Vec                xglobal,yglobal;
1089:   PetscMPIInt        rank;

1092:   MatShellGetContext(A,(void**)&mpack);
1093:   MPI_Comm_rank(mpack->right->comm,&rank);
1094:   xnext = mpack->right->next;
1095:   ynext = mpack->left->next;
1096:   anext = mpack->next;

1098:   while (xnext) {
1099:     if (xnext->type == VECPACK_ARRAY) {
1100:       if (rank == xnext->rank) {
1101:         VecGetArray(x,&xarray);
1102:         VecGetArray(y,&yarray);
1103:         if (add) {
1104:           for (i=0; i<xnext->n; i++) {
1105:             yarray[ynext->rstart+i] += xarray[xnext->rstart+i];
1106:           }
1107:         } else {
1108:           PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1109:         }
1110:         VecRestoreArray(x,&xarray);
1111:         VecRestoreArray(y,&yarray);
1112:       }
1113:     } else if (xnext->type == VECPACK_DA) {
1114:       VecGetArray(x,&xarray);
1115:       VecGetArray(y,&yarray);
1116:       DAGetGlobalVector(xnext->da,&xglobal);
1117:       DAGetGlobalVector(ynext->da,&yglobal);
1118:       VecPlaceArray(xglobal,xarray+xnext->rstart);
1119:       VecPlaceArray(yglobal,yarray+ynext->rstart);
1120:       if (add) {
1121:         MatMultAdd(anext->A,xglobal,yglobal,yglobal);
1122:       } else {
1123:         MatMult(anext->A,xglobal,yglobal);
1124:       }
1125:       VecRestoreArray(x,&xarray);
1126:       VecRestoreArray(y,&yarray);
1127:       VecResetArray(xglobal);
1128:       VecResetArray(yglobal);
1129:       DARestoreGlobalVector(xnext->da,&xglobal);
1130:       DARestoreGlobalVector(ynext->da,&yglobal);
1131:       anext = anext->next;
1132:     } else {
1133:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1134:     }
1135:     xnext = xnext->next;
1136:     ynext = ynext->next;
1137:   }
1138:   return(0);
1139: }

1143: PetscErrorCode MatMultAdd_Shell_Pack(Mat A,Vec x,Vec y,Vec z)
1144: {
1147:   if (z != y) SETERRQ(PETSC_ERR_SUP,"Handles y == z only");
1148:   MatMultBoth_Shell_Pack(A,x,y,PETSC_TRUE);
1149:   return(0);
1150: }

1154: PetscErrorCode MatMult_Shell_Pack(Mat A,Vec x,Vec y)
1155: {
1158:   MatMultBoth_Shell_Pack(A,x,y,PETSC_FALSE);
1159:   return(0);
1160: }

1164: PetscErrorCode MatMultTranspose_Shell_Pack(Mat A,Vec x,Vec y)
1165: {
1166:   struct MatPack     *mpack;
1167:   struct VecPackLink *xnext,*ynext;
1168:   struct MatPackLink *anext;
1169:   PetscScalar        *xarray,*yarray;
1170:   PetscErrorCode     ierr;
1171:   Vec                xglobal,yglobal;
1172:   PetscMPIInt        rank;

1175:   MatShellGetContext(A,(void**)&mpack);
1176:   MPI_Comm_rank(mpack->right->comm,&rank);
1177:   xnext = mpack->left->next;
1178:   ynext = mpack->right->next;
1179:   anext = mpack->next;

1181:   while (xnext) {
1182:     if (xnext->type == VECPACK_ARRAY) {
1183:       if (rank == ynext->rank) {
1184:         VecGetArray(x,&xarray);
1185:         VecGetArray(y,&yarray);
1186:         PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1187:         VecRestoreArray(x,&xarray);
1188:         VecRestoreArray(y,&yarray);
1189:       }
1190:     } else if (xnext->type == VECPACK_DA) {
1191:       VecGetArray(x,&xarray);
1192:       VecGetArray(y,&yarray);
1193:       DAGetGlobalVector(xnext->da,&xglobal);
1194:       DAGetGlobalVector(ynext->da,&yglobal);
1195:       VecPlaceArray(xglobal,xarray+xnext->rstart);
1196:       VecPlaceArray(yglobal,yarray+ynext->rstart);
1197:       MatMultTranspose(anext->A,xglobal,yglobal);
1198:       VecRestoreArray(x,&xarray);
1199:       VecRestoreArray(y,&yarray);
1200:       VecResetArray(xglobal);
1201:       VecResetArray(yglobal);
1202:       DARestoreGlobalVector(xnext->da,&xglobal);
1203:       DARestoreGlobalVector(ynext->da,&yglobal);
1204:       anext = anext->next;
1205:     } else {
1206:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1207:     }
1208:     xnext = xnext->next;
1209:     ynext = ynext->next;
1210:   }
1211:   return(0);
1212: }

1216: PetscErrorCode MatDestroy_Shell_Pack(Mat A)
1217: {
1218:   struct MatPack     *mpack;
1219:   struct MatPackLink *anext,*oldanext;
1220:   PetscErrorCode     ierr;

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

1226:   while (anext) {
1227:     MatDestroy(anext->A);
1228:     oldanext = anext;
1229:     anext    = anext->next;
1230:     PetscFree(oldanext);
1231:   }
1232:   PetscFree(mpack);
1233:   PetscObjectChangeTypeName((PetscObject)A,0);
1234:   return(0);
1235: }

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

1242:     Collective on VecPack

1244:     Input Parameters:
1245: +    coarse - coarse grid packer
1246: -    fine - fine grid packer

1248:     Output Parameter:
1249: +    A - interpolation matrix
1250: -    v - scaling vector
1251:  
1252:     Level: advanced

1254: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
1255:          VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
1256:          VecPackGetLocalVectors(), VecPackRestoreLocalVectors()

1258: @*/
1259: PetscErrorCode  VecPackGetInterpolation(VecPack coarse,VecPack fine,Mat *A,Vec *v)
1260: {
1261:   PetscErrorCode     ierr;
1262:   PetscInt           m,n,M,N;
1263:   struct VecPackLink *nextc  = coarse->next;
1264:   struct VecPackLink *nextf = fine->next;
1265:   struct MatPackLink *nextmat,*pnextmat = 0;
1266:   struct MatPack     *mpack;
1267:   Vec                gcoarse,gfine;

1270:   /* use global vectors only for determining matrix layout */
1271:   VecPackCreateGlobalVector(coarse,&gcoarse);
1272:   VecPackCreateGlobalVector(fine,&gfine);
1273:   VecGetLocalSize(gcoarse,&n);
1274:   VecGetLocalSize(gfine,&m);
1275:   VecGetSize(gcoarse,&N);
1276:   VecGetSize(gfine,&M);
1277:   VecDestroy(gcoarse);
1278:   VecDestroy(gfine);

1280:   PetscNew(struct MatPack,&mpack);
1281:   mpack->right = coarse;
1282:   mpack->left  = fine;
1283:   MatCreate(fine->comm,A);
1284:   MatSetSizes(*A,m,n,M,N);
1285:   MatSetType(*A,MATSHELL);
1286:   MatShellSetContext(*A,mpack);
1287:   MatShellSetOperation(*A,MATOP_MULT,(void(*)(void))MatMult_Shell_Pack);
1288:   MatShellSetOperation(*A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell_Pack);
1289:   MatShellSetOperation(*A,MATOP_MULT_ADD,(void(*)(void))MatMultAdd_Shell_Pack);
1290:   MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell_Pack);

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

1296:     if (nextc->type == VECPACK_ARRAY) {
1297:       ;
1298:     } else if (nextc->type == VECPACK_DA) {
1299:       PetscNew(struct MatPackLink,&nextmat);
1300:       nextmat->next = 0;
1301:       if (pnextmat) {
1302:         pnextmat->next = nextmat;
1303:         pnextmat       = nextmat;
1304:       } else {
1305:         pnextmat    = nextmat;
1306:         mpack->next = nextmat;
1307:       }
1308:       DAGetInterpolation(nextc->da,nextf->da,&nextmat->A,PETSC_NULL);
1309:     } else {
1310:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1311:     }
1312:     nextc = nextc->next;
1313:     nextf = nextf->next;
1314:   }
1315:   return(0);
1316: }

1320: /*@C
1321:     VecPackGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for 
1322:       computing the Jacobian on a function defined using the stencils set in the DA's and coupling in the array variables

1324:     Collective on DA

1326:     Input Parameter:
1327: +   da - the distributed array
1328: -   mtype - Supported types are MATSEQAIJ, MATMPIAIJ

1330:     Output Parameters:
1331: .   J  - matrix with the correct nonzero structure
1332:         (obviously without the correct Jacobian values)

1334:     Level: advanced

1336:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you 
1337:        do not need to do it yourself. 


1340: .seealso DAGetMatrix(), VecPackCreate()

1342: @*/
1343: PetscErrorCode  VecPackGetMatrix(VecPack packer, MatType mtype,Mat *J)
1344: {
1345:   PetscErrorCode     ierr;
1346:   struct VecPackLink *next = packer->next;
1347:   Vec                gvec;
1348:   PetscInt           m,*dnz,*onz,i,j,mA;
1349:   Mat                Atmp;
1350:   PetscMPIInt        rank;
1351:   PetscScalar        zero = 0.0;

1354:   /* use global vector to determine layout needed for matrix */
1355:   VecPackCreateGlobalVector(packer,&gvec);
1356:   VecGetLocalSize(gvec,&m);
1357:   VecDestroy(gvec);

1359:   MPI_Comm_rank(packer->comm,&rank);
1360:   MatCreate(packer->comm,J);
1361:   MatSetSizes(*J,m,m,PETSC_DETERMINE,PETSC_DETERMINE);
1362:   MatSetType(*J,mtype);

1364:   MatPreallocateInitialize(packer->comm,m,m,dnz,onz);
1365:   /* loop over packed objects, handling one at at time */
1366:   while (next) {
1367:     if (next->type == VECPACK_ARRAY) {
1368:       if (rank == next->rank) {  /* zero the entire row */
1369:         for (j=packer->rstart+next->rstart; j<packer->rstart+next->rstart+next->n; j++) {
1370:           for (i=0; i<packer->N; i++) {
1371:             MatPreallocateSet(j,1,&i,dnz,onz);
1372:           }
1373:         }
1374:       }
1375:       for (j=next->grstart; j<next->grstart+next->n; j++) {
1376:         for (i=packer->rstart; i<packer->rstart+m; i++) { /* zero the entire local column */
1377:           if (j != i) { /* don't count diagonal twice */
1378:             MatPreallocateSet(i,1,&j,dnz,onz);
1379:           }
1380:         }
1381:       }
1382:     } else if (next->type == VECPACK_DA) {
1383:       PetscInt       nc,rstart,*ccols,maxnc;
1384:       const PetscInt *cols,*rstarts;
1385:       PetscMPIInt    proc;

1387:       DAGetMatrix(next->da,mtype,&Atmp);
1388:       MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);
1389:       MatGetOwnershipRanges(Atmp,&rstarts);
1390:       MatGetLocalSize(Atmp,&mA,PETSC_NULL);

1392:       maxnc = 0;
1393:       for (i=0; i<mA; i++) {
1394:         MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1395:         MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1396:         maxnc = PetscMax(nc,maxnc);
1397:       }
1398:       PetscMalloc(maxnc*sizeof(PetscInt),&ccols);
1399:       for (i=0; i<mA; i++) {
1400:         MatGetRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);
1401:         /* remap the columns taking into how much they are shifted on each process */
1402:         for (j=0; j<nc; j++) {
1403:           proc = 0;
1404:           while (cols[j] >= rstarts[proc+1]) proc++;
1405:           ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
1406:         }
1407:         MatPreallocateSet(packer->rstart+next->rstart+i,nc,ccols,dnz,onz);
1408:         MatRestoreRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);
1409:       }
1410:       PetscFree(ccols);
1411:     } else {
1412:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1413:     }
1414:     next = next->next;
1415:   }
1416:   MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
1417:   MatSeqAIJSetPreallocation(*J,0,dnz);
1418:   MatPreallocateFinalize(dnz,onz);

1420:   next = packer->next;
1421:   while (next) {
1422:     if (next->type == VECPACK_ARRAY) {
1423:       if (rank == next->rank) {
1424:         for (j=packer->rstart+next->rstart; j<packer->rstart+next->rstart+next->n; j++) {
1425:           for (i=0; i<packer->N; i++) {
1426:             MatSetValues(*J,1,&j,1,&i,&zero,INSERT_VALUES);
1427:           }
1428:         }
1429:       }
1430:       for (j=next->grstart; j<next->grstart+next->n; j++) {
1431:         for (i=packer->rstart; i<packer->rstart+m; i++) {
1432:           MatSetValues(*J,1,&i,1,&j,&zero,INSERT_VALUES);
1433:         }
1434:       }
1435:     } else if (next->type == VECPACK_DA) {
1436:       PetscInt          nc,rstart,row,maxnc,*ccols;
1437:       const PetscInt    *cols,*rstarts;
1438:       const PetscScalar *values;
1439:       PetscMPIInt       proc;

1441:       DAGetMatrix(next->da,mtype,&Atmp);
1442:       MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);
1443:       MatGetOwnershipRanges(Atmp,&rstarts);
1444:       MatGetLocalSize(Atmp,&mA,PETSC_NULL);
1445:       maxnc = 0;
1446:       for (i=0; i<mA; i++) {
1447:         MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1448:         MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1449:         maxnc = PetscMax(nc,maxnc);
1450:       }
1451:       PetscMalloc(maxnc*sizeof(PetscInt),&ccols);
1452:       for (i=0; i<mA; i++) {
1453:         MatGetRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);
1454:         for (j=0; j<nc; j++) {
1455:           proc = 0;
1456:           while (cols[j] >= rstarts[proc+1]) proc++;
1457:           ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
1458:         }
1459:         row  = packer->rstart+next->rstart+i;
1460:         MatSetValues(*J,1,&row,nc,ccols,values,INSERT_VALUES);
1461:         MatRestoreRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);
1462:       }
1463:       PetscFree(ccols);
1464:     } else {
1465:       SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1466:     }
1467:     next = next->next;
1468:   }
1469:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1470:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1471:   return(0);
1472: }

1476: /*@
1477:     VecPackGetColoring - Gets the coloring required for computing the Jacobian via
1478:     finite differences on a function defined using a VecPack "grid"

1480:     Collective on DA

1482:     Input Parameter:
1483: +   vecpack - the VecPack object
1484: -   ctype - IS_COLORING_GLOBAL or IS_COLORING_GHOSTED

1486:     Output Parameters:
1487: .   coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)

1489:     Level: advanced

1491:     Notes: This currentlu uses one color per column so is very slow.

1493:     Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used 
1494:    for efficient (parallel or thread based) triangular solves etc is NOT yet 
1495:    available. 


1498: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring, DAGetColoring()

1500: @*/
1501: PetscErrorCode  VecPackGetColoring(VecPack vecpack,ISColoringType ctype,ISColoring *coloring)
1502: {
1503:   PetscErrorCode  ierr;
1504:   PetscInt        n,i;
1505:   ISColoringValue *colors;

1508:   if (ctype == IS_COLORING_GHOSTED) {
1509:     SETERRQ(PETSC_ERR_SUP,"Lazy Barry");
1510:   } else if (ctype == IS_COLORING_GLOBAL) {
1511:     n = vecpack->n;
1512:   } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");

1514:   PetscMalloc(n*sizeof(ISColoringValue),&colors); /* freed in ISColoringDestroy() */
1515:   for (i=0; i<n; i++) {
1516:     colors[i] = (ISColoringValue)(vecpack->rstart + i);
1517:   }
1518:   ISColoringCreate(vecpack->comm,vecpack->N,n,colors,coloring);
1519:   return(0);
1520: }