Actual source code: pack.c
1: #define PETSCDM_DLL
2:
3: #include petscda.h
4: #include src/dm/dmimpl.h
5: #include petscmat.h
7: typedef struct _DMCompositeOps *DMCompositeOps;
8: struct _DMCompositeOps {
9: DMOPS(DMComposite)
10: };
12: /*
13: rstart is where an array/subvector starts in the global parallel vector, so arrays
14: rstarts are meaningless (and set to the previous one) except on processor 0
15: */
17: typedef enum {DMCOMPOSITE_ARRAY, DMCOMPOSITE_DA, DMCOMPOSITE_VECSCATTER} DMCompositeLinkType;
19: struct DMCompositeLink {
20: DMCompositeLinkType type;
21: struct DMCompositeLink *next;
22: PetscInt n,rstart; /* rstart is relative to this processor */
24: /* only used for DMCOMPOSITE_DA */
25: PetscInt *grstarts; /* global row for first unknown of this DA on each process */
26: DA da;
28: /* only used for DMCOMPOSITE_ARRAY */
29: PetscInt grstart; /* global row for first array unknown */
30: PetscMPIInt rank; /* process where array unknowns live */
31: };
33: struct _p_DMComposite {
34: PETSCHEADER(struct _DMCompositeOps);
35: PetscInt n,N,rstart; /* rstart is relative to all processors, n unknowns owned by this process, N is total unknowns */
36: PetscInt nghost; /* number of all local entries include DA ghost points and any shared redundant arrays */
37: PetscInt nDA,nredundant; /* how many DA's and seperate redundant arrays used to build DMComposite */
38: PetscTruth setup; /* after this is set, cannot add new links to the DMComposite */
39: struct DMCompositeLink *next;
40: };
45: /*@C
46: DMCompositeCreate - Creates a vector packer, used to generate "composite"
47: vectors made up of several subvectors.
49: Collective on MPI_Comm
51: Input Parameter:
52: . comm - the processors that will share the global vector
54: Output Parameters:
55: . packer - the packer object
57: Level: advanced
59: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeScatter(),
60: DMCompositeGather(), DMCompositeCreateGlobalVector(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess()
61: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
63: @*/
64: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DMComposite *packer)
65: {
67: DMComposite p;
71: *packer = PETSC_NULL;
72: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
73: DMInitializePackage(PETSC_NULL);
74: #endif
76: PetscHeaderCreate(p,_p_DMComposite,struct _DMCompositeOps,DA_COOKIE,0,"DMComposite",comm,DMCompositeDestroy,0);
77: p->n = 0;
78: p->next = PETSC_NULL;
79: p->comm = comm;
80: p->nredundant = 0;
81: p->nDA = 0;
83: p->ops->createglobalvector = DMCompositeCreateGlobalVector;
84: p->ops->refine = DMCompositeRefine;
85: p->ops->getinterpolation = DMCompositeGetInterpolation;
86: p->ops->getmatrix = DMCompositeGetMatrix;
87: p->ops->getcoloring = DMCompositeGetColoring;
89: *packer = p;
90: return(0);
91: }
95: /*@C
96: DMCompositeDestroy - Destroys a vector packer.
98: Collective on DMComposite
100: Input Parameter:
101: . packer - the packer object
103: Level: advanced
105: .seealso DMCompositeCreate(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeScatter(),DMCompositeGetEntries()
106: DMCompositeGather(), DMCompositeCreateGlobalVector(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess()
108: @*/
109: PetscErrorCode DMCompositeDestroy(DMComposite packer)
110: {
111: PetscErrorCode ierr;
112: struct DMCompositeLink *next, *prev;
116: next = packer->next;
117: if (--packer->refct > 0) return(0);
118: while (next) {
119: prev = next;
120: next = next->next;
121: if (prev->type == DMCOMPOSITE_DA) {
122: DADestroy(prev->da);
123: }
124: if (prev->grstarts) {
125: PetscFree(prev->grstarts);
126: }
127: PetscFree(prev);
128: }
129: PetscHeaderDestroy(packer);
130: return(0);
131: }
133: /* --------------------------------------------------------------------------------------*/
136: PetscErrorCode DMCompositeSetUp(DMComposite packer)
137: {
138: PetscErrorCode ierr;
139: PetscInt nprev = 0;
140: PetscMPIInt rank,size;
141: struct DMCompositeLink *next = packer->next;
142: PetscMap map;
145: if (packer->setup) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
146: PetscMapInitialize(packer->comm,&map);
147: PetscMapSetLocalSize(&map,packer->n);
148: PetscMapSetSize(&map,PETSC_DETERMINE);
149: PetscMapSetBlockSize(&map,1);
150: PetscMapSetUp(&map);
151: PetscMapGetSize(&map,&packer->N);
152: PetscMapGetLocalRange(&map,&packer->rstart,PETSC_NULL);
153: PetscFree(map.range);
154:
155: /* now set the rstart for each linked array/vector */
156: MPI_Comm_rank(packer->comm,&rank);
157: MPI_Comm_size(packer->comm,&size);
158: while (next) {
159: next->rstart = nprev;
160: if ((rank == next->rank) || next->type != DMCOMPOSITE_ARRAY) nprev += next->n;
161: next->grstart = packer->rstart + next->rstart;
162: if (next->type == DMCOMPOSITE_ARRAY) {
163: MPI_Bcast(&next->grstart,1,MPIU_INT,next->rank,packer->comm);
164: } else {
165: PetscMalloc(size*sizeof(PetscInt),&next->grstarts);
166: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,packer->comm);
167: }
168: next = next->next;
169: }
170: packer->setup = PETSC_TRUE;
171: return(0);
172: }
177: PetscErrorCode DMCompositeGetAccess_Array(DMComposite packer,struct DMCompositeLink *mine,Vec vec,PetscScalar **array)
178: {
180: PetscScalar *varray;
181: PetscMPIInt rank;
184: MPI_Comm_rank(packer->comm,&rank);
185: if (array) {
186: if (rank == mine->rank) {
187: VecGetArray(vec,&varray);
188: *array = varray + mine->rstart;
189: VecRestoreArray(vec,&varray);
190: } else {
191: *array = 0;
192: }
193: }
194: return(0);
195: }
199: PetscErrorCode DMCompositeGetAccess_DA(DMComposite packer,struct DMCompositeLink *mine,Vec vec,Vec *global)
200: {
202: PetscScalar *array;
205: if (global) {
206: DAGetGlobalVector(mine->da,global);
207: VecGetArray(vec,&array);
208: VecPlaceArray(*global,array+mine->rstart);
209: VecRestoreArray(vec,&array);
210: }
211: return(0);
212: }
216: PetscErrorCode DMCompositeRestoreAccess_Array(DMComposite packer,struct DMCompositeLink *mine,Vec vec,PetscScalar **array)
217: {
219: return(0);
220: }
224: PetscErrorCode DMCompositeRestoreAccess_DA(DMComposite packer,struct DMCompositeLink *mine,Vec vec,Vec *global)
225: {
229: if (global) {
230: VecResetArray(*global);
231: DARestoreGlobalVector(mine->da,global);
232: }
233: return(0);
234: }
238: PetscErrorCode DMCompositeScatter_Array(DMComposite packer,struct DMCompositeLink *mine,Vec vec,PetscScalar *array)
239: {
241: PetscScalar *varray;
242: PetscMPIInt rank;
245: MPI_Comm_rank(packer->comm,&rank);
246: if (rank == mine->rank) {
247: VecGetArray(vec,&varray);
248: PetscMemcpy(array,varray+mine->rstart,mine->n*sizeof(PetscScalar));
249: VecRestoreArray(vec,&varray);
250: }
251: MPI_Bcast(array,mine->n,MPIU_SCALAR,mine->rank,packer->comm);
252: return(0);
253: }
257: PetscErrorCode DMCompositeScatter_DA(DMComposite packer,struct DMCompositeLink *mine,Vec vec,Vec local)
258: {
260: PetscScalar *array;
261: Vec global;
264: DAGetGlobalVector(mine->da,&global);
265: VecGetArray(vec,&array);
266: VecPlaceArray(global,array+mine->rstart);
267: DAGlobalToLocalBegin(mine->da,global,INSERT_VALUES,local);
268: DAGlobalToLocalEnd(mine->da,global,INSERT_VALUES,local);
269: VecRestoreArray(vec,&array);
270: VecResetArray(global);
271: DARestoreGlobalVector(mine->da,&global);
272: return(0);
273: }
277: PetscErrorCode DMCompositeGather_Array(DMComposite packer,struct DMCompositeLink *mine,Vec vec,PetscScalar *array)
278: {
280: PetscScalar *varray;
281: PetscMPIInt rank;
284: MPI_Comm_rank(packer->comm,&rank);
285: if (rank == mine->rank) {
286: VecGetArray(vec,&varray);
287: if (varray+mine->rstart == array) SETERRQ(PETSC_ERR_ARG_WRONG,"You need not DMCompositeGather() into objects obtained via DMCompositeGetAccess()");
288: PetscMemcpy(varray+mine->rstart,array,mine->n*sizeof(PetscScalar));
289: VecRestoreArray(vec,&varray);
290: }
291: return(0);
292: }
296: PetscErrorCode DMCompositeGather_DA(DMComposite packer,struct DMCompositeLink *mine,Vec vec,Vec local)
297: {
299: PetscScalar *array;
300: Vec global;
303: DAGetGlobalVector(mine->da,&global);
304: VecGetArray(vec,&array);
305: VecPlaceArray(global,array+mine->rstart);
306: DALocalToGlobal(mine->da,local,INSERT_VALUES,global);
307: VecRestoreArray(vec,&array);
308: VecResetArray(global);
309: DARestoreGlobalVector(mine->da,&global);
310: return(0);
311: }
313: /* ----------------------------------------------------------------------------------*/
315: #include <stdarg.h>
319: /*@C
320: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
321: representation.
323: Collective on DMComposite
325: Input Parameter:
326: + packer - the packer object
327: . gvec - the global vector
328: - ... - the individual sequential or parallel objects (arrays or vectors)
330: Notes: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
331:
332: Level: advanced
334: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
335: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeScatter(),
336: DMCompositeRestoreAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),
337: DMCompositeGetEntries()
339: @*/
340: PetscErrorCode DMCompositeGetAccess(DMComposite packer,Vec gvec,...)
341: {
342: va_list Argp;
343: PetscErrorCode ierr;
344: struct DMCompositeLink *next;
349: next = packer->next;
350: if (!packer->setup) {
351: DMCompositeSetUp(packer);
352: }
354: /* loop over packed objects, handling one at at time */
355: va_start(Argp,gvec);
356: while (next) {
357: if (next->type == DMCOMPOSITE_ARRAY) {
358: PetscScalar **array;
359: array = va_arg(Argp, PetscScalar**);
360: DMCompositeGetAccess_Array(packer,next,gvec,array);
361: } else if (next->type == DMCOMPOSITE_DA) {
362: Vec *vec;
363: vec = va_arg(Argp, Vec*);
364: DMCompositeGetAccess_DA(packer,next,gvec,vec);
365: } else {
366: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
367: }
368: next = next->next;
369: }
370: va_end(Argp);
371: return(0);
372: }
376: /*@C
377: DMCompositeRestoreAccess - Returns the vectors obtained with DACompositeGetAccess()
378: representation.
380: Collective on DMComposite
382: Input Parameter:
383: + packer - the packer object
384: . gvec - the global vector
385: - ... - the individual sequential or parallel objects (arrays or vectors)
386:
387: Level: advanced
389: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
390: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeScatter(),
391: DMCompositeRestoreAccess(), DACompositeGetAccess()
393: @*/
394: PetscErrorCode DMCompositeRestoreAccess(DMComposite packer,Vec gvec,...)
395: {
396: va_list Argp;
397: PetscErrorCode ierr;
398: struct DMCompositeLink *next;
403: next = packer->next;
404: if (!packer->setup) {
405: DMCompositeSetUp(packer);
406: }
408: /* loop over packed objects, handling one at at time */
409: va_start(Argp,gvec);
410: while (next) {
411: if (next->type == DMCOMPOSITE_ARRAY) {
412: PetscScalar **array;
413: array = va_arg(Argp, PetscScalar**);
414: DMCompositeRestoreAccess_Array(packer,next,gvec,array);
415: } else if (next->type == DMCOMPOSITE_DA) {
416: Vec *vec;
417: vec = va_arg(Argp, Vec*);
418: DMCompositeRestoreAccess_DA(packer,next,gvec,vec);
419: } else {
420: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
421: }
422: next = next->next;
423: }
424: va_end(Argp);
425: return(0);
426: }
430: /*@C
431: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
433: Collective on DMComposite
435: Input Parameter:
436: + packer - the packer object
437: . gvec - the global vector
438: - ... - the individual sequential objects (arrays or vectors)
439:
440: Level: advanced
442: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
443: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
444: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
446: @*/
447: PetscErrorCode DMCompositeScatter(DMComposite packer,Vec gvec,...)
448: {
449: va_list Argp;
450: PetscErrorCode ierr;
451: struct DMCompositeLink *next;
452: PetscInt cnt = 3;
457: next = packer->next;
458: if (!packer->setup) {
459: DMCompositeSetUp(packer);
460: }
462: /* loop over packed objects, handling one at at time */
463: va_start(Argp,gvec);
464: while (next) {
465: if (next->type == DMCOMPOSITE_ARRAY) {
466: PetscScalar *array;
467: array = va_arg(Argp, PetscScalar*);
468: DMCompositeScatter_Array(packer,next,gvec,array);
469: } else if (next->type == DMCOMPOSITE_DA) {
470: Vec vec;
471: vec = va_arg(Argp, Vec);
473: DMCompositeScatter_DA(packer,next,gvec,vec);
474: } else {
475: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
476: }
477: cnt++;
478: next = next->next;
479: }
480: va_end(Argp);
481: return(0);
482: }
486: /*@C
487: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
489: Collective on DMComposite
491: Input Parameter:
492: + packer - the packer object
493: . gvec - the global vector
494: - ... - the individual sequential objects (arrays or vectors)
495:
496: Level: advanced
498: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
499: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
500: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
502: @*/
503: PetscErrorCode DMCompositeGather(DMComposite packer,Vec gvec,...)
504: {
505: va_list Argp;
506: PetscErrorCode ierr;
507: struct DMCompositeLink *next;
512: next = packer->next;
513: if (!packer->setup) {
514: DMCompositeSetUp(packer);
515: }
517: /* loop over packed objects, handling one at at time */
518: va_start(Argp,gvec);
519: while (next) {
520: if (next->type == DMCOMPOSITE_ARRAY) {
521: PetscScalar *array;
522: array = va_arg(Argp, PetscScalar*);
523: DMCompositeGather_Array(packer,next,gvec,array);
524: } else if (next->type == DMCOMPOSITE_DA) {
525: Vec vec;
526: vec = va_arg(Argp, Vec);
528: DMCompositeGather_DA(packer,next,gvec,vec);
529: } else {
530: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
531: }
532: next = next->next;
533: }
534: va_end(Argp);
535: return(0);
536: }
540: /*@C
541: DMCompositeAddArray - adds an "redundant" array to a DMComposite. The array values will
542: be stored in part of the array on process orank.
544: Collective on DMComposite
546: Input Parameter:
547: + packer - the packer object
548: . orank - the process on which the array entries officially live, this number must be
549: the same on all processes.
550: - n - the length of the array
551:
552: Level: advanced
554: .seealso DMCompositeDestroy(), DMCompositeGather(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
555: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
556: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
558: @*/
559: PetscErrorCode DMCompositeAddArray(DMComposite packer,PetscMPIInt orank,PetscInt n)
560: {
561: struct DMCompositeLink *mine,*next;
562: PetscErrorCode ierr;
563: PetscMPIInt rank;
567: next = packer->next;
568: if (packer->setup) {
569: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add an array once you have used the DMComposite");
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 DMCompositeLink,&mine);
582: mine->n = n;
583: mine->rank = orank;
584: mine->da = PETSC_NULL;
585: mine->type = DMCOMPOSITE_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: DMCompositeAddDA - adds a DA vector to a DMComposite
605: Collective on DMComposite
607: Input Parameter:
608: + packer - the packer object
609: - da - the DA object
610:
611: Level: advanced
613: .seealso DMCompositeDestroy(), DMCompositeGather(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
614: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
615: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
617: @*/
618: PetscErrorCode DMCompositeAddDA(DMComposite packer,DA da)
619: {
620: PetscErrorCode ierr;
621: PetscInt n;
622: struct DMCompositeLink *mine,*next;
623: Vec global;
628: next = packer->next;
629: if (packer->setup) {
630: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DA once you have used the DMComposite");
631: }
633: /* create new link */
634: PetscNew(struct DMCompositeLink,&mine);
635: PetscObjectReference((PetscObject)da);
636: DAGetGlobalVector(da,&global);
637: VecGetLocalSize(global,&n);
638: DARestoreGlobalVector(da,&global);
639: mine->n = n;
640: mine->da = da;
641: mine->type = DMCOMPOSITE_DA;
642: mine->next = PETSC_NULL;
643: packer->n += n;
645: /* add to end of list */
646: if (!next) {
647: packer->next = mine;
648: } else {
649: while (next->next) next = next->next;
650: next->next = mine;
651: }
652: packer->nDA++;
653: return(0);
654: }
660: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
661: {
662: DMComposite packer;
663: PetscErrorCode ierr;
664: struct DMCompositeLink *next;
665: PetscTruth isdraw;
668: PetscObjectQuery((PetscObject)gvec,"DMComposite",(PetscObject*)&packer);
669: if (!packer) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
670: next = packer->next;
672: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
673: if (!isdraw) {
674: /* do I really want to call this? */
675: VecView_MPI(gvec,viewer);
676: } else {
677: PetscInt cnt = 0;
679: /* loop over packed objects, handling one at at time */
680: while (next) {
681: if (next->type == DMCOMPOSITE_ARRAY) {
682: PetscScalar *array;
683: DMCompositeGetAccess_Array(packer,next,gvec,&array);
685: /*skip it for now */
686: } else if (next->type == DMCOMPOSITE_DA) {
687: Vec vec;
688: PetscInt bs;
690: DMCompositeGetAccess_DA(packer,next,gvec,&vec);
691: VecView(vec,viewer);
692: VecGetBlockSize(vec,&bs);
693: DMCompositeRestoreAccess_DA(packer,next,gvec,&vec);
694: PetscViewerDrawBaseAdd(viewer,bs);
695: cnt += bs;
696: } else {
697: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
698: }
699: next = next->next;
700: }
701: PetscViewerDrawBaseAdd(viewer,-cnt);
702: }
703: return(0);
704: }
709: /*@C
710: DMCompositeCreateGlobalVector - Creates a vector of the correct size to be gathered into
711: by the packer.
713: Collective on DMComposite
715: Input Parameter:
716: . packer - the packer object
718: Output Parameters:
719: . gvec - the global vector
721: Level: advanced
723: Notes: Once this has been created you cannot add additional arrays or vectors to be packed.
725: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeScatter(),
726: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
727: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
729: @*/
730: PetscErrorCode DMCompositeCreateGlobalVector(DMComposite packer,Vec *gvec)
731: {
732: PetscErrorCode ierr;
736: if (!packer->setup) {
737: DMCompositeSetUp(packer);
738: }
739: VecCreateMPI(packer->comm,packer->n,packer->N,gvec);
740: PetscObjectCompose((PetscObject)*gvec,"DMComposite",(PetscObject)packer);
741: VecSetOperation(*gvec,VECOP_VIEW,(void(*)(void))VecView_DMComposite);
742: return(0);
743: }
747: /*@C
748: DMCompositeGetGlobalIndices - Gets the global indices for all the entries in the packed
749: vectors.
751: Collective on DMComposite
753: Input Parameter:
754: . packer - the packer object
756: Output Parameters:
757: . idx - the individual indices for each packed vector/array. Note that this includes
758: all the ghost points that individual ghosted DA's may have.
759:
760: Level: advanced
762: Notes:
763: The idx parameters should be freed by the calling routine with PetscFree()
765: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
766: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
767: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
769: @*/
770: PetscErrorCode DMCompositeGetGlobalIndices(DMComposite packer,...)
771: {
772: va_list Argp;
773: PetscErrorCode ierr;
774: PetscInt i,**idx,n;
775: struct DMCompositeLink *next;
776: Vec global,dglobal;
777: PF pf;
778: PetscScalar *array;
779: PetscMPIInt rank;
783: next = packer->next;
784: DMCompositeCreateGlobalVector(packer,&global);
785: MPI_Comm_rank(packer->comm,&rank);
787: /* put 0 to N-1 into the global vector */
788: PFCreate(PETSC_COMM_WORLD,1,1,&pf);
789: PFSetType(pf,PFIDENTITY,PETSC_NULL);
790: PFApplyVec(pf,PETSC_NULL,global);
791: PFDestroy(pf);
793: /* loop over packed objects, handling one at at time */
794: va_start(Argp,packer);
795: while (next) {
796: idx = va_arg(Argp, PetscInt**);
798: if (next->type == DMCOMPOSITE_ARRAY) {
799:
800: PetscMalloc(next->n*sizeof(PetscInt),idx);
801: if (rank == next->rank) {
802: VecGetArray(global,&array);
803: array += next->rstart;
804: for (i=0; i<next->n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
805: array -= next->rstart;
806: VecRestoreArray(global,&array);
807: }
808: MPI_Bcast(*idx,next->n,MPIU_INT,next->rank,packer->comm);
810: } else if (next->type == DMCOMPOSITE_DA) {
811: Vec local;
813: DACreateLocalVector(next->da,&local);
814: VecGetArray(global,&array);
815: array += next->rstart;
816: DAGetGlobalVector(next->da,&dglobal);
817: VecPlaceArray(dglobal,array);
818: DAGlobalToLocalBegin(next->da,dglobal,INSERT_VALUES,local);
819: DAGlobalToLocalEnd(next->da,dglobal,INSERT_VALUES,local);
820: array -= next->rstart;
821: VecRestoreArray(global,&array);
822: VecResetArray(dglobal);
823: DARestoreGlobalVector(next->da,&dglobal);
825: VecGetArray(local,&array);
826: VecGetSize(local,&n);
827: PetscMalloc(n*sizeof(PetscInt),idx);
828: for (i=0; i<n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
829: VecRestoreArray(local,&array);
830: VecDestroy(local);
832: } else {
833: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
834: }
835: next = next->next;
836: }
837: va_end(Argp);
838: VecDestroy(global);
839: return(0);
840: }
842: /* -------------------------------------------------------------------------------------*/
845: PetscErrorCode DMCompositeGetLocalVectors_Array(DMComposite packer,struct DMCompositeLink *mine,PetscScalar **array)
846: {
849: if (array) {
850: PetscMalloc(mine->n*sizeof(PetscScalar),array);
851: }
852: return(0);
853: }
857: PetscErrorCode DMCompositeGetLocalVectors_DA(DMComposite packer,struct DMCompositeLink *mine,Vec *local)
858: {
861: if (local) {
862: DAGetLocalVector(mine->da,local);
863: }
864: return(0);
865: }
869: PetscErrorCode DMCompositeRestoreLocalVectors_Array(DMComposite packer,struct DMCompositeLink *mine,PetscScalar **array)
870: {
873: if (array) {
874: PetscFree(*array);
875: }
876: return(0);
877: }
881: PetscErrorCode DMCompositeRestoreLocalVectors_DA(DMComposite packer,struct DMCompositeLink *mine,Vec *local)
882: {
885: if (local) {
886: DARestoreLocalVector(mine->da,local);
887: }
888: return(0);
889: }
893: /*@C
894: DMCompositeGetLocalVectors - Gets local vectors and arrays for each part of a DMComposite.'
895: Use VecPakcRestoreLocalVectors() to return them.
897: Collective on DMComposite
899: Input Parameter:
900: . packer - the packer object
901:
902: Output Parameter:
903: . ... - the individual sequential objects (arrays or vectors)
904:
905: Level: advanced
907: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
908: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
909: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
911: @*/
912: PetscErrorCode DMCompositeGetLocalVectors(DMComposite packer,...)
913: {
914: va_list Argp;
915: PetscErrorCode ierr;
916: struct DMCompositeLink *next;
920: next = packer->next;
921: /* loop over packed objects, handling one at at time */
922: va_start(Argp,packer);
923: while (next) {
924: if (next->type == DMCOMPOSITE_ARRAY) {
925: PetscScalar **array;
926: array = va_arg(Argp, PetscScalar**);
927: DMCompositeGetLocalVectors_Array(packer,next,array);
928: } else if (next->type == DMCOMPOSITE_DA) {
929: Vec *vec;
930: vec = va_arg(Argp, Vec*);
931: DMCompositeGetLocalVectors_DA(packer,next,vec);
932: } else {
933: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
934: }
935: next = next->next;
936: }
937: va_end(Argp);
938: return(0);
939: }
943: /*@C
944: DMCompositeRestoreLocalVectors - Restores local vectors and arrays for each part of a DMComposite.'
945: Use VecPakcRestoreLocalVectors() to return them.
947: Collective on DMComposite
949: Input Parameter:
950: . packer - the packer object
951:
952: Output Parameter:
953: . ... - the individual sequential objects (arrays or vectors)
954:
955: Level: advanced
957: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
958: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
959: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
961: @*/
962: PetscErrorCode DMCompositeRestoreLocalVectors(DMComposite packer,...)
963: {
964: va_list Argp;
965: PetscErrorCode ierr;
966: struct DMCompositeLink *next;
970: next = packer->next;
971: /* loop over packed objects, handling one at at time */
972: va_start(Argp,packer);
973: while (next) {
974: if (next->type == DMCOMPOSITE_ARRAY) {
975: PetscScalar **array;
976: array = va_arg(Argp, PetscScalar**);
977: DMCompositeRestoreLocalVectors_Array(packer,next,array);
978: } else if (next->type == DMCOMPOSITE_DA) {
979: Vec *vec;
980: vec = va_arg(Argp, Vec*);
981: DMCompositeRestoreLocalVectors_DA(packer,next,vec);
982: } else {
983: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
984: }
985: next = next->next;
986: }
987: va_end(Argp);
988: return(0);
989: }
991: /* -------------------------------------------------------------------------------------*/
994: PetscErrorCode DMCompositeGetEntries_Array(DMComposite packer,struct DMCompositeLink *mine,PetscInt *n)
995: {
997: if (n) *n = mine->n;
998: return(0);
999: }
1003: PetscErrorCode DMCompositeGetEntries_DA(DMComposite packer,struct DMCompositeLink *mine,DA *da)
1004: {
1006: if (da) *da = mine->da;
1007: return(0);
1008: }
1012: /*@C
1013: DMCompositeGetEntries - Gets the DA, redundant size, etc for each entry in a DMComposite.
1015: Collective on DMComposite
1017: Input Parameter:
1018: . packer - the packer object
1019:
1020: Output Parameter:
1021: . ... - the individual entries, DAs or integer sizes)
1022:
1023: Level: advanced
1025: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
1026: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
1027: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1028: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1030: @*/
1031: PetscErrorCode DMCompositeGetEntries(DMComposite packer,...)
1032: {
1033: va_list Argp;
1034: PetscErrorCode ierr;
1035: struct DMCompositeLink *next;
1039: next = packer->next;
1040: /* loop over packed objects, handling one at at time */
1041: va_start(Argp,packer);
1042: while (next) {
1043: if (next->type == DMCOMPOSITE_ARRAY) {
1044: PetscInt *n;
1045: n = va_arg(Argp, PetscInt*);
1046: DMCompositeGetEntries_Array(packer,next,n);
1047: } else if (next->type == DMCOMPOSITE_DA) {
1048: DA *da;
1049: da = va_arg(Argp, DA*);
1050: DMCompositeGetEntries_DA(packer,next,da);
1051: } else {
1052: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1053: }
1054: next = next->next;
1055: }
1056: va_end(Argp);
1057: return(0);
1058: }
1062: /*@C
1063: DMCompositeRefine - Refines a DMComposite by refining all of its DAs
1065: Collective on DMComposite
1067: Input Parameters:
1068: + packer - the packer object
1069: - comm - communicator to contain the new DM object, usually PETSC_NULL
1071: Output Parameter:
1072: . fine - new packer
1073:
1074: Level: advanced
1076: .seealso DMCompositeDestroy(), DMCompositeAddArray(), DMCompositeAddDA(), DMCompositeCreateGlobalVector(),
1077: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetGlobalIndices(), DMCompositeGetAccess(),
1078: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeScatter(),
1079: DMCompositeGetEntries()
1081: @*/
1082: PetscErrorCode DMCompositeRefine(DMComposite packer,MPI_Comm comm,DMComposite *fine)
1083: {
1084: PetscErrorCode ierr;
1085: struct DMCompositeLink *next;
1086: DA da;
1090: next = packer->next;
1091: DMCompositeCreate(comm,fine);
1093: /* loop over packed objects, handling one at at time */
1094: while (next) {
1095: if (next->type == DMCOMPOSITE_ARRAY) {
1096: DMCompositeAddArray(*fine,next->rank,next->n);
1097: } else if (next->type == DMCOMPOSITE_DA) {
1098: DARefine(next->da,comm,&da);
1099: DMCompositeAddDA(*fine,da);
1100: PetscObjectDereference((PetscObject)da);
1101: } else {
1102: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1103: }
1104: next = next->next;
1105: }
1106: return(0);
1107: }
1109: #include petscmat.h
1111: struct MatPackLink {
1112: Mat A;
1113: struct MatPackLink *next;
1114: };
1116: struct MatPack {
1117: DMComposite right,left;
1118: struct MatPackLink *next;
1119: };
1123: PetscErrorCode MatMultBoth_Shell_Pack(Mat A,Vec x,Vec y,PetscTruth add)
1124: {
1125: struct MatPack *mpack;
1126: struct DMCompositeLink *xnext,*ynext;
1127: struct MatPackLink *anext;
1128: PetscScalar *xarray,*yarray;
1129: PetscErrorCode ierr;
1130: PetscInt i;
1131: Vec xglobal,yglobal;
1132: PetscMPIInt rank;
1135: MatShellGetContext(A,(void**)&mpack);
1136: MPI_Comm_rank(mpack->right->comm,&rank);
1137: xnext = mpack->right->next;
1138: ynext = mpack->left->next;
1139: anext = mpack->next;
1141: while (xnext) {
1142: if (xnext->type == DMCOMPOSITE_ARRAY) {
1143: if (rank == xnext->rank) {
1144: VecGetArray(x,&xarray);
1145: VecGetArray(y,&yarray);
1146: if (add) {
1147: for (i=0; i<xnext->n; i++) {
1148: yarray[ynext->rstart+i] += xarray[xnext->rstart+i];
1149: }
1150: } else {
1151: PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1152: }
1153: VecRestoreArray(x,&xarray);
1154: VecRestoreArray(y,&yarray);
1155: }
1156: } else if (xnext->type == DMCOMPOSITE_DA) {
1157: VecGetArray(x,&xarray);
1158: VecGetArray(y,&yarray);
1159: DAGetGlobalVector(xnext->da,&xglobal);
1160: DAGetGlobalVector(ynext->da,&yglobal);
1161: VecPlaceArray(xglobal,xarray+xnext->rstart);
1162: VecPlaceArray(yglobal,yarray+ynext->rstart);
1163: if (add) {
1164: MatMultAdd(anext->A,xglobal,yglobal,yglobal);
1165: } else {
1166: MatMult(anext->A,xglobal,yglobal);
1167: }
1168: VecRestoreArray(x,&xarray);
1169: VecRestoreArray(y,&yarray);
1170: VecResetArray(xglobal);
1171: VecResetArray(yglobal);
1172: DARestoreGlobalVector(xnext->da,&xglobal);
1173: DARestoreGlobalVector(ynext->da,&yglobal);
1174: anext = anext->next;
1175: } else {
1176: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1177: }
1178: xnext = xnext->next;
1179: ynext = ynext->next;
1180: }
1181: return(0);
1182: }
1186: PetscErrorCode MatMultAdd_Shell_Pack(Mat A,Vec x,Vec y,Vec z)
1187: {
1190: if (z != y) SETERRQ(PETSC_ERR_SUP,"Handles y == z only");
1191: MatMultBoth_Shell_Pack(A,x,y,PETSC_TRUE);
1192: return(0);
1193: }
1197: PetscErrorCode MatMult_Shell_Pack(Mat A,Vec x,Vec y)
1198: {
1201: MatMultBoth_Shell_Pack(A,x,y,PETSC_FALSE);
1202: return(0);
1203: }
1207: PetscErrorCode MatMultTranspose_Shell_Pack(Mat A,Vec x,Vec y)
1208: {
1209: struct MatPack *mpack;
1210: struct DMCompositeLink *xnext,*ynext;
1211: struct MatPackLink *anext;
1212: PetscScalar *xarray,*yarray;
1213: PetscErrorCode ierr;
1214: Vec xglobal,yglobal;
1215: PetscMPIInt rank;
1218: MatShellGetContext(A,(void**)&mpack);
1219: MPI_Comm_rank(mpack->right->comm,&rank);
1220: xnext = mpack->left->next;
1221: ynext = mpack->right->next;
1222: anext = mpack->next;
1224: while (xnext) {
1225: if (xnext->type == DMCOMPOSITE_ARRAY) {
1226: if (rank == ynext->rank) {
1227: VecGetArray(x,&xarray);
1228: VecGetArray(y,&yarray);
1229: PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1230: VecRestoreArray(x,&xarray);
1231: