Actual source code: pack.c
petsc-3.8.4 2018-03-24
2: #include <../src/dm/impls/composite/packimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petscds.h>
6: /*@C
7: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
8: separate 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: }
32: PetscErrorCode DMDestroy_Composite(DM dm)
33: {
34: PetscErrorCode ierr;
35: struct DMCompositeLink *next, *prev;
36: DM_Composite *com = (DM_Composite*)dm->data;
39: next = com->next;
40: while (next) {
41: prev = next;
42: next = next->next;
43: DMDestroy(&prev->dm);
44: PetscFree(prev->grstarts);
45: PetscFree(prev);
46: }
47: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
48: PetscFree(com);
49: return(0);
50: }
52: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
53: {
55: PetscBool iascii;
56: DM_Composite *com = (DM_Composite*)dm->data;
59: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
60: if (iascii) {
61: struct DMCompositeLink *lnk = com->next;
62: PetscInt i;
64: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
65: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
66: PetscViewerASCIIPushTab(v);
67: for (i=0; lnk; lnk=lnk->next,i++) {
68: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
69: PetscViewerASCIIPushTab(v);
70: DMView(lnk->dm,v);
71: PetscViewerASCIIPopTab(v);
72: }
73: PetscViewerASCIIPopTab(v);
74: }
75: return(0);
76: }
78: /* --------------------------------------------------------------------------------------*/
79: PetscErrorCode DMSetUp_Composite(DM dm)
80: {
81: PetscErrorCode ierr;
82: PetscInt nprev = 0;
83: PetscMPIInt rank,size;
84: DM_Composite *com = (DM_Composite*)dm->data;
85: struct DMCompositeLink *next = com->next;
86: PetscLayout map;
89: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
90: PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
91: PetscLayoutSetLocalSize(map,com->n);
92: PetscLayoutSetSize(map,PETSC_DETERMINE);
93: PetscLayoutSetBlockSize(map,1);
94: PetscLayoutSetUp(map);
95: PetscLayoutGetSize(map,&com->N);
96: PetscLayoutGetRange(map,&com->rstart,NULL);
97: PetscLayoutDestroy(&map);
99: /* now set the rstart for each linked vector */
100: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
101: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
102: while (next) {
103: next->rstart = nprev;
104: nprev += next->n;
105: next->grstart = com->rstart + next->rstart;
106: PetscMalloc1(size,&next->grstarts);
107: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
108: next = next->next;
109: }
110: com->setup = PETSC_TRUE;
111: return(0);
112: }
114: /* ----------------------------------------------------------------------------------*/
116: /*@
117: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
118: representation.
120: Not Collective
122: Input Parameter:
123: . dm - the packer object
125: Output Parameter:
126: . nDM - the number of DMs
128: Level: beginner
130: @*/
131: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
132: {
133: DM_Composite *com = (DM_Composite*)dm->data;
137: *nDM = com->nDM;
138: return(0);
139: }
142: /*@C
143: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
144: representation.
146: Collective on DMComposite
148: Input Parameters:
149: + dm - the packer object
150: - gvec - the global vector
152: Output Parameters:
153: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
155: Notes: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
157: Fortran Notes:
159: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
160: or use the alternative interface DMCompositeGetAccessArray().
162: Level: advanced
164: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
165: @*/
166: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
167: {
168: va_list Argp;
169: PetscErrorCode ierr;
170: struct DMCompositeLink *next;
171: DM_Composite *com = (DM_Composite*)dm->data;
172: PetscInt readonly;
177: next = com->next;
178: if (!com->setup) {
179: DMSetUp(dm);
180: }
182: VecLockGet(gvec,&readonly);
183: /* loop over packed objects, handling one at at time */
184: va_start(Argp,gvec);
185: while (next) {
186: Vec *vec;
187: vec = va_arg(Argp, Vec*);
188: if (vec) {
189: DMGetGlobalVector(next->dm,vec);
190: if (readonly) {
191: const PetscScalar *array;
192: VecGetArrayRead(gvec,&array);
193: VecPlaceArray(*vec,array+next->rstart);
194: VecLockPush(*vec);
195: VecRestoreArrayRead(gvec,&array);
196: } else {
197: PetscScalar *array;
198: VecGetArray(gvec,&array);
199: VecPlaceArray(*vec,array+next->rstart);
200: VecRestoreArray(gvec,&array);
201: }
202: }
203: next = next->next;
204: }
205: va_end(Argp);
206: return(0);
207: }
209: /*@C
210: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
211: representation.
213: Collective on DMComposite
215: Input Parameters:
216: + dm - the packer object
217: . pvec - packed vector
218: . nwanted - number of vectors wanted
219: - wanted - sorted array of vectors wanted, or NULL to get all vectors
221: Output Parameters:
222: . vecs - array of requested global vectors (must be allocated)
224: Notes: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
226: Level: advanced
228: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
229: @*/
230: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
231: {
232: PetscErrorCode ierr;
233: struct DMCompositeLink *link;
234: PetscInt i,wnum;
235: DM_Composite *com = (DM_Composite*)dm->data;
236: PetscInt readonly;
241: if (!com->setup) {
242: DMSetUp(dm);
243: }
245: VecLockGet(pvec,&readonly);
246: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
247: if (!wanted || i == wanted[wnum]) {
248: Vec v;
249: DMGetGlobalVector(link->dm,&v);
250: if (readonly) {
251: const PetscScalar *array;
252: VecGetArrayRead(pvec,&array);
253: VecPlaceArray(v,array+link->rstart);
254: VecLockPush(v);
255: VecRestoreArrayRead(pvec,&array);
256: } else {
257: PetscScalar *array;
258: VecGetArray(pvec,&array);
259: VecPlaceArray(v,array+link->rstart);
260: VecRestoreArray(pvec,&array);
261: }
262: vecs[wnum++] = v;
263: }
264: }
265: return(0);
266: }
268: /*@C
269: DMCompositeGetLocalAccessArray - Allows one to access the individual
270: packed vectors in their local representation.
272: Collective on DMComposite.
274: Input Parameters:
275: + dm - the packer object
276: . pvec - packed vector
277: . nwanted - number of vectors wanted
278: - wanted - sorted array of vectors wanted, or NULL to get all vectors
280: Output Parameters:
281: . vecs - array of requested local vectors (must be allocated)
283: Notes: Use DMCompositeRestoreLocalAccessArray() to return the vectors
284: when you no longer need them.
286: Level: advanced
288: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
289: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
290: @*/
291: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
292: {
293: PetscErrorCode ierr;
294: struct DMCompositeLink *link;
295: PetscInt i,wnum;
296: DM_Composite *com = (DM_Composite*)dm->data;
297: PetscInt readonly;
298: PetscInt nlocal = 0;
303: if (!com->setup) {
304: DMSetUp(dm);
305: }
307: VecLockGet(pvec,&readonly);
308: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
309: if (!wanted || i == wanted[wnum]) {
310: Vec v;
311: DMGetLocalVector(link->dm,&v);
312: if (readonly) {
313: const PetscScalar *array;
314: VecGetArrayRead(pvec,&array);
315: VecPlaceArray(v,array+nlocal);
316: VecLockPush(v);
317: VecRestoreArrayRead(pvec,&array);
318: } else {
319: PetscScalar *array;
320: VecGetArray(pvec,&array);
321: VecPlaceArray(v,array+nlocal);
322: VecRestoreArray(pvec,&array);
323: }
324: vecs[wnum++] = v;
325: }
327: nlocal += link->nlocal;
328: }
330: return(0);
331: }
333: /*@C
334: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
335: representation.
337: Collective on DMComposite
339: Input Parameters:
340: + dm - the packer object
341: . gvec - the global vector
342: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
344: Level: advanced
346: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
347: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
348: DMCompositeRestoreAccess(), DMCompositeGetAccess()
350: @*/
351: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
352: {
353: va_list Argp;
354: PetscErrorCode ierr;
355: struct DMCompositeLink *next;
356: DM_Composite *com = (DM_Composite*)dm->data;
357: PetscInt readonly;
362: next = com->next;
363: if (!com->setup) {
364: DMSetUp(dm);
365: }
367: VecLockGet(gvec,&readonly);
368: /* loop over packed objects, handling one at at time */
369: va_start(Argp,gvec);
370: while (next) {
371: Vec *vec;
372: vec = va_arg(Argp, Vec*);
373: if (vec) {
374: VecResetArray(*vec);
375: if (readonly) {
376: VecLockPop(*vec);
377: }
378: DMRestoreGlobalVector(next->dm,vec);
379: }
380: next = next->next;
381: }
382: va_end(Argp);
383: return(0);
384: }
386: /*@C
387: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
389: Collective on DMComposite
391: Input Parameters:
392: + dm - the packer object
393: . pvec - packed vector
394: . nwanted - number of vectors wanted
395: . wanted - sorted array of vectors wanted, or NULL to get all vectors
396: - vecs - array of global vectors to return
398: Level: advanced
400: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
401: @*/
402: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
403: {
404: PetscErrorCode ierr;
405: struct DMCompositeLink *link;
406: PetscInt i,wnum;
407: DM_Composite *com = (DM_Composite*)dm->data;
408: PetscInt readonly;
413: if (!com->setup) {
414: DMSetUp(dm);
415: }
417: VecLockGet(pvec,&readonly);
418: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
419: if (!wanted || i == wanted[wnum]) {
420: VecResetArray(vecs[wnum]);
421: if (readonly) {
422: VecLockPop(vecs[wnum]);
423: }
424: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
425: wnum++;
426: }
427: }
428: return(0);
429: }
431: /*@C
432: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
434: Collective on DMComposite.
436: Input Parameters:
437: + dm - the packer object
438: . pvec - packed vector
439: . nwanted - number of vectors wanted
440: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
441: - vecs - array of local vectors to return
443: Level: advanced
445: Notes:
446: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
447: otherwise the call will fail.
449: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
450: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
451: DMCompositeScatter(), DMCompositeGather()
452: @*/
453: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
454: {
455: PetscErrorCode ierr;
456: struct DMCompositeLink *link;
457: PetscInt i,wnum;
458: DM_Composite *com = (DM_Composite*)dm->data;
459: PetscInt readonly;
464: if (!com->setup) {
465: DMSetUp(dm);
466: }
468: VecLockGet(pvec,&readonly);
469: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
470: if (!wanted || i == wanted[wnum]) {
471: VecResetArray(vecs[wnum]);
472: if (readonly) {
473: VecLockPop(vecs[wnum]);
474: }
475: DMRestoreLocalVector(link->dm,&vecs[wnum]);
476: wnum++;
477: }
478: }
479: return(0);
480: }
482: /*@C
483: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
485: Collective on DMComposite
487: Input Parameters:
488: + dm - the packer object
489: . gvec - the global vector
490: - Vec ... - the individual sequential vectors, NULL for those that are not needed
492: Level: advanced
494: Notes:
495: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
496: accessible from Fortran.
498: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
499: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
500: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
501: DMCompositeScatterArray()
503: @*/
504: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
505: {
506: va_list Argp;
507: PetscErrorCode ierr;
508: struct DMCompositeLink *next;
509: PetscInt cnt;
510: DM_Composite *com = (DM_Composite*)dm->data;
515: if (!com->setup) {
516: DMSetUp(dm);
517: }
519: /* loop over packed objects, handling one at at time */
520: va_start(Argp,gvec);
521: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
522: Vec local;
523: local = va_arg(Argp, Vec);
524: if (local) {
525: Vec global;
526: const PetscScalar *array;
528: DMGetGlobalVector(next->dm,&global);
529: VecGetArrayRead(gvec,&array);
530: VecPlaceArray(global,array+next->rstart);
531: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
532: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
533: VecRestoreArrayRead(gvec,&array);
534: VecResetArray(global);
535: DMRestoreGlobalVector(next->dm,&global);
536: }
537: }
538: va_end(Argp);
539: return(0);
540: }
542: /*@
543: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
545: Collective on DMComposite
547: Input Parameters:
548: + dm - the packer object
549: . gvec - the global vector
550: . lvecs - array of local vectors, NULL for any that are not needed
552: Level: advanced
554: Note:
555: This is a non-variadic alternative to DMCompositeScatter()
557: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
558: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
559: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
561: @*/
562: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
563: {
564: PetscErrorCode ierr;
565: struct DMCompositeLink *next;
566: PetscInt i;
567: DM_Composite *com = (DM_Composite*)dm->data;
572: if (!com->setup) {
573: DMSetUp(dm);
574: }
576: /* loop over packed objects, handling one at at time */
577: for (i=0,next=com->next; next; next=next->next,i++) {
578: if (lvecs[i]) {
579: Vec global;
580: const PetscScalar *array;
582: DMGetGlobalVector(next->dm,&global);
583: VecGetArrayRead(gvec,&array);
584: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
585: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
586: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
587: VecRestoreArrayRead(gvec,&array);
588: VecResetArray(global);
589: DMRestoreGlobalVector(next->dm,&global);
590: }
591: }
592: return(0);
593: }
595: /*@C
596: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
598: Collective on DMComposite
600: Input Parameter:
601: + dm - the packer object
602: . gvec - the global vector
603: . imode - INSERT_VALUES or ADD_VALUES
604: - Vec ... - the individual sequential vectors, NULL for any that are not needed
606: Level: advanced
608: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
609: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
610: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
612: @*/
613: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
614: {
615: va_list Argp;
616: PetscErrorCode ierr;
617: struct DMCompositeLink *next;
618: DM_Composite *com = (DM_Composite*)dm->data;
619: PetscInt cnt;
624: if (!com->setup) {
625: DMSetUp(dm);
626: }
628: /* loop over packed objects, handling one at at time */
629: va_start(Argp,gvec);
630: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
631: Vec local;
632: local = va_arg(Argp, Vec);
633: if (local) {
634: PetscScalar *array;
635: Vec global;
637: DMGetGlobalVector(next->dm,&global);
638: VecGetArray(gvec,&array);
639: VecPlaceArray(global,array+next->rstart);
640: DMLocalToGlobalBegin(next->dm,local,imode,global);
641: DMLocalToGlobalEnd(next->dm,local,imode,global);
642: VecRestoreArray(gvec,&array);
643: VecResetArray(global);
644: DMRestoreGlobalVector(next->dm,&global);
645: }
646: }
647: va_end(Argp);
648: return(0);
649: }
651: /*@
652: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
654: Collective on DMComposite
656: Input Parameter:
657: + dm - the packer object
658: . gvec - the global vector
659: . imode - INSERT_VALUES or ADD_VALUES
660: - lvecs - the individual sequential vectors, NULL for any that are not needed
662: Level: advanced
664: Notes:
665: This is a non-variadic alternative to DMCompositeGather().
667: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
668: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
669: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
670: @*/
671: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
672: {
673: PetscErrorCode ierr;
674: struct DMCompositeLink *next;
675: DM_Composite *com = (DM_Composite*)dm->data;
676: PetscInt i;
681: if (!com->setup) {
682: DMSetUp(dm);
683: }
685: /* loop over packed objects, handling one at at time */
686: for (next=com->next,i=0; next; next=next->next,i++) {
687: if (lvecs[i]) {
688: PetscScalar *array;
689: Vec global;
691: DMGetGlobalVector(next->dm,&global);
692: VecGetArray(gvec,&array);
693: VecPlaceArray(global,array+next->rstart);
694: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
695: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
696: VecRestoreArray(gvec,&array);
697: VecResetArray(global);
698: DMRestoreGlobalVector(next->dm,&global);
699: }
700: }
701: return(0);
702: }
704: /*@C
705: DMCompositeAddDM - adds a DM vector to a DMComposite
707: Collective on DMComposite
709: Input Parameter:
710: + dmc - the DMComposite (packer) object
711: - dm - the DM object; if the DM is a DMDA you will need to cast it with a (DM)
713: Level: advanced
715: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
716: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
717: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
719: @*/
720: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
721: {
722: PetscErrorCode ierr;
723: PetscInt n,nlocal;
724: struct DMCompositeLink *mine,*next;
725: Vec global,local;
726: DM_Composite *com = (DM_Composite*)dmc->data;
731: next = com->next;
732: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
734: /* create new link */
735: PetscNew(&mine);
736: PetscObjectReference((PetscObject)dm);
737: DMGetGlobalVector(dm,&global);
738: VecGetLocalSize(global,&n);
739: DMRestoreGlobalVector(dm,&global);
740: DMGetLocalVector(dm,&local);
741: VecGetSize(local,&nlocal);
742: DMRestoreLocalVector(dm,&local);
744: mine->n = n;
745: mine->nlocal = nlocal;
746: mine->dm = dm;
747: mine->next = NULL;
748: com->n += n;
749: com->nghost += nlocal;
751: /* add to end of list */
752: if (!next) com->next = mine;
753: else {
754: while (next->next) next = next->next;
755: next->next = mine;
756: }
757: com->nDM++;
758: com->nmine++;
759: return(0);
760: }
762: #include <petscdraw.h>
763: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
764: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
765: {
766: DM dm;
767: PetscErrorCode ierr;
768: struct DMCompositeLink *next;
769: PetscBool isdraw;
770: DM_Composite *com;
773: VecGetDM(gvec, &dm);
774: if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
775: com = (DM_Composite*)dm->data;
776: next = com->next;
778: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
779: if (!isdraw) {
780: /* do I really want to call this? */
781: VecView_MPI(gvec,viewer);
782: } else {
783: PetscInt cnt = 0;
785: /* loop over packed objects, handling one at at time */
786: while (next) {
787: Vec vec;
788: const PetscScalar *array;
789: PetscInt bs;
791: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
792: DMGetGlobalVector(next->dm,&vec);
793: VecGetArrayRead(gvec,&array);
794: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
795: VecRestoreArrayRead(gvec,&array);
796: VecView(vec,viewer);
797: VecResetArray(vec);
798: VecGetBlockSize(vec,&bs);
799: DMRestoreGlobalVector(next->dm,&vec);
800: PetscViewerDrawBaseAdd(viewer,bs);
801: cnt += bs;
802: next = next->next;
803: }
804: PetscViewerDrawBaseAdd(viewer,-cnt);
805: }
806: return(0);
807: }
809: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
810: {
812: DM_Composite *com = (DM_Composite*)dm->data;
816: DMSetUp(dm);
817: VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
818: VecSetDM(*gvec, dm);
819: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
820: return(0);
821: }
823: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
824: {
826: DM_Composite *com = (DM_Composite*)dm->data;
830: if (!com->setup) {
831: DMSetUp(dm);
832: }
833: VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
834: VecSetDM(*lvec, dm);
835: return(0);
836: }
838: /*@C
839: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
841: Collective on DM
843: Input Parameter:
844: . dm - the packer object
846: Output Parameters:
847: . ltogs - the individual mappings for each packed vector. Note that this includes
848: all the ghost points that individual ghosted DMDA's may have.
850: Level: advanced
852: Notes:
853: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
855: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
856: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
857: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
859: @*/
860: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
861: {
862: PetscErrorCode ierr;
863: PetscInt i,*idx,n,cnt;
864: struct DMCompositeLink *next;
865: PetscMPIInt rank;
866: DM_Composite *com = (DM_Composite*)dm->data;
870: DMSetUp(dm);
871: PetscMalloc1(com->nDM,ltogs);
872: next = com->next;
873: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
875: /* loop over packed objects, handling one at at time */
876: cnt = 0;
877: while (next) {
878: ISLocalToGlobalMapping ltog;
879: PetscMPIInt size;
880: const PetscInt *suboff,*indices;
881: Vec global;
883: /* Get sub-DM global indices for each local dof */
884: DMGetLocalToGlobalMapping(next->dm,<og);
885: ISLocalToGlobalMappingGetSize(ltog,&n);
886: ISLocalToGlobalMappingGetIndices(ltog,&indices);
887: PetscMalloc1(n,&idx);
889: /* Get the offsets for the sub-DM global vector */
890: DMGetGlobalVector(next->dm,&global);
891: VecGetOwnershipRanges(global,&suboff);
892: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
894: /* Shift the sub-DM definition of the global space to the composite global space */
895: for (i=0; i<n; i++) {
896: PetscInt subi = indices[i],lo = 0,hi = size,t;
897: /* Binary search to find which rank owns subi */
898: while (hi-lo > 1) {
899: t = lo + (hi-lo)/2;
900: if (suboff[t] > subi) hi = t;
901: else lo = t;
902: }
903: idx[i] = subi - suboff[lo] + next->grstarts[lo];
904: }
905: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
906: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
907: DMRestoreGlobalVector(next->dm,&global);
908: next = next->next;
909: cnt++;
910: }
911: return(0);
912: }
914: /*@C
915: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
917: Not Collective
919: Input Arguments:
920: . dm - composite DM
922: Output Arguments:
923: . is - array of serial index sets for each each component of the DMComposite
925: Level: intermediate
927: Notes:
928: At present, a composite local vector does not normally exist. This function is used to provide index sets for
929: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
930: scatter to a composite local vector. The user should not typically need to know which is being done.
932: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
934: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
936: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
938: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
939: @*/
940: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
941: {
942: PetscErrorCode ierr;
943: DM_Composite *com = (DM_Composite*)dm->data;
944: struct DMCompositeLink *link;
945: PetscInt cnt,start;
950: PetscMalloc1(com->nmine,is);
951: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
952: PetscInt bs;
953: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
954: DMGetBlockSize(link->dm,&bs);
955: ISSetBlockSize((*is)[cnt],bs);
956: }
957: return(0);
958: }
960: /*@C
961: DMCompositeGetGlobalISs - Gets the index sets for each composed object
963: Collective on DMComposite
965: Input Parameter:
966: . dm - the packer object
968: Output Parameters:
969: . is - the array of index sets
971: Level: advanced
973: Notes:
974: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
976: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
978: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
979: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
980: indices.
982: Fortran Notes:
984: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
986: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
987: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
988: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
990: @*/
991: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
992: {
993: PetscErrorCode ierr;
994: PetscInt cnt = 0;
995: struct DMCompositeLink *next;
996: PetscMPIInt rank;
997: DM_Composite *com = (DM_Composite*)dm->data;
1001: PetscMalloc1(com->nDM,is);
1002: next = com->next;
1003: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1005: /* loop over packed objects, handling one at at time */
1006: while (next) {
1007: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1008: if (dm->prob) {
1009: MatNullSpace space;
1010: Mat pmat;
1011: PetscObject disc;
1012: PetscInt Nf;
1014: PetscDSGetNumFields(dm->prob, &Nf);
1015: if (cnt < Nf) {
1016: PetscDSGetDiscretization(dm->prob, cnt, &disc);
1017: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1018: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1019: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1020: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1021: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1022: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1023: }
1024: }
1025: cnt++;
1026: next = next->next;
1027: }
1028: return(0);
1029: }
1031: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1032: {
1033: PetscInt nDM;
1034: DM *dms;
1035: PetscInt i;
1039: DMCompositeGetNumberDM(dm, &nDM);
1040: if (numFields) *numFields = nDM;
1041: DMCompositeGetGlobalISs(dm, fields);
1042: if (fieldNames) {
1043: PetscMalloc1(nDM, &dms);
1044: PetscMalloc1(nDM, fieldNames);
1045: DMCompositeGetEntriesArray(dm, dms);
1046: for (i=0; i<nDM; i++) {
1047: char buf[256];
1048: const char *splitname;
1050: /* Split naming precedence: object name, prefix, number */
1051: splitname = ((PetscObject) dm)->name;
1052: if (!splitname) {
1053: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1054: if (splitname) {
1055: size_t len;
1056: PetscStrncpy(buf,splitname,sizeof(buf));
1057: buf[sizeof(buf) - 1] = 0;
1058: PetscStrlen(buf,&len);
1059: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1060: splitname = buf;
1061: }
1062: }
1063: if (!splitname) {
1064: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1065: splitname = buf;
1066: }
1067: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1068: }
1069: PetscFree(dms);
1070: }
1071: return(0);
1072: }
1074: /*
1075: This could take over from DMCreateFieldIS(), as it is more general,
1076: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1077: At this point it's probably best to be less intrusive, however.
1078: */
1079: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1080: {
1081: PetscInt nDM;
1082: PetscInt i;
1086: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1087: if (dmlist) {
1088: DMCompositeGetNumberDM(dm, &nDM);
1089: PetscMalloc1(nDM, dmlist);
1090: DMCompositeGetEntriesArray(dm, *dmlist);
1091: for (i=0; i<nDM; i++) {
1092: PetscObjectReference((PetscObject)((*dmlist)[i]));
1093: }
1094: }
1095: return(0);
1096: }
1100: /* -------------------------------------------------------------------------------------*/
1101: /*@C
1102: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1103: Use DMCompositeRestoreLocalVectors() to return them.
1105: Not Collective
1107: Input Parameter:
1108: . dm - the packer object
1110: Output Parameter:
1111: . Vec ... - the individual sequential Vecs
1113: Level: advanced
1115: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1116: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1117: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1119: @*/
1120: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1121: {
1122: va_list Argp;
1123: PetscErrorCode ierr;
1124: struct DMCompositeLink *next;
1125: DM_Composite *com = (DM_Composite*)dm->data;
1129: next = com->next;
1130: /* loop over packed objects, handling one at at time */
1131: va_start(Argp,dm);
1132: while (next) {
1133: Vec *vec;
1134: vec = va_arg(Argp, Vec*);
1135: if (vec) {DMGetLocalVector(next->dm,vec);}
1136: next = next->next;
1137: }
1138: va_end(Argp);
1139: return(0);
1140: }
1142: /*@C
1143: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1145: Not Collective
1147: Input Parameter:
1148: . dm - the packer object
1150: Output Parameter:
1151: . Vec ... - the individual sequential Vecs
1153: Level: advanced
1155: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1156: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1157: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1159: @*/
1160: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1161: {
1162: va_list Argp;
1163: PetscErrorCode ierr;
1164: struct DMCompositeLink *next;
1165: DM_Composite *com = (DM_Composite*)dm->data;
1169: next = com->next;
1170: /* loop over packed objects, handling one at at time */
1171: va_start(Argp,dm);
1172: while (next) {
1173: Vec *vec;
1174: vec = va_arg(Argp, Vec*);
1175: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1176: next = next->next;
1177: }
1178: va_end(Argp);
1179: return(0);
1180: }
1182: /* -------------------------------------------------------------------------------------*/
1183: /*@C
1184: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1186: Not Collective
1188: Input Parameter:
1189: . dm - the packer object
1191: Output Parameter:
1192: . DM ... - the individual entries (DMs)
1194: Level: advanced
1196: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1197: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1198: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1199: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1201: @*/
1202: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1203: {
1204: va_list Argp;
1205: struct DMCompositeLink *next;
1206: DM_Composite *com = (DM_Composite*)dm->data;
1210: next = com->next;
1211: /* loop over packed objects, handling one at at time */
1212: va_start(Argp,dm);
1213: while (next) {
1214: DM *dmn;
1215: dmn = va_arg(Argp, DM*);
1216: if (dmn) *dmn = next->dm;
1217: next = next->next;
1218: }
1219: va_end(Argp);
1220: return(0);
1221: }
1223: /*@C
1224: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1226: Not Collective
1228: Input Parameter:
1229: . dm - the packer object
1231: Output Parameter:
1232: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1234: Level: advanced
1236: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1237: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1238: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1239: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1241: @*/
1242: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1243: {
1244: struct DMCompositeLink *next;
1245: DM_Composite *com = (DM_Composite*)dm->data;
1246: PetscInt i;
1250: /* loop over packed objects, handling one at at time */
1251: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1252: return(0);
1253: }
1255: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1256: {
1257: PetscErrorCode ierr;
1258: struct DMCompositeLink *next;
1259: DM_Composite *com = (DM_Composite*)dmi->data;
1260: DM dm;
1264: if (comm == MPI_COMM_NULL) {
1265: PetscObjectGetComm((PetscObject)dmi,&comm);
1266: }
1267: DMSetUp(dmi);
1268: next = com->next;
1269: DMCompositeCreate(comm,fine);
1271: /* loop over packed objects, handling one at at time */
1272: while (next) {
1273: DMRefine(next->dm,comm,&dm);
1274: DMCompositeAddDM(*fine,dm);
1275: PetscObjectDereference((PetscObject)dm);
1276: next = next->next;
1277: }
1278: return(0);
1279: }
1281: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1282: {
1283: PetscErrorCode ierr;
1284: struct DMCompositeLink *next;
1285: DM_Composite *com = (DM_Composite*)dmi->data;
1286: DM dm;
1290: DMSetUp(dmi);
1291: if (comm == MPI_COMM_NULL) {
1292: PetscObjectGetComm((PetscObject)dmi,&comm);
1293: }
1294: next = com->next;
1295: DMCompositeCreate(comm,fine);
1297: /* loop over packed objects, handling one at at time */
1298: while (next) {
1299: DMCoarsen(next->dm,comm,&dm);
1300: DMCompositeAddDM(*fine,dm);
1301: PetscObjectDereference((PetscObject)dm);
1302: next = next->next;
1303: }
1304: return(0);
1305: }
1307: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1308: {
1309: PetscErrorCode ierr;
1310: PetscInt m,n,M,N,nDM,i;
1311: struct DMCompositeLink *nextc;
1312: struct DMCompositeLink *nextf;
1313: Vec gcoarse,gfine,*vecs;
1314: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1315: DM_Composite *comfine = (DM_Composite*)fine->data;
1316: Mat *mats;
1321: DMSetUp(coarse);
1322: DMSetUp(fine);
1323: /* use global vectors only for determining matrix layout */
1324: DMGetGlobalVector(coarse,&gcoarse);
1325: DMGetGlobalVector(fine,&gfine);
1326: VecGetLocalSize(gcoarse,&n);
1327: VecGetLocalSize(gfine,&m);
1328: VecGetSize(gcoarse,&N);
1329: VecGetSize(gfine,&M);
1330: DMRestoreGlobalVector(coarse,&gcoarse);
1331: DMRestoreGlobalVector(fine,&gfine);
1333: nDM = comfine->nDM;
1334: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1335: PetscCalloc1(nDM*nDM,&mats);
1336: if (v) {
1337: PetscCalloc1(nDM,&vecs);
1338: }
1340: /* loop over packed objects, handling one at at time */
1341: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1342: if (!v) {
1343: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1344: } else {
1345: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1346: }
1347: }
1348: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1349: if (v) {
1350: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1351: }
1352: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1353: PetscFree(mats);
1354: if (v) {
1355: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1356: PetscFree(vecs);
1357: }
1358: return(0);
1359: }
1361: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1362: {
1363: DM_Composite *com = (DM_Composite*)dm->data;
1364: ISLocalToGlobalMapping *ltogs;
1365: PetscInt i;
1366: PetscErrorCode ierr;
1369: /* Set the ISLocalToGlobalMapping on the new matrix */
1370: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1371: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1372: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1373: PetscFree(ltogs);
1374: return(0);
1375: }
1378: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1379: {
1380: PetscErrorCode ierr;
1381: PetscInt n,i,cnt;
1382: ISColoringValue *colors;
1383: PetscBool dense = PETSC_FALSE;
1384: ISColoringValue maxcol = 0;
1385: DM_Composite *com = (DM_Composite*)dm->data;
1389: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1390: else if (ctype == IS_COLORING_GLOBAL) {
1391: n = com->n;
1392: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1393: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1395: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1396: if (dense) {
1397: for (i=0; i<n; i++) {
1398: colors[i] = (ISColoringValue)(com->rstart + i);
1399: }
1400: maxcol = com->N;
1401: } else {
1402: struct DMCompositeLink *next = com->next;
1403: PetscMPIInt rank;
1405: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1406: cnt = 0;
1407: while (next) {
1408: ISColoring lcoloring;
1410: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1411: for (i=0; i<lcoloring->N; i++) {
1412: colors[cnt++] = maxcol + lcoloring->colors[i];
1413: }
1414: maxcol += lcoloring->n;
1415: ISColoringDestroy(&lcoloring);
1416: next = next->next;
1417: }
1418: }
1419: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1420: return(0);
1421: }
1423: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1424: {
1425: PetscErrorCode ierr;
1426: struct DMCompositeLink *next;
1427: PetscScalar *garray,*larray;
1428: DM_Composite *com = (DM_Composite*)dm->data;
1434: if (!com->setup) {
1435: DMSetUp(dm);
1436: }
1438: VecGetArray(gvec,&garray);
1439: VecGetArray(lvec,&larray);
1441: /* loop over packed objects, handling one at at time */
1442: next = com->next;
1443: while (next) {
1444: Vec local,global;
1445: PetscInt N;
1447: DMGetGlobalVector(next->dm,&global);
1448: VecGetLocalSize(global,&N);
1449: VecPlaceArray(global,garray);
1450: DMGetLocalVector(next->dm,&local);
1451: VecPlaceArray(local,larray);
1452: DMGlobalToLocalBegin(next->dm,global,mode,local);
1453: DMGlobalToLocalEnd(next->dm,global,mode,local);
1454: VecResetArray(global);
1455: VecResetArray(local);
1456: DMRestoreGlobalVector(next->dm,&global);
1457: DMRestoreLocalVector(next->dm,&local);
1459: larray += next->nlocal;
1460: garray += next->n;
1461: next = next->next;
1462: }
1464: VecRestoreArray(gvec,NULL);
1465: VecRestoreArray(lvec,NULL);
1466: return(0);
1467: }
1469: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1470: {
1475: return(0);
1476: }
1478: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1479: {
1480: PetscErrorCode ierr;
1481: struct DMCompositeLink *next;
1482: PetscScalar *larray,*garray;
1483: DM_Composite *com = (DM_Composite*)dm->data;
1490: if (!com->setup) {
1491: DMSetUp(dm);
1492: }
1494: VecGetArray(lvec,&larray);
1495: VecGetArray(gvec,&garray);
1497: /* loop over packed objects, handling one at at time */
1498: next = com->next;
1499: while (next) {
1500: Vec global,local;
1502: DMGetLocalVector(next->dm,&local);
1503: VecPlaceArray(local,larray);
1504: DMGetGlobalVector(next->dm,&global);
1505: VecPlaceArray(global,garray);
1506: DMLocalToGlobalBegin(next->dm,local,mode,global);
1507: DMLocalToGlobalEnd(next->dm,local,mode,global);
1508: VecResetArray(local);
1509: VecResetArray(global);
1510: DMRestoreGlobalVector(next->dm,&global);
1511: DMRestoreLocalVector(next->dm,&local);
1513: garray += next->n;
1514: larray += next->nlocal;
1515: next = next->next;
1516: }
1518: VecRestoreArray(gvec,NULL);
1519: VecRestoreArray(lvec,NULL);
1521: return(0);
1522: }
1524: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1525: {
1530: return(0);
1531: }
1533: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1534: {
1535: PetscErrorCode ierr;
1536: struct DMCompositeLink *next;
1537: PetscScalar *array1,*array2;
1538: DM_Composite *com = (DM_Composite*)dm->data;
1545: if (!com->setup) {
1546: DMSetUp(dm);
1547: }
1549: VecGetArray(vec1,&array1);
1550: VecGetArray(vec2,&array2);
1552: /* loop over packed objects, handling one at at time */
1553: next = com->next;
1554: while (next) {
1555: Vec local1,local2;
1557: DMGetLocalVector(next->dm,&local1);
1558: VecPlaceArray(local1,array1);
1559: DMGetLocalVector(next->dm,&local2);
1560: VecPlaceArray(local2,array2);
1561: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1562: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1563: VecResetArray(local2);
1564: DMRestoreLocalVector(next->dm,&local2);
1565: VecResetArray(local1);
1566: DMRestoreLocalVector(next->dm,&local1);
1568: array1 += next->nlocal;
1569: array2 += next->nlocal;
1570: next = next->next;
1571: }
1573: VecRestoreArray(vec1,NULL);
1574: VecRestoreArray(vec2,NULL);
1576: return(0);
1577: }
1579: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1580: {
1585: return(0);
1586: }
1588: /*MC
1589: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1591: Level: intermediate
1593: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1594: M*/
1597: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1598: {
1600: DM_Composite *com;
1603: PetscNewLog(p,&com);
1604: p->data = com;
1605: PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1606: com->n = 0;
1607: com->nghost = 0;
1608: com->next = NULL;
1609: com->nDM = 0;
1611: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1612: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1613: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1614: p->ops->createfieldis = DMCreateFieldIS_Composite;
1615: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1616: p->ops->refine = DMRefine_Composite;
1617: p->ops->coarsen = DMCoarsen_Composite;
1618: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1619: p->ops->creatematrix = DMCreateMatrix_Composite;
1620: p->ops->getcoloring = DMCreateColoring_Composite;
1621: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1622: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1623: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1624: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1625: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1626: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1627: p->ops->destroy = DMDestroy_Composite;
1628: p->ops->view = DMView_Composite;
1629: p->ops->setup = DMSetUp_Composite;
1630: return(0);
1631: }
1633: /*@C
1634: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1635: vectors made up of several subvectors.
1637: Collective on MPI_Comm
1639: Input Parameter:
1640: . comm - the processors that will share the global vector
1642: Output Parameters:
1643: . packer - the packer object
1645: Level: advanced
1647: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1648: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1649: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1651: @*/
1652: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1653: {
1658: DMCreate(comm,packer);
1659: DMSetType(*packer,DMCOMPOSITE);
1660: return(0);
1661: }