Actual source code: pack.c
1: /*$Id: pack.c,v 1.19 2001/08/07 03:04:45 balay Exp $*/
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: int 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: int (*view)(VecPack,PetscViewer);
23: int (*createglobalvector)(VecPack,Vec*);
24: int (*getcoloring)(VecPack,ISColoringType,ISColoring*);
25: int (*getmatrix)(VecPack,MatType,Mat*);
26: int (*getinterpolation)(VecPack,VecPack,Mat*,Vec*);
27: int (*refine)(VecPack,MPI_Comm,VecPack*);
28: };
30: struct _p_VecPack {
31: PETSCHEADER(struct _VecPackOps)
32: int rank;
33: int n,N,rstart; /* rstart is relative to all processors */
34: Vec globalvector;
35: int nDA,nredundant;
36: struct VecPackLink *next;
37: };
39: #undef __FUNCT__
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()
58: @*/
59: int VecPackCreate(MPI_Comm comm,VecPack *packer)
60: {
61: int ierr;
62: VecPack p;
66: *packer = PETSC_NULL;
67: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
68: DMInitializePackage(PETSC_NULL);
69: #endif
71: PetscHeaderCreate(p,_p_VecPack,struct _VecPackOps,DA_COOKIE,0,"VecPack",comm,VecPackDestroy,0);
72: PetscLogObjectCreate(p);
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: ierr = 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: }
88: #undef __FUNCT__
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: int VecPackDestroy(VecPack packer)
105: {
106: int 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: /* --------------------------------------------------------------------------------------*/
128: #undef __FUNCT__
130: int VecPackGetAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
131: {
132: int ierr;
133: PetscScalar *varray;
136: if (array) {
137: if (!packer->rank) {
138: ierr = VecGetArray(vec,&varray);
139: *array = varray + mine->rstart;
140: ierr = VecRestoreArray(vec,&varray);
141: } else {
142: *array = 0;
143: }
144: }
145: return(0);
146: }
148: #undef __FUNCT__
150: int VecPackGetAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
151: {
152: int ierr;
153: PetscScalar *array;
156: if (global) {
157: ierr = DAGetGlobalVector(mine->da,global);
158: ierr = VecGetArray(vec,&array);
159: ierr = VecPlaceArray(*global,array+mine->rstart);
160: ierr = VecRestoreArray(vec,&array);
161: }
162: return(0);
163: }
165: #undef __FUNCT__
167: int VecPackRestoreAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
168: {
170: return(0);
171: }
173: #undef __FUNCT__
175: int VecPackRestoreAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
176: {
177: int ierr;
180: if (global) {
181: VecResetArray(*global);
182: DARestoreGlobalVector(mine->da,global);
183: }
184: return(0);
185: }
187: #undef __FUNCT__
189: int VecPackScatter_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
190: {
191: int ierr;
192: PetscScalar *varray;
196: if (!packer->rank) {
197: ierr = VecGetArray(vec,&varray);
198: ierr = PetscMemcpy(array,varray+mine->rstart,mine->n*sizeof(PetscScalar));
199: ierr = VecRestoreArray(vec,&varray);
200: }
201: ierr = MPI_Bcast(array,mine->n,MPIU_SCALAR,0,packer->comm);
202: return(0);
203: }
205: #undef __FUNCT__
207: int VecPackScatter_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
208: {
209: int ierr;
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: }
225: #undef __FUNCT__
227: int VecPackGather_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
228: {
229: int ierr;
230: PetscScalar *varray;
233: if (!packer->rank) {
234: ierr = VecGetArray(vec,&varray);
235: if (varray+mine->rstart == array) SETERRQ(1,"You need not VecPackGather() into objects obtained via VecPackGetAccess()");
236: ierr = PetscMemcpy(varray+mine->rstart,array,mine->n*sizeof(PetscScalar));
237: ierr = VecRestoreArray(vec,&varray);
238: }
239: return(0);
240: }
242: #undef __FUNCT__
244: int VecPackGather_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
245: {
246: int ierr;
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>
265: #undef __FUNCT__
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()
284: @*/
285: int VecPackGetAccess(VecPack packer,Vec gvec,...)
286: {
287: va_list Argp;
288: int ierr;
289: struct VecPackLink *next = packer->next;
292: if (!packer->globalvector) {
293: SETERRQ(1,"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: ierr = 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(1,"Cannot handle that object type yet");
309: }
310: next = next->next;
311: }
312: va_end(Argp);
313: return(0);
314: }
316: #undef __FUNCT__
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: int VecPackRestoreAccess(VecPack packer,Vec gvec,...)
337: {
338: va_list Argp;
339: int ierr;
340: struct VecPackLink *next = packer->next;
343: if (!packer->globalvector) {
344: SETERRQ(1,"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: ierr = 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(1,"Cannot handle that object type yet");
360: }
361: next = next->next;
362: }
363: va_end(Argp);
364: return(0);
365: }
367: #undef __FUNCT__
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()
384: @*/
385: int VecPackScatter(VecPack packer,Vec gvec,...)
386: {
387: va_list Argp;
388: int ierr;
389: struct VecPackLink *next = packer->next;
392: if (!packer->globalvector) {
393: SETERRQ(1,"Must first create global vector with VecPackCreateGlobalVector()");
394: }
396: /* loop over packed objects, handling one at at time */
397: va_start(Argp,gvec);
398: while (next) {
399: if (next->type == VECPACK_ARRAY) {
400: PetscScalar *array;
401: array = va_arg(Argp, PetscScalar*);
402: VecPackScatter_Array(packer,next,gvec,array);
403: } else if (next->type == VECPACK_DA) {
404: Vec vec;
405: vec = va_arg(Argp, Vec);
407: VecPackScatter_DA(packer,next,gvec,vec);
408: } else {
409: SETERRQ(1,"Cannot handle that object type yet");
410: }
411: next = next->next;
412: }
413: va_end(Argp);
414: return(0);
415: }
417: #undef __FUNCT__
419: /*@C
420: VecPackGather - Gathers into a global packed vector from its individual local vectors
422: Collective on VecPack
424: Input Parameter:
425: + packer - the packer object
426: . gvec - the global vector
427: - ... - the individual sequential objects (arrays or vectors)
428:
429: Level: advanced
431: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
432: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
434: @*/
435: int VecPackGather(VecPack packer,Vec gvec,...)
436: {
437: va_list Argp;
438: int ierr;
439: struct VecPackLink *next = packer->next;
442: if (!packer->globalvector) {
443: SETERRQ(1,"Must first create global vector with VecPackCreateGlobalVector()");
444: }
446: /* loop over packed objects, handling one at at time */
447: va_start(Argp,gvec);
448: while (next) {
449: if (next->type == VECPACK_ARRAY) {
450: PetscScalar *array;
451: array = va_arg(Argp, PetscScalar*);
452: ierr = VecPackGather_Array(packer,next,gvec,array);
453: } else if (next->type == VECPACK_DA) {
454: Vec vec;
455: vec = va_arg(Argp, Vec);
457: VecPackGather_DA(packer,next,gvec,vec);
458: } else {
459: SETERRQ(1,"Cannot handle that object type yet");
460: }
461: next = next->next;
462: }
463: va_end(Argp);
464: return(0);
465: }
467: #undef __FUNCT__
469: /*@C
470: VecPackAddArray - adds an "redundant" array to a VecPack. The array values will
471: be stored in part of the array on processor 0.
473: Collective on VecPack
475: Input Parameter:
476: + packer - the packer object
477: - n - the length of the array
478:
479: Level: advanced
481: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
482: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
484: @*/
485: int VecPackAddArray(VecPack packer,int n)
486: {
487: struct VecPackLink *mine,*next = packer->next;
491: if (packer->globalvector) {
492: SETERRQ(1,"Cannot add an array once you have called VecPackCreateGlobalVector()");
493: }
495: /* create new link */
496: ierr = PetscNew(struct VecPackLink,&mine);
497: mine->n = n;
498: mine->da = PETSC_NULL;
499: mine->type = VECPACK_ARRAY;
500: mine->next = PETSC_NULL;
501: if (!packer->rank) packer->n += n;
503: /* add to end of list */
504: if (!next) {
505: packer->next = mine;
506: } else {
507: while (next->next) next = next->next;
508: next->next = mine;
509: }
510: packer->nredundant++;
511: return(0);
512: }
514: #undef __FUNCT__
516: /*@C
517: VecPackAddDA - adds a DA vector to a VecPack
519: Collective on VecPack
521: Input Parameter:
522: + packer - the packer object
523: - da - the DA object
524:
525: Level: advanced
527: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
528: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
530: @*/
531: int VecPackAddDA(VecPack packer,DA da)
532: {
533: int ierr,n;
534: struct VecPackLink *mine,*next = packer->next;
535: Vec global;
538: if (packer->globalvector) {
539: SETERRQ(1,"Cannot add a DA once you have called VecPackCreateGlobalVector()");
540: }
542: /* create new link */
543: PetscNew(struct VecPackLink,&mine);
544: PetscObjectReference((PetscObject)da);
545: DAGetGlobalVector(da,&global);
546: VecGetLocalSize(global,&n);
547: DARestoreGlobalVector(da,&global);
548: mine->n = n;
549: mine->da = da;
550: mine->type = VECPACK_DA;
551: mine->next = PETSC_NULL;
552: packer->n += n;
554: /* add to end of list */
555: if (!next) {
556: packer->next = mine;
557: } else {
558: while (next->next) next = next->next;
559: next->next = mine;
560: }
561: packer->nDA++;
562: return(0);
563: }
565: #undef __FUNCT__
567: /*@C
568: VecPackCreateGlobalVector - Creates a vector of the correct size to be gathered into
569: by the packer.
571: Collective on VecPack
573: Input Parameter:
574: . packer - the packer object
576: Output Parameters:
577: . gvec - the global vector
579: Level: advanced
581: Notes: Once this has been created you cannot add additional arrays or vectors to be packed.
583: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
584: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
586: @*/
587: int VecPackCreateGlobalVector(VecPack packer,Vec *gvec)
588: {
589: int ierr,nprev = 0,rank;
590: struct VecPackLink *next = packer->next;
593: if (packer->globalvector) {
594: VecDuplicate(packer->globalvector,gvec);
595: } else {
596: VecCreateMPI(packer->comm,packer->n,PETSC_DETERMINE,gvec);
597: PetscObjectReference((PetscObject)*gvec);
598: packer->globalvector = *gvec;
600: VecGetSize(*gvec,&packer->N);
601: VecGetOwnershipRange(*gvec,&packer->rstart,PETSC_NULL);
602:
603: /* now set the rstart for each linked array/vector */
604: MPI_Comm_rank(packer->comm,&rank);
605: while (next) {
606: next->rstart = nprev;
607: if (!rank || next->type != VECPACK_ARRAY) nprev += next->n;
608: next = next->next;
609: }
610: }
611: return(0);
612: }
614: #undef __FUNCT__
616: /*@C
617: VecPackGetGlobalIndices - Gets the global indices for all the entries in the packed
618: vectors.
620: Collective on VecPack
622: Input Parameter:
623: . packer - the packer object
625: Output Parameters:
626: . idx - the individual indices for each packed vector/array
627:
628: Level: advanced
630: Notes:
631: The idx parameters should be freed by the calling routine with PetscFree()
633: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
634: VecPackGather(), VecPackCreate(), VecPackGetAccess()
636: @*/
637: int VecPackGetGlobalIndices(VecPack packer,...)
638: {
639: va_list Argp;
640: int ierr,i,**idx,n;
641: struct VecPackLink *next = packer->next;
642: Vec global,dglobal;
643: PF pf;
644: PetscScalar *array;
647: VecPackCreateGlobalVector(packer,&global);
649: /* put 0 to N-1 into the global vector */
650: PFCreate(PETSC_COMM_WORLD,1,1,&pf);
651: PFSetType(pf,PFIDENTITY,PETSC_NULL);
652: PFApplyVec(pf,PETSC_NULL,global);
653: PFDestroy(pf);
655: /* loop over packed objects, handling one at at time */
656: va_start(Argp,packer);
657: while (next) {
658: idx = va_arg(Argp, int**);
660: if (next->type == VECPACK_ARRAY) {
661:
662: PetscMalloc(next->n*sizeof(int),idx);
663: if (!packer->rank) {
664: ierr = VecGetArray(global,&array);
665: array += next->rstart;
666: for (i=0; i<next->n; i++) (*idx)[i] = (int)PetscRealPart(array[i]);
667: array -= next->rstart;
668: ierr = VecRestoreArray(global,&array);
669: }
670: MPI_Bcast(*idx,next->n,MPI_INT,0,packer->comm);
672: } else if (next->type == VECPACK_DA) {
673: Vec local;
675: ierr = DACreateLocalVector(next->da,&local);
676: ierr = VecGetArray(global,&array);
677: array += next->rstart;
678: ierr = DAGetGlobalVector(next->da,&dglobal);
679: ierr = VecPlaceArray(dglobal,array);
680: ierr = DAGlobalToLocalBegin(next->da,dglobal,INSERT_VALUES,local);
681: ierr = DAGlobalToLocalEnd(next->da,dglobal,INSERT_VALUES,local);
682: array -= next->rstart;
683: ierr = VecRestoreArray(global,&array);
684: ierr = VecResetArray(dglobal);
685: ierr = DARestoreGlobalVector(next->da,&dglobal);
687: ierr = VecGetArray(local,&array);
688: ierr = VecGetSize(local,&n);
689: ierr = PetscMalloc(n*sizeof(int),idx);
690: for (i=0; i<n; i++) (*idx)[i] = (int)PetscRealPart(array[i]);
691: ierr = VecRestoreArray(local,&array);
692: ierr = VecDestroy(local);
694: } else {
695: SETERRQ(1,"Cannot handle that object type yet");
696: }
697: next = next->next;
698: }
699: va_end(Argp);
700: VecDestroy(global);
701: return(0);
702: }
704: /* -------------------------------------------------------------------------------------*/
705: #undef __FUNCT__
707: int VecPackGetLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
708: {
712: PetscMalloc(mine->n*sizeof(PetscScalar),array);
713: return(0);
714: }
716: #undef __FUNCT__
718: int VecPackGetLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
719: {
720: int ierr;
722: DAGetLocalVector(mine->da,local);
723: return(0);
724: }
726: #undef __FUNCT__
728: int VecPackRestoreLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
729: {
732: PetscFree(*array);
733: return(0);
734: }
736: #undef __FUNCT__
738: int VecPackRestoreLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
739: {
740: int ierr;
742: DARestoreLocalVector(mine->da,local);
743: return(0);
744: }
746: #undef __FUNCT__
748: /*@C
749: VecPackGetLocalVectors - Gets local vectors and arrays for each part of a VecPack.'
750: Use VecPakcRestoreLocalVectors() to return them.
752: Collective on VecPack
754: Input Parameter:
755: . packer - the packer object
756:
757: Output Parameter:
758: . ... - the individual sequential objects (arrays or vectors)
759:
760: Level: advanced
762: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
763: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
764: VecPackRestoreLocalVectors()
766: @*/
767: int VecPackGetLocalVectors(VecPack packer,...)
768: {
769: va_list Argp;
770: int ierr;
771: struct VecPackLink *next = packer->next;
775: /* loop over packed objects, handling one at at time */
776: va_start(Argp,packer);
777: while (next) {
778: if (next->type == VECPACK_ARRAY) {
779: PetscScalar **array;
780: array = va_arg(Argp, PetscScalar**);
781: VecPackGetLocalVectors_Array(packer,next,array);
782: } else if (next->type == VECPACK_DA) {
783: Vec *vec;
784: vec = va_arg(Argp, Vec*);
785: VecPackGetLocalVectors_DA(packer,next,vec);
786: } else {
787: SETERRQ(1,"Cannot handle that object type yet");
788: }
789: next = next->next;
790: }
791: va_end(Argp);
792: return(0);
793: }
795: #undef __FUNCT__
797: /*@C
798: VecPackRestoreLocalVectors - Restores local vectors and arrays for each part of a VecPack.'
799: Use VecPakcRestoreLocalVectors() to return them.
801: Collective on VecPack
803: Input Parameter:
804: . packer - the packer object
805:
806: Output Parameter:
807: . ... - the individual sequential objects (arrays or vectors)
808:
809: Level: advanced
811: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
812: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
813: VecPackGetLocalVectors()
815: @*/
816: int VecPackRestoreLocalVectors(VecPack packer,...)
817: {
818: va_list Argp;
819: int ierr;
820: struct VecPackLink *next = packer->next;
824: /* loop over packed objects, handling one at at time */
825: va_start(Argp,packer);
826: while (next) {
827: if (next->type == VECPACK_ARRAY) {
828: PetscScalar **array;
829: array = va_arg(Argp, PetscScalar**);
830: VecPackRestoreLocalVectors_Array(packer,next,array);
831: } else if (next->type == VECPACK_DA) {
832: Vec *vec;
833: vec = va_arg(Argp, Vec*);
834: VecPackRestoreLocalVectors_DA(packer,next,vec);
835: } else {
836: SETERRQ(1,"Cannot handle that object type yet");
837: }
838: next = next->next;
839: }
840: va_end(Argp);
841: return(0);
842: }
844: /* -------------------------------------------------------------------------------------*/
845: #undef __FUNCT__
847: int VecPackGetEntries_Array(VecPack packer,struct VecPackLink *mine,int *n)
848: {
850: if (n) *n = mine->n;
851: return(0);
852: }
854: #undef __FUNCT__
856: int VecPackGetEntries_DA(VecPack packer,struct VecPackLink *mine,DA *da)
857: {
859: if (da) *da = mine->da;
860: return(0);
861: }
863: #undef __FUNCT__
865: /*@C
866: VecPackGetEntries - Gets the DA, redundant size, etc for each entry in a VecPack.
867: Use VecPackRestoreEntries() to return them.
869: Collective on VecPack
871: Input Parameter:
872: . packer - the packer object
873:
874: Output Parameter:
875: . ... - the individual entries, DAs or integer sizes)
876:
877: Level: advanced
879: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
880: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
881: VecPackRestoreLocalVectors(), VecPackGetLocalVectors(), VecPackRestoreEntries()
883: @*/
884: int VecPackGetEntries(VecPack packer,...)
885: {
886: va_list Argp;
887: int ierr;
888: struct VecPackLink *next = packer->next;
892: /* loop over packed objects, handling one at at time */
893: va_start(Argp,packer);
894: while (next) {
895: if (next->type == VECPACK_ARRAY) {
896: int *n;
897: n = va_arg(Argp, int*);
898: VecPackGetEntries_Array(packer,next,n);
899: } else if (next->type == VECPACK_DA) {
900: DA *da;
901: da = va_arg(Argp, DA*);
902: VecPackGetEntries_DA(packer,next,da);
903: } else {
904: SETERRQ(1,"Cannot handle that object type yet");
905: }
906: next = next->next;
907: }
908: va_end(Argp);
909: return(0);
910: }
912: #undef __FUNCT__
914: /*@C
915: VecPackRefine - Refines a VecPack by refining all of its DAs
917: Collective on VecPack
919: Input Parameters:
920: + packer - the packer object
921: - comm - communicator to contain the new DM object, usually PETSC_NULL
923: Output Parameter:
924: . fine - new packer
925:
926: Level: advanced
928: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
929: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
931: @*/
932: int VecPackRefine(VecPack packer,MPI_Comm comm,VecPack *fine)
933: {
934: int ierr;
935: struct VecPackLink *next = packer->next;
936: DA da;
939: VecPackCreate(comm,fine);
941: /* loop over packed objects, handling one at at time */
942: while (next) {
943: if (next->type == VECPACK_ARRAY) {
944: VecPackAddArray(*fine,next->n);
945: } else if (next->type == VECPACK_DA) {
946: DARefine(next->da,comm,&da);
947: VecPackAddDA(*fine,da);
948: PetscObjectDereference((PetscObject)da);
949: } else {
950: SETERRQ(1,"Cannot handle that object type yet");
951: }
952: next = next->next;
953: }
954: return(0);
955: }
957: #include petscmat.h
959: struct MatPackLink {
960: Mat A;
961: struct MatPackLink *next;
962: };
964: struct MatPack {
965: VecPack right,left;
966: struct MatPackLink *next;
967: };
969: #undef __FUNCT__
971: int MatMultBoth_Shell_Pack(Mat A,Vec x,Vec y,PetscTruth add)
972: {
973: struct MatPack *mpack;
974: struct VecPackLink *xnext,*ynext;
975: struct MatPackLink *anext;
976: PetscScalar *xarray,*yarray;
977: int ierr,i;
978: Vec xglobal,yglobal;
981: MatShellGetContext(A,(void**)&mpack);
982: xnext = mpack->right->next;
983: ynext = mpack->left->next;
984: anext = mpack->next;
986: while (xnext) {
987: if (xnext->type == VECPACK_ARRAY) {
988: if (!mpack->right->rank) {
989: ierr = VecGetArray(x,&xarray);
990: ierr = VecGetArray(y,&yarray);
991: if (add) {
992: for (i=0; i<xnext->n; i++) {
993: yarray[ynext->rstart+i] += xarray[xnext->rstart+i];
994: }
995: } else {
996: ierr = PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
997: }
998: ierr = VecRestoreArray(x,&xarray);
999: ierr = VecRestoreArray(y,&yarray);
1000: }
1001: } else if (xnext->type == VECPACK_DA) {
1002: ierr = VecGetArray(x,&xarray);
1003: ierr = VecGetArray(y,&yarray);
1004: ierr = DAGetGlobalVector(xnext->da,&xglobal);
1005: ierr = DAGetGlobalVector(ynext->da,&yglobal);
1006: ierr = VecPlaceArray(xglobal,xarray+xnext->rstart);
1007: ierr = VecPlaceArray(yglobal,yarray+ynext->rstart);
1008: if (add) {
1009: ierr = MatMultAdd(anext->A,xglobal,yglobal,yglobal);
1010: } else {
1011: ierr = MatMult(anext->A,xglobal,yglobal);
1012: }
1013: ierr = VecRestoreArray(x,&xarray);
1014: ierr = VecRestoreArray(y,&yarray);
1015: ierr = VecResetArray(xglobal);
1016: ierr = VecResetArray(yglobal);
1017: ierr = DARestoreGlobalVector(xnext->da,&xglobal);
1018: ierr = DARestoreGlobalVector(ynext->da,&yglobal);
1019: anext = anext->next;
1020: } else {
1021: SETERRQ(1,"Cannot handle that object type yet");
1022: }
1023: xnext = xnext->next;
1024: ynext = ynext->next;
1025: }
1026: return(0);
1027: }
1029: #undef __FUNCT__
1031: int MatMultAdd_Shell_Pack(Mat A,Vec x,Vec y,Vec z)
1032: {
1035: if (z != y) SETERRQ(1,"Handles y == z only");
1036: MatMultBoth_Shell_Pack(A,x,y,PETSC_TRUE);
1037: return(0);
1038: }
1040: #undef __FUNCT__
1042: int MatMult_Shell_Pack(Mat A,Vec x,Vec y)
1043: {
1046: MatMultBoth_Shell_Pack(A,x,y,PETSC_FALSE);
1047: return(0);
1048: }
1050: #undef __FUNCT__
1052: int MatMultTranspose_Shell_Pack(Mat A,Vec x,Vec y)
1053: {
1054: struct MatPack *mpack;
1055: struct VecPackLink *xnext,*ynext;
1056: struct MatPackLink *anext;
1057: PetscScalar *xarray,*yarray;
1058: int ierr;
1059: Vec xglobal,yglobal;
1062: ierr = MatShellGetContext(A,(void**)&mpack);
1063: xnext = mpack->left->next;
1064: ynext = mpack->right->next;
1065: anext = mpack->next;
1067: while (xnext) {
1068: if (xnext->type == VECPACK_ARRAY) {
1069: if (!mpack->right->rank) {
1070: ierr = VecGetArray(x,&xarray);
1071: ierr = VecGetArray(y,&yarray);
1072: ierr = PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1073: ierr = VecRestoreArray(x,&xarray);
1074: ierr = VecRestoreArray(y,&yarray);
1075: }
1076: } else if (xnext->type == VECPACK_DA) {
1077: ierr = VecGetArray(x,&xarray);
1078: ierr = VecGetArray(y,&yarray);
1079: ierr = DAGetGlobalVector(xnext->da,&xglobal);
1080: ierr = DAGetGlobalVector(ynext->da,&yglobal);
1081: ierr = VecPlaceArray(xglobal,xarray+xnext->rstart);
1082: ierr = VecPlaceArray(yglobal,yarray+ynext->rstart);
1083: ierr = MatMultTranspose(anext->A,xglobal,yglobal);
1084: ierr = VecRestoreArray(x,&xarray);
1085: ierr = VecRestoreArray(y,&yarray);
1086: ierr = VecResetArray(xglobal);
1087: ierr = VecResetArray(yglobal);
1088: ierr = DARestoreGlobalVector(xnext->da,&xglobal);
1089: ierr = DARestoreGlobalVector(ynext->da,&yglobal);
1090: anext = anext->next;
1091: } else {
1092: SETERRQ(1,"Cannot handle that object type yet");
1093: }
1094: xnext = xnext->next;
1095: ynext = ynext->next;
1096: }
1097: return(0);
1098: }
1100: #undef __FUNCT__
1102: int MatDestroy_Shell_Pack(Mat A)
1103: {
1104: struct MatPack *mpack;
1105: struct MatPackLink *anext,*oldanext;
1106: int ierr;
1109: ierr = MatShellGetContext(A,(void**)&mpack);
1110: anext = mpack->next;
1112: while (anext) {
1113: ierr = MatDestroy(anext->A);
1114: oldanext = anext;
1115: anext = anext->next;
1116: ierr = PetscFree(oldanext);
1117: }
1118: PetscFree(mpack);
1119: return(0);
1120: }
1122: #undef __FUNCT__
1124: /*@C
1125: VecPackGetInterpolation - GetInterpolations a VecPack by refining all of its DAs
1127: Collective on VecPack
1129: Input Parameters:
1130: + coarse - coarse grid packer
1131: - fine - fine grid packer
1133: Output Parameter:
1134: + A - interpolation matrix
1135: - v - scaling vector
1136:
1137: Level: advanced
1139: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
1140: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess()
1142: @*/
1143: int VecPackGetInterpolation(VecPack coarse,VecPack fine,Mat *A,Vec *v)
1144: {
1145: int ierr,m,n,M,N;
1146: struct VecPackLink *nextc = coarse->next;
1147: struct VecPackLink *nextf = fine->next;
1148: struct MatPackLink *nextmat,*pnextmat = 0;
1149: struct MatPack *mpack;
1150: Vec gcoarse,gfine;
1153: /* use global vectors only for determining matrix layout */
1154: VecPackCreateGlobalVector(coarse,&gcoarse);
1155: VecPackCreateGlobalVector(fine,&gfine);
1156: VecGetLocalSize(gcoarse,&n);
1157: VecGetLocalSize(gfine,&m);
1158: VecGetSize(gcoarse,&N);
1159: VecGetSize(gfine,&M);
1160: VecDestroy(gcoarse);
1161: VecDestroy(gfine);
1163: ierr = PetscNew(struct MatPack,&mpack);
1164: mpack->right = coarse;
1165: mpack->left = fine;
1166: ierr = MatCreateShell(fine->comm,m,n,M,N,mpack,A);
1167: ierr = MatShellSetOperation(*A,MATOP_MULT,(void(*)(void))MatMult_Shell_Pack);
1168: ierr = MatShellSetOperation(*A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell_Pack);
1169: ierr = MatShellSetOperation(*A,MATOP_MULT_ADD,(void(*)(void))MatMultAdd_Shell_Pack);
1170: ierr = MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell_Pack);
1172: /* loop over packed objects, handling one at at time */
1173: while (nextc) {
1174: if (nextc->type != nextf->type) SETERRQ(1,"Two VecPack have different layout");
1176: if (nextc->type == VECPACK_ARRAY) {
1177: ;
1178: } else if (nextc->type == VECPACK_DA) {
1179: ierr = PetscNew(struct MatPackLink,&nextmat);
1180: nextmat->next = 0;
1181: if (pnextmat) {
1182: pnextmat->next = nextmat;
1183: pnextmat = nextmat;
1184: } else {
1185: pnextmat = nextmat;
1186: mpack->next = nextmat;
1187: }
1188: DAGetInterpolation(nextc->da,nextf->da,&nextmat->A,PETSC_NULL);
1189: } else {
1190: SETERRQ(1,"Cannot handle that object type yet");
1191: }
1192: nextc = nextc->next;
1193: nextf = nextf->next;
1194: }
1195: return(0);
1196: }