Actual source code: pack.c
petsc-3.3-p7 2013-05-11
2: #include packimpl.h
6: /*@C
7: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
8: seperate components (DMs) in a DMto build the correct matrix nonzero structure.
11: Logically Collective on MPI_Comm
13: Input Parameter:
14: + dm - the composite object
15: - formcouplelocations - routine to set the nonzero locations in the matrix
17: Level: advanced
19: Notes: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
20: this routine
22: @*/
23: PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt))
24: {
25: DM_Composite *com = (DM_Composite*)dm->data;
28: com->FormCoupleLocations = FormCoupleLocations;
29: return(0);
30: }
34: PetscErrorCode DMDestroy_Composite(DM dm)
35: {
36: PetscErrorCode ierr;
37: struct DMCompositeLink *next, *prev;
38: DM_Composite *com = (DM_Composite *)dm->data;
41: next = com->next;
42: while (next) {
43: prev = next;
44: next = next->next;
45: DMDestroy(&prev->dm);
46: PetscFree(prev->grstarts);
47: PetscFree(prev);
48: }
49: return(0);
50: }
54: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
55: {
57: PetscBool iascii;
58: DM_Composite *com = (DM_Composite *)dm->data;
61: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
62: if (iascii) {
63: struct DMCompositeLink *lnk = com->next;
64: PetscInt i;
66: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix?((PetscObject)dm)->prefix:"no prefix");
67: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
68: PetscViewerASCIIPushTab(v);
69: for (i=0; lnk; lnk=lnk->next,i++) {
70: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
71: PetscViewerASCIIPushTab(v);
72: DMView(lnk->dm,v);
73: PetscViewerASCIIPopTab(v);
74: }
75: PetscViewerASCIIPopTab(v);
76: }
77: return(0);
78: }
80: /* --------------------------------------------------------------------------------------*/
83: PetscErrorCode DMSetUp_Composite(DM dm)
84: {
85: PetscErrorCode ierr;
86: PetscInt nprev = 0;
87: PetscMPIInt rank,size;
88: DM_Composite *com = (DM_Composite*)dm->data;
89: struct DMCompositeLink *next = com->next;
90: PetscLayout map;
93: if (com->setup) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
94: PetscLayoutCreate(((PetscObject)dm)->comm,&map);
95: PetscLayoutSetLocalSize(map,com->n);
96: PetscLayoutSetSize(map,PETSC_DETERMINE);
97: PetscLayoutSetBlockSize(map,1);
98: PetscLayoutSetUp(map);
99: PetscLayoutGetSize(map,&com->N);
100: PetscLayoutGetRange(map,&com->rstart,PETSC_NULL);
101: PetscLayoutDestroy(&map);
102:
103: /* now set the rstart for each linked vector */
104: MPI_Comm_rank(((PetscObject)dm)->comm,&rank);
105: MPI_Comm_size(((PetscObject)dm)->comm,&size);
106: while (next) {
107: next->rstart = nprev;
108: nprev += next->n;
109: next->grstart = com->rstart + next->rstart;
110: PetscMalloc(size*sizeof(PetscInt),&next->grstarts);
111: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,((PetscObject)dm)->comm);
112: next = next->next;
113: }
114: com->setup = PETSC_TRUE;
115: return(0);
116: }
118: /* ----------------------------------------------------------------------------------*/
120: #include <stdarg.h>
124: /*@C
125: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
126: representation.
128: Not Collective
130: Input Parameter:
131: . dm - the packer object
133: Output Parameter:
134: . nDM - the number of DMs
136: Level: beginner
138: @*/
139: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
140: {
141: DM_Composite *com = (DM_Composite*)dm->data;
144: *nDM = com->nDM;
145: return(0);
146: }
151: /*@C
152: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
153: representation.
155: Collective on DMComposite
157: Input Parameters:
158: + dm - the packer object
159: - gvec - the global vector
161: Output Parameters:
162: . Vec* ... - the packed parallel vectors, PETSC_NULL for those that are not needed
164: Notes: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
165:
166: Level: advanced
168: @*/
169: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
170: {
171: va_list Argp;
172: PetscErrorCode ierr;
173: struct DMCompositeLink *next;
174: DM_Composite *com = (DM_Composite*)dm->data;
179: next = com->next;
180: if (!com->setup) {
181: DMSetUp(dm);
182: }
184: /* loop over packed objects, handling one at at time */
185: va_start(Argp,gvec);
186: while (next) {
187: Vec *vec;
188: vec = va_arg(Argp, Vec*);
189: if (vec) {
190: PetscScalar *array;
191: DMGetGlobalVector(next->dm,vec);
192: VecGetArray(gvec,&array);
193: VecPlaceArray(*vec,array+next->rstart);
194: VecRestoreArray(gvec,&array);
195: }
196: next = next->next;
197: }
198: va_end(Argp);
199: return(0);
200: }
204: /*@C
205: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
206: representation.
208: Collective on DMComposite
210: Input Parameters:
211: + dm - the packer object
212: . gvec - the global vector
213: - Vec* ... - the individual parallel vectors, PETSC_NULL for those that are not needed
214:
215: Level: advanced
217: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
218: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
219: DMCompositeRestoreAccess(), DMCompositeGetAccess()
221: @*/
222: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
223: {
224: va_list Argp;
225: PetscErrorCode ierr;
226: struct DMCompositeLink *next;
227: DM_Composite *com = (DM_Composite*)dm->data;
232: next = com->next;
233: if (!com->setup) {
234: DMSetUp(dm);
235: }
237: /* loop over packed objects, handling one at at time */
238: va_start(Argp,gvec);
239: while (next) {
240: Vec *vec;
241: vec = va_arg(Argp, Vec*);
242: if (vec) {
243: VecResetArray(*vec);
244: DMRestoreGlobalVector(next->dm,vec);
245: }
246: next = next->next;
247: }
248: va_end(Argp);
249: return(0);
250: }
254: /*@C
255: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
257: Collective on DMComposite
259: Input Parameters:
260: + dm - the packer object
261: . gvec - the global vector
262: - Vec ... - the individual sequential vectors, PETSC_NULL for those that are not needed
264: Level: advanced
266: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
267: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
268: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
270: @*/
271: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
272: {
273: va_list Argp;
274: PetscErrorCode ierr;
275: struct DMCompositeLink *next;
276: PetscInt cnt;
277: DM_Composite *com = (DM_Composite*)dm->data;
282: if (!com->setup) {
283: DMSetUp(dm);
284: }
286: /* loop over packed objects, handling one at at time */
287: va_start(Argp,gvec);
288: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
289: Vec local;
290: local = va_arg(Argp, Vec);
291: if (local) {
292: Vec global;
293: PetscScalar *array;
295: DMGetGlobalVector(next->dm,&global);
296: VecGetArray(gvec,&array);
297: VecPlaceArray(global,array+next->rstart);
298: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
299: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
300: VecRestoreArray(gvec,&array);
301: VecResetArray(global);
302: DMRestoreGlobalVector(next->dm,&global);
303: }
304: }
305: va_end(Argp);
306: return(0);
307: }
311: /*@C
312: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
314: Collective on DMComposite
316: Input Parameter:
317: + dm - the packer object
318: . gvec - the global vector
319: - Vec ... - the individual sequential vectors, PETSC_NULL for any that are not needed
321: Level: advanced
323: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
324: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
325: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
327: @*/
328: PetscErrorCode DMCompositeGather(DM dm,Vec gvec,InsertMode imode,...)
329: {
330: va_list Argp;
331: PetscErrorCode ierr;
332: struct DMCompositeLink *next;
333: DM_Composite *com = (DM_Composite*)dm->data;
334: PetscInt cnt;
339: if (!com->setup) {
340: DMSetUp(dm);
341: }
343: /* loop over packed objects, handling one at at time */
344: va_start(Argp,imode);
345: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
346: Vec local;
347: local = va_arg(Argp, Vec);
348: if (local) {
349: PetscScalar *array;
350: Vec global;
352: DMGetGlobalVector(next->dm,&global);
353: VecGetArray(gvec,&array);
354: VecPlaceArray(global,array+next->rstart);
355: DMLocalToGlobalBegin(next->dm,local,imode,global);
356: DMLocalToGlobalEnd(next->dm,local,imode,global);
357: VecRestoreArray(gvec,&array);
358: VecResetArray(global);
359: DMRestoreGlobalVector(next->dm,&global);
360: }
361: }
362: va_end(Argp);
363: return(0);
364: }
368: /*@C
369: DMCompositeAddDM - adds a DM vector to a DMComposite
371: Collective on DMComposite
373: Input Parameter:
374: + dm - the packer object
375: - dm - the DM object, if the DM is a da you will need to caste it with a (DM)
376:
377: Level: advanced
379: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
380: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
381: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
383: @*/
384: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
385: {
386: PetscErrorCode ierr;
387: PetscInt n,nlocal;
388: struct DMCompositeLink *mine,*next;
389: Vec global,local;
390: DM_Composite *com = (DM_Composite*)dmc->data;
395: next = com->next;
396: if (com->setup) SETERRQ(((PetscObject)dmc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
398: /* create new link */
399: PetscNew(struct DMCompositeLink,&mine);
400: PetscObjectReference((PetscObject)dm);
401: DMGetGlobalVector(dm,&global);
402: VecGetLocalSize(global,&n);
403: DMRestoreGlobalVector(dm,&global);
404: DMGetLocalVector(dm,&local);
405: VecGetSize(local,&nlocal);
406: DMRestoreLocalVector(dm,&local);
407: mine->n = n;
408: mine->nlocal = nlocal;
409: mine->dm = dm;
410: mine->next = PETSC_NULL;
411: com->n += n;
413: /* add to end of list */
414: if (!next) {
415: com->next = mine;
416: } else {
417: while (next->next) next = next->next;
418: next->next = mine;
419: }
420: com->nDM++;
421: com->nmine++;
422: return(0);
423: }
425: extern PetscErrorCode VecView_MPI(Vec,PetscViewer);
426: EXTERN_C_BEGIN
429: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
430: {
431: DM dm;
432: PetscErrorCode ierr;
433: struct DMCompositeLink *next;
434: PetscBool isdraw;
435: DM_Composite *com;
438: PetscObjectQuery((PetscObject)gvec,"DM",(PetscObject*)&dm);
439: if (!dm) SETERRQ(((PetscObject)gvec)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
440: com = (DM_Composite*)dm->data;
441: next = com->next;
443: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
444: if (!isdraw) {
445: /* do I really want to call this? */
446: VecView_MPI(gvec,viewer);
447: } else {
448: PetscInt cnt = 0;
450: /* loop over packed objects, handling one at at time */
451: while (next) {
452: Vec vec;
453: PetscScalar *array;
454: PetscInt bs;
456: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
457: DMGetGlobalVector(next->dm,&vec);
458: VecGetArray(gvec,&array);
459: VecPlaceArray(vec,array+next->rstart);
460: VecRestoreArray(gvec,&array);
461: VecView(vec,viewer);
462: VecGetBlockSize(vec,&bs);
463: VecResetArray(vec);
464: DMRestoreGlobalVector(next->dm,&vec);
465: PetscViewerDrawBaseAdd(viewer,bs);
466: cnt += bs;
467: next = next->next;
468: }
469: PetscViewerDrawBaseAdd(viewer,-cnt);
470: }
471: return(0);
472: }
473: EXTERN_C_END
478: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
479: {
480: PetscErrorCode ierr;
481: DM_Composite *com = (DM_Composite*)dm->data;
485: DMSetUp(dm);
486: VecCreateMPI(((PetscObject)dm)->comm,com->n,com->N,gvec);
487: PetscObjectCompose((PetscObject)*gvec,"DM",(PetscObject)dm);
488: VecSetOperation(*gvec,VECOP_VIEW,(void(*)(void))VecView_DMComposite);
489: return(0);
490: }
494: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
495: {
496: PetscErrorCode ierr;
497: DM_Composite *com = (DM_Composite*)dm->data;
501: if (!com->setup) {
502: DMSetUp(dm);
503: }
504: VecCreateSeq(((PetscObject)dm)->comm,com->nghost,lvec);
505: PetscObjectCompose((PetscObject)*lvec,"DM",(PetscObject)dm);
506: return(0);
507: }
511: /*@C
512: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
514: Collective on DM
516: Input Parameter:
517: . dm - the packer object
519: Output Parameters:
520: . ltogs - the individual mappings for each packed vector. Note that this includes
521: all the ghost points that individual ghosted DMDA's may have.
523: Level: advanced
525: Notes:
526: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
528: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
529: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
530: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
532: @*/
533: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
534: {
535: PetscErrorCode ierr;
536: PetscInt i,*idx,n,cnt;
537: struct DMCompositeLink *next;
538: PetscMPIInt rank;
539: DM_Composite *com = (DM_Composite*)dm->data;
543: DMSetUp(dm);
544: PetscMalloc((com->nDM)*sizeof(ISLocalToGlobalMapping),ltogs);
545: next = com->next;
546: MPI_Comm_rank(((PetscObject)dm)->comm,&rank);
548: /* loop over packed objects, handling one at at time */
549: cnt = 0;
550: while (next) {
551: ISLocalToGlobalMapping ltog;
552: PetscMPIInt size;
553: const PetscInt *suboff,*indices;
554: Vec global;
556: /* Get sub-DM global indices for each local dof */
557: DMGetLocalToGlobalMapping(next->dm,<og);
558: ISLocalToGlobalMappingGetSize(ltog,&n);
559: ISLocalToGlobalMappingGetIndices(ltog,&indices);
560: PetscMalloc(n*sizeof(PetscInt),&idx);
562: /* Get the offsets for the sub-DM global vector */
563: DMGetGlobalVector(next->dm,&global);
564: VecGetOwnershipRanges(global,&suboff);
565: MPI_Comm_size(((PetscObject)global)->comm,&size);
567: /* Shift the sub-DM definition of the global space to the composite global space */
568: for (i=0; i<n; i++) {
569: PetscInt subi = indices[i],lo = 0,hi = size,t;
570: /* Binary search to find which rank owns subi */
571: while (hi-lo > 1) {
572: t = lo + (hi-lo)/2;
573: if (suboff[t] > subi) hi = t;
574: else lo = t;
575: }
576: idx[i] = subi - suboff[lo] + next->grstarts[lo];
577: }
578: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
579: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
580: DMRestoreGlobalVector(next->dm,&global);
581: next = next->next;
582: cnt++;
583: }
584: return(0);
585: }
589: /*@C
590: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
592: Not Collective
594: Input Arguments:
595: . dm - composite DM
597: Output Arguments:
598: . is - array of serial index sets for each each component of the DMComposite
600: Level: intermediate
602: Notes:
603: At present, a composite local vector does not normally exist. This function is used to provide index sets for
604: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
605: scatter to a composite local vector. The user should not typically need to know which is being done.
607: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
609: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
611: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
613: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
614: @*/
615: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
616: {
617: PetscErrorCode ierr;
618: DM_Composite *com = (DM_Composite*)dm->data;
619: struct DMCompositeLink *link;
620: PetscInt cnt,start;
625: PetscMalloc(com->nmine*sizeof(IS),is);
626: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
627: PetscInt bs;
628: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
629: DMGetBlockSize(link->dm,&bs);
630: ISSetBlockSize((*is)[cnt],bs);
631: }
632: return(0);
633: }
637: /*@C
638: DMCompositeGetGlobalISs - Gets the index sets for each composed object
640: Collective on DMComposite
642: Input Parameter:
643: . dm - the packer object
645: Output Parameters:
646: . is - the array of index sets
648: Level: advanced
650: Notes:
651: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
653: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
655: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
656: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
657: indices.
659: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
660: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
661: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
663: @*/
665: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
666: {
667: PetscErrorCode ierr;
668: PetscInt cnt = 0,*idx,i;
669: struct DMCompositeLink *next;
670: PetscMPIInt rank;
671: DM_Composite *com = (DM_Composite*)dm->data;
675: PetscMalloc((com->nDM)*sizeof(IS),is);
676: next = com->next;
677: MPI_Comm_rank(((PetscObject)dm)->comm,&rank);
679: /* loop over packed objects, handling one at at time */
680: while (next) {
681: PetscMalloc(next->n*sizeof(PetscInt),&idx);
682: for (i=0; i<next->n; i++) idx[i] = next->grstart + i;
683: ISCreateGeneral(((PetscObject)dm)->comm,next->n,idx,PETSC_OWN_POINTER,&(*is)[cnt]);
684: cnt++;
685: next = next->next;
686: }
687: return(0);
688: }
692: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
693: {
694: PetscInt nDM;
695: DM *dms;
696: PetscInt i;
700: DMCompositeGetNumberDM(dm, &nDM);
701: if (numFields) {*numFields = nDM;}
702: DMCompositeGetGlobalISs(dm, fields);
703: if (fieldNames) {
704: PetscMalloc(nDM*sizeof(DM), &dms);
705: PetscMalloc(nDM*sizeof(const char *), fieldNames);
706: DMCompositeGetEntriesArray(dm, dms);
707: for (i=0; i<nDM; i++) {
708: char buf[256];
709: const char *splitname;
711: /* Split naming precedence: object name, prefix, number */
712: splitname = ((PetscObject) dm)->name;
713: if (!splitname) {
714: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
715: if (splitname) {
716: size_t len;
717: PetscStrncpy(buf,splitname,sizeof buf);
718: buf[sizeof buf - 1] = 0;
719: PetscStrlen(buf,&len);
720: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
721: splitname = buf;
722: }
723: }
724: if (!splitname) {
725: PetscSNPrintf(buf,sizeof buf,"%D",i);
726: splitname = buf;
727: }
728: PetscStrallocpy(splitname,&(*fieldNames)[i]);
729: }
730: PetscFree(dms);
731: }
732: return(0);
733: }
735: /*
736: This could take over from DMCreateFieldIS(), as it is more general,
737: making DMCreateFieldIS() a special case -- calling with dmlist == PETSC_NULL;
738: At this point it's probably best to be less intrusive, however.
739: */
742: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist)
743: {
744: PetscInt nDM;
745: PetscInt i;
749: DMCreateFieldIS_Composite(dm, len, namelist, islist);
750: if(dmlist) {
751: DMCompositeGetNumberDM(dm, &nDM);
752: PetscMalloc(nDM*sizeof(DM), dmlist);
753: DMCompositeGetEntriesArray(dm, *dmlist);
754: for (i=0; i<nDM; i++) {
755: PetscObjectReference((PetscObject)((*dmlist)[i]));
756: }
757: }
758: return(0);
759: }
763: /* -------------------------------------------------------------------------------------*/
766: /*@C
767: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
768: Use DMCompositeRestoreLocalVectors() to return them.
770: Not Collective
772: Input Parameter:
773: . dm - the packer object
774:
775: Output Parameter:
776: . Vec ... - the individual sequential Vecs
777:
778: Level: advanced
780: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
781: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
782: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
784: @*/
785: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
786: {
787: va_list Argp;
788: PetscErrorCode ierr;
789: struct DMCompositeLink *next;
790: DM_Composite *com = (DM_Composite*)dm->data;
794: next = com->next;
795: /* loop over packed objects, handling one at at time */
796: va_start(Argp,dm);
797: while (next) {
798: Vec *vec;
799: vec = va_arg(Argp, Vec*);
800: if (vec) {DMGetLocalVector(next->dm,vec);}
801: next = next->next;
802: }
803: va_end(Argp);
804: return(0);
805: }
809: /*@C
810: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
812: Not Collective
814: Input Parameter:
815: . dm - the packer object
816:
817: Output Parameter:
818: . Vec ... - the individual sequential Vecs
819:
820: Level: advanced
822: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
823: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
824: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
826: @*/
827: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
828: {
829: va_list Argp;
830: PetscErrorCode ierr;
831: struct DMCompositeLink *next;
832: DM_Composite *com = (DM_Composite*)dm->data;
836: next = com->next;
837: /* loop over packed objects, handling one at at time */
838: va_start(Argp,dm);
839: while (next) {
840: Vec *vec;
841: vec = va_arg(Argp, Vec*);
842: if (vec) {DMRestoreLocalVector(next->dm,vec);}
843: next = next->next;
844: }
845: va_end(Argp);
846: return(0);
847: }
849: /* -------------------------------------------------------------------------------------*/
852: /*@C
853: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
855: Not Collective
857: Input Parameter:
858: . dm - the packer object
859:
860: Output Parameter:
861: . DM ... - the individual entries (DMs)
862:
863: Level: advanced
865: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
866: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
867: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
868: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
870: @*/
871: PetscErrorCode DMCompositeGetEntries(DM dm,...)
872: {
873: va_list Argp;
874: struct DMCompositeLink *next;
875: DM_Composite *com = (DM_Composite*)dm->data;
879: next = com->next;
880: /* loop over packed objects, handling one at at time */
881: va_start(Argp,dm);
882: while (next) {
883: DM *dmn;
884: dmn = va_arg(Argp, DM*);
885: if (dmn) *dmn = next->dm;
886: next = next->next;
887: }
888: va_end(Argp);
889: return(0);
890: }
894: /*@
895: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
897: Not Collective
899: Input Parameter:
900: + dm - the packer object
901: - dms - array of sufficient length (see DMCompositeGetNumberDM()), holds the DMs on output
903: Level: advanced
905: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
906: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
907: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
908: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
910: @*/
911: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
912: {
913: struct DMCompositeLink *next;
914: DM_Composite *com = (DM_Composite*)dm->data;
915: PetscInt i;
919: /* loop over packed objects, handling one at at time */
920: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
921: return(0);
922: }
926: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
927: {
928: PetscErrorCode ierr;
929: struct DMCompositeLink *next;
930: DM_Composite *com = (DM_Composite*)dmi->data;
931: DM dm;
935: if (comm == MPI_COMM_NULL) comm = ((PetscObject)dmi)->comm;
936: DMSetUp(dmi);
937: next = com->next;
938: DMCompositeCreate(comm,fine);
940: /* loop over packed objects, handling one at at time */
941: while (next) {
942: DMRefine(next->dm,comm,&dm);
943: DMCompositeAddDM(*fine,dm);
944: PetscObjectDereference((PetscObject)dm);
945: next = next->next;
946: }
947: return(0);
948: }
952: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
953: {
954: PetscErrorCode ierr;
955: struct DMCompositeLink *next;
956: DM_Composite *com = (DM_Composite*)dmi->data;
957: DM dm;
961: DMSetUp(dmi);
962: if (comm == MPI_COMM_NULL) {
963: PetscObjectGetComm((PetscObject)dmi,&comm);
964: }
965: next = com->next;
966: DMCompositeCreate(comm,fine);
968: /* loop over packed objects, handling one at at time */
969: while (next) {
970: DMCoarsen(next->dm,comm,&dm);
971: DMCompositeAddDM(*fine,dm);
972: PetscObjectDereference((PetscObject)dm);
973: next = next->next;
974: }
975: return(0);
976: }
980: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
981: {
982: PetscErrorCode ierr;
983: PetscInt m,n,M,N,nDM,i;
984: struct DMCompositeLink *nextc;
985: struct DMCompositeLink *nextf;
986: Vec gcoarse,gfine,*vecs;
987: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
988: DM_Composite *comfine = (DM_Composite*)fine->data;
989: Mat *mats;
994: DMSetUp(coarse);
995: DMSetUp(fine);
996: /* use global vectors only for determining matrix layout */
997: DMGetGlobalVector(coarse,&gcoarse);
998: DMGetGlobalVector(fine,&gfine);
999: VecGetLocalSize(gcoarse,&n);
1000: VecGetLocalSize(gfine,&m);
1001: VecGetSize(gcoarse,&N);
1002: VecGetSize(gfine,&M);
1003: DMRestoreGlobalVector(coarse,&gcoarse);
1004: DMRestoreGlobalVector(fine,&gfine);
1006: nDM = comfine->nDM;
1007: if (nDM != comcoarse->nDM) SETERRQ2(((PetscObject)fine)->comm,PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1008: PetscMalloc(nDM*nDM*sizeof(Mat),&mats);
1009: PetscMemzero(mats,nDM*nDM*sizeof(Mat));
1010: if (v) {
1011: PetscMalloc(nDM*sizeof(Vec),&vecs);
1012: PetscMemzero(vecs,nDM*sizeof(Vec));
1013: }
1015: /* loop over packed objects, handling one at at time */
1016: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1017: if (!v) {
1018: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],PETSC_NULL);
1019: } else {
1020: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1021: }
1022: }
1023: MatCreateNest(((PetscObject)fine)->comm,nDM,PETSC_NULL,nDM,PETSC_NULL,mats,A);
1024: if (v) {
1025: VecCreateNest(((PetscObject)fine)->comm,nDM,PETSC_NULL,vecs,v);
1026: }
1027: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1028: PetscFree(mats);
1029: if (v) {
1030: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1031: PetscFree(vecs);
1032: }
1033: return(0);
1034: }
1038: static PetscErrorCode DMCreateLocalToGlobalMapping_Composite(DM dm)
1039: {
1040: DM_Composite *com = (DM_Composite*)dm->data;
1041: ISLocalToGlobalMapping *ltogs;
1042: PetscInt i;
1043: PetscErrorCode ierr;
1046: /* Set the ISLocalToGlobalMapping on the new matrix */
1047: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1048: ISLocalToGlobalMappingConcatenate(((PetscObject)dm)->comm,com->nDM,ltogs,&dm->ltogmap);
1049: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1050: PetscFree(ltogs);
1051: return(0);
1052: }
1057: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,const MatType mtype,ISColoring *coloring)
1058: {
1059: PetscErrorCode ierr;
1060: PetscInt n,i,cnt;
1061: ISColoringValue *colors;
1062: PetscBool dense = PETSC_FALSE;
1063: ISColoringValue maxcol = 0;
1064: DM_Composite *com = (DM_Composite*)dm->data;
1068: if (ctype == IS_COLORING_GHOSTED) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Only global coloring supported" );
1069: else if (ctype == IS_COLORING_GLOBAL) {
1070: n = com->n;
1071: } else SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1072: PetscMalloc(n*sizeof(ISColoringValue),&colors); /* freed in ISColoringDestroy() */
1074: PetscOptionsGetBool(PETSC_NULL,"-dmcomposite_dense_jacobian",&dense,PETSC_NULL);
1075: if (dense) {
1076: for (i=0; i<n; i++) {
1077: colors[i] = (ISColoringValue)(com->rstart + i);
1078: }
1079: maxcol = com->N;
1080: } else {
1081: struct DMCompositeLink *next = com->next;
1082: PetscMPIInt rank;
1083:
1084: MPI_Comm_rank(((PetscObject)dm)->comm,&rank);
1085: cnt = 0;
1086: while (next) {
1087: ISColoring lcoloring;
1089: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,mtype,&lcoloring);
1090: for (i=0; i<lcoloring->N; i++) {
1091: colors[cnt++] = maxcol + lcoloring->colors[i];
1092: }
1093: maxcol += lcoloring->n;
1094: ISColoringDestroy(&lcoloring);
1095: next = next->next;
1096: }
1097: }
1098: ISColoringCreate(((PetscObject)dm)->comm,maxcol,n,colors,coloring);
1099: return(0);
1100: }
1104: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1105: {
1106: PetscErrorCode ierr;
1107: struct DMCompositeLink *next;
1108: PetscInt cnt = 3;
1109: PetscMPIInt rank;
1110: PetscScalar *garray,*larray;
1111: DM_Composite *com = (DM_Composite*)dm->data;
1116: next = com->next;
1117: if (!com->setup) {
1118: DMSetUp(dm);
1119: }
1120: MPI_Comm_rank(((PetscObject)dm)->comm,&rank);
1121: VecGetArray(gvec,&garray);
1122: VecGetArray(lvec,&larray);
1124: /* loop over packed objects, handling one at at time */
1125: while (next) {
1126: Vec local,global;
1127: PetscInt N;
1129: DMGetGlobalVector(next->dm,&global);
1130: VecGetLocalSize(global,&N);
1131: VecPlaceArray(global,garray);
1132: DMGetLocalVector(next->dm,&local);
1133: VecPlaceArray(local,larray);
1134: DMGlobalToLocalBegin(next->dm,global,mode,local);
1135: DMGlobalToLocalEnd(next->dm,global,mode,local);
1136: VecResetArray(global);
1137: VecResetArray(local);
1138: DMRestoreGlobalVector(next->dm,&global);
1139: DMRestoreGlobalVector(next->dm,&local);
1140: cnt++;
1141: larray += next->nlocal;
1142: next = next->next;
1143: }
1145: VecRestoreArray(gvec,PETSC_NULL);
1146: VecRestoreArray(lvec,PETSC_NULL);
1147: return(0);
1148: }
1152: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1153: {
1155: return(0);
1156: }
1158: EXTERN_C_BEGIN
1161: PetscErrorCode DMCreate_Composite(DM p)
1162: {
1164: DM_Composite *com;
1167: PetscNewLog(p,DM_Composite,&com);
1168: p->data = com;
1169: PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1170: com->n = 0;
1171: com->next = PETSC_NULL;
1172: com->nDM = 0;
1174: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1175: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1176: p->ops->createlocaltoglobalmapping = DMCreateLocalToGlobalMapping_Composite;
1177: p->ops->createlocaltoglobalmappingblock = 0;
1178: p->ops->createfieldis = DMCreateFieldIS_Composite;
1179: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1180: p->ops->refine = DMRefine_Composite;
1181: p->ops->coarsen = DMCoarsen_Composite;
1182: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1183: p->ops->creatematrix = DMCreateMatrix_Composite;
1184: p->ops->getcoloring = DMCreateColoring_Composite;
1185: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1186: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1187: p->ops->destroy = DMDestroy_Composite;
1188: p->ops->view = DMView_Composite;
1189: p->ops->setup = DMSetUp_Composite;
1190: return(0);
1191: }
1192: EXTERN_C_END
1196: /*@C
1197: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1198: vectors made up of several subvectors.
1200: Collective on MPI_Comm
1202: Input Parameter:
1203: . comm - the processors that will share the global vector
1205: Output Parameters:
1206: . packer - the packer object
1208: Level: advanced
1210: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(),
1211: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1212: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1214: @*/
1215: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1216: {
1221: DMCreate(comm,packer);
1222: DMSetType(*packer,DMCOMPOSITE);
1223: return(0);
1224: }