Actual source code: pack.c
petsc-3.12.5 2020-03-29
2: #include <../src/dm/impls/composite/packimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petsc/private/glvisviewerimpl.h>
5: #include <petscds.h>
7: /*@C
8: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
9: separate components (DMs) in a DMto build the correct matrix nonzero structure.
12: Logically Collective
14: Input Parameter:
15: + dm - the composite object
16: - formcouplelocations - routine to set the nonzero locations in the matrix
18: Level: advanced
20: Not available from Fortran
22: Notes:
23: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
24: this routine
26: @*/
27: PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt))
28: {
29: DM_Composite *com = (DM_Composite*)dm->data;
30: PetscBool flg;
34: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
35: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
36: com->FormCoupleLocations = FormCoupleLocations;
37: return(0);
38: }
40: PetscErrorCode DMDestroy_Composite(DM dm)
41: {
42: PetscErrorCode ierr;
43: struct DMCompositeLink *next, *prev;
44: DM_Composite *com = (DM_Composite*)dm->data;
47: next = com->next;
48: while (next) {
49: prev = next;
50: next = next->next;
51: DMDestroy(&prev->dm);
52: PetscFree(prev->grstarts);
53: PetscFree(prev);
54: }
55: PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);
56: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
57: PetscFree(com);
58: return(0);
59: }
61: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
62: {
64: PetscBool iascii;
65: DM_Composite *com = (DM_Composite*)dm->data;
68: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
69: if (iascii) {
70: struct DMCompositeLink *lnk = com->next;
71: PetscInt i;
73: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
74: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
75: PetscViewerASCIIPushTab(v);
76: for (i=0; lnk; lnk=lnk->next,i++) {
77: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
78: PetscViewerASCIIPushTab(v);
79: DMView(lnk->dm,v);
80: PetscViewerASCIIPopTab(v);
81: }
82: PetscViewerASCIIPopTab(v);
83: }
84: return(0);
85: }
87: /* --------------------------------------------------------------------------------------*/
88: PetscErrorCode DMSetUp_Composite(DM dm)
89: {
90: PetscErrorCode ierr;
91: PetscInt nprev = 0;
92: PetscMPIInt rank,size;
93: DM_Composite *com = (DM_Composite*)dm->data;
94: struct DMCompositeLink *next = com->next;
95: PetscLayout map;
98: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
99: PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
100: PetscLayoutSetLocalSize(map,com->n);
101: PetscLayoutSetSize(map,PETSC_DETERMINE);
102: PetscLayoutSetBlockSize(map,1);
103: PetscLayoutSetUp(map);
104: PetscLayoutGetSize(map,&com->N);
105: PetscLayoutGetRange(map,&com->rstart,NULL);
106: PetscLayoutDestroy(&map);
108: /* now set the rstart for each linked vector */
109: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
110: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
111: while (next) {
112: next->rstart = nprev;
113: nprev += next->n;
114: next->grstart = com->rstart + next->rstart;
115: PetscMalloc1(size,&next->grstarts);
116: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
117: next = next->next;
118: }
119: com->setup = PETSC_TRUE;
120: return(0);
121: }
123: /* ----------------------------------------------------------------------------------*/
125: /*@
126: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
127: representation.
129: Not Collective
131: Input Parameter:
132: . dm - the packer object
134: Output Parameter:
135: . nDM - the number of DMs
137: Level: beginner
139: @*/
140: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
141: {
142: DM_Composite *com = (DM_Composite*)dm->data;
143: PetscBool flg;
148: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
149: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
150: *nDM = com->nDM;
151: return(0);
152: }
155: /*@C
156: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
157: representation.
159: Collective on dm
161: Input Parameters:
162: + dm - the packer object
163: - gvec - the global vector
165: Output Parameters:
166: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
168: Notes:
169: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
171: Fortran Notes:
173: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
174: or use the alternative interface DMCompositeGetAccessArray().
176: Level: advanced
178: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
179: @*/
180: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
181: {
182: va_list Argp;
183: PetscErrorCode ierr;
184: struct DMCompositeLink *next;
185: DM_Composite *com = (DM_Composite*)dm->data;
186: PetscInt readonly;
187: PetscBool flg;
192: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
193: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
194: next = com->next;
195: if (!com->setup) {
196: DMSetUp(dm);
197: }
199: VecLockGet(gvec,&readonly);
200: /* loop over packed objects, handling one at at time */
201: va_start(Argp,gvec);
202: while (next) {
203: Vec *vec;
204: vec = va_arg(Argp, Vec*);
205: if (vec) {
206: DMGetGlobalVector(next->dm,vec);
207: if (readonly) {
208: const PetscScalar *array;
209: VecGetArrayRead(gvec,&array);
210: VecPlaceArray(*vec,array+next->rstart);
211: VecLockReadPush(*vec);
212: VecRestoreArrayRead(gvec,&array);
213: } else {
214: PetscScalar *array;
215: VecGetArray(gvec,&array);
216: VecPlaceArray(*vec,array+next->rstart);
217: VecRestoreArray(gvec,&array);
218: }
219: }
220: next = next->next;
221: }
222: va_end(Argp);
223: return(0);
224: }
226: /*@C
227: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
228: representation.
230: Collective on dm
232: Input Parameters:
233: + dm - the packer object
234: . pvec - packed vector
235: . nwanted - number of vectors wanted
236: - wanted - sorted array of vectors wanted, or NULL to get all vectors
238: Output Parameters:
239: . vecs - array of requested global vectors (must be allocated)
241: Notes:
242: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
244: Level: advanced
246: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
247: @*/
248: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
249: {
250: PetscErrorCode ierr;
251: struct DMCompositeLink *link;
252: PetscInt i,wnum;
253: DM_Composite *com = (DM_Composite*)dm->data;
254: PetscInt readonly;
255: PetscBool flg;
260: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
261: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
262: if (!com->setup) {
263: DMSetUp(dm);
264: }
266: VecLockGet(pvec,&readonly);
267: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
268: if (!wanted || i == wanted[wnum]) {
269: Vec v;
270: DMGetGlobalVector(link->dm,&v);
271: if (readonly) {
272: const PetscScalar *array;
273: VecGetArrayRead(pvec,&array);
274: VecPlaceArray(v,array+link->rstart);
275: VecLockReadPush(v);
276: VecRestoreArrayRead(pvec,&array);
277: } else {
278: PetscScalar *array;
279: VecGetArray(pvec,&array);
280: VecPlaceArray(v,array+link->rstart);
281: VecRestoreArray(pvec,&array);
282: }
283: vecs[wnum++] = v;
284: }
285: }
286: return(0);
287: }
289: /*@C
290: DMCompositeGetLocalAccessArray - Allows one to access the individual
291: packed vectors in their local representation.
293: Collective on dm.
295: Input Parameters:
296: + dm - the packer object
297: . pvec - packed vector
298: . nwanted - number of vectors wanted
299: - wanted - sorted array of vectors wanted, or NULL to get all vectors
301: Output Parameters:
302: . vecs - array of requested local vectors (must be allocated)
304: Notes:
305: Use DMCompositeRestoreLocalAccessArray() to return the vectors
306: when you no longer need them.
308: Level: advanced
310: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
311: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
312: @*/
313: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
314: {
315: PetscErrorCode ierr;
316: struct DMCompositeLink *link;
317: PetscInt i,wnum;
318: DM_Composite *com = (DM_Composite*)dm->data;
319: PetscInt readonly;
320: PetscInt nlocal = 0;
321: PetscBool flg;
326: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
327: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
328: if (!com->setup) {
329: DMSetUp(dm);
330: }
332: VecLockGet(pvec,&readonly);
333: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
334: if (!wanted || i == wanted[wnum]) {
335: Vec v;
336: DMGetLocalVector(link->dm,&v);
337: if (readonly) {
338: const PetscScalar *array;
339: VecGetArrayRead(pvec,&array);
340: VecPlaceArray(v,array+nlocal);
341: VecLockReadPush(v);
342: VecRestoreArrayRead(pvec,&array);
343: } else {
344: PetscScalar *array;
345: VecGetArray(pvec,&array);
346: VecPlaceArray(v,array+nlocal);
347: VecRestoreArray(pvec,&array);
348: }
349: vecs[wnum++] = v;
350: }
352: nlocal += link->nlocal;
353: }
355: return(0);
356: }
358: /*@C
359: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
360: representation.
362: Collective on dm
364: Input Parameters:
365: + dm - the packer object
366: . gvec - the global vector
367: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
369: Level: advanced
371: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
372: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
373: DMCompositeRestoreAccess(), DMCompositeGetAccess()
375: @*/
376: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
377: {
378: va_list Argp;
379: PetscErrorCode ierr;
380: struct DMCompositeLink *next;
381: DM_Composite *com = (DM_Composite*)dm->data;
382: PetscInt readonly;
383: PetscBool flg;
388: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
389: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
390: next = com->next;
391: if (!com->setup) {
392: DMSetUp(dm);
393: }
395: VecLockGet(gvec,&readonly);
396: /* loop over packed objects, handling one at at time */
397: va_start(Argp,gvec);
398: while (next) {
399: Vec *vec;
400: vec = va_arg(Argp, Vec*);
401: if (vec) {
402: VecResetArray(*vec);
403: if (readonly) {
404: VecLockReadPop(*vec);
405: }
406: DMRestoreGlobalVector(next->dm,vec);
407: }
408: next = next->next;
409: }
410: va_end(Argp);
411: return(0);
412: }
414: /*@C
415: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
417: Collective on dm
419: Input Parameters:
420: + dm - the packer object
421: . pvec - packed vector
422: . nwanted - number of vectors wanted
423: . wanted - sorted array of vectors wanted, or NULL to get all vectors
424: - vecs - array of global vectors to return
426: Level: advanced
428: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
429: @*/
430: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
431: {
432: PetscErrorCode ierr;
433: struct DMCompositeLink *link;
434: PetscInt i,wnum;
435: DM_Composite *com = (DM_Composite*)dm->data;
436: PetscInt readonly;
437: PetscBool flg;
442: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
443: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
444: if (!com->setup) {
445: DMSetUp(dm);
446: }
448: VecLockGet(pvec,&readonly);
449: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
450: if (!wanted || i == wanted[wnum]) {
451: VecResetArray(vecs[wnum]);
452: if (readonly) {
453: VecLockReadPop(vecs[wnum]);
454: }
455: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
456: wnum++;
457: }
458: }
459: return(0);
460: }
462: /*@C
463: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
465: Collective on dm.
467: Input Parameters:
468: + dm - the packer object
469: . pvec - packed vector
470: . nwanted - number of vectors wanted
471: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
472: - vecs - array of local vectors to return
474: Level: advanced
476: Notes:
477: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
478: otherwise the call will fail.
480: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
481: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
482: DMCompositeScatter(), DMCompositeGather()
483: @*/
484: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
485: {
486: PetscErrorCode ierr;
487: struct DMCompositeLink *link;
488: PetscInt i,wnum;
489: DM_Composite *com = (DM_Composite*)dm->data;
490: PetscInt readonly;
491: PetscBool flg;
496: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
497: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
498: if (!com->setup) {
499: DMSetUp(dm);
500: }
502: VecLockGet(pvec,&readonly);
503: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
504: if (!wanted || i == wanted[wnum]) {
505: VecResetArray(vecs[wnum]);
506: if (readonly) {
507: VecLockReadPop(vecs[wnum]);
508: }
509: DMRestoreLocalVector(link->dm,&vecs[wnum]);
510: wnum++;
511: }
512: }
513: return(0);
514: }
516: /*@C
517: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
519: Collective on dm
521: Input Parameters:
522: + dm - the packer object
523: . gvec - the global vector
524: - Vec ... - the individual sequential vectors, NULL for those that are not needed
526: Level: advanced
528: Notes:
529: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
530: accessible from Fortran.
532: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
533: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
534: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
535: DMCompositeScatterArray()
537: @*/
538: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
539: {
540: va_list Argp;
541: PetscErrorCode ierr;
542: struct DMCompositeLink *next;
543: PetscInt cnt;
544: DM_Composite *com = (DM_Composite*)dm->data;
545: PetscBool flg;
550: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
551: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
552: if (!com->setup) {
553: DMSetUp(dm);
554: }
556: /* loop over packed objects, handling one at at time */
557: va_start(Argp,gvec);
558: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
559: Vec local;
560: local = va_arg(Argp, Vec);
561: if (local) {
562: Vec global;
563: const PetscScalar *array;
565: DMGetGlobalVector(next->dm,&global);
566: VecGetArrayRead(gvec,&array);
567: VecPlaceArray(global,array+next->rstart);
568: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
569: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
570: VecRestoreArrayRead(gvec,&array);
571: VecResetArray(global);
572: DMRestoreGlobalVector(next->dm,&global);
573: }
574: }
575: va_end(Argp);
576: return(0);
577: }
579: /*@
580: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
582: Collective on dm
584: Input Parameters:
585: + dm - the packer object
586: . gvec - the global vector
587: - lvecs - array of local vectors, NULL for any that are not needed
589: Level: advanced
591: Note:
592: This is a non-variadic alternative to DMCompositeScatter()
594: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
595: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
596: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
598: @*/
599: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
600: {
601: PetscErrorCode ierr;
602: struct DMCompositeLink *next;
603: PetscInt i;
604: DM_Composite *com = (DM_Composite*)dm->data;
605: PetscBool flg;
610: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
611: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
612: if (!com->setup) {
613: DMSetUp(dm);
614: }
616: /* loop over packed objects, handling one at at time */
617: for (i=0,next=com->next; next; next=next->next,i++) {
618: if (lvecs[i]) {
619: Vec global;
620: const PetscScalar *array;
622: DMGetGlobalVector(next->dm,&global);
623: VecGetArrayRead(gvec,&array);
624: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
625: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
626: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
627: VecRestoreArrayRead(gvec,&array);
628: VecResetArray(global);
629: DMRestoreGlobalVector(next->dm,&global);
630: }
631: }
632: return(0);
633: }
635: /*@C
636: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
638: Collective on dm
640: Input Parameter:
641: + dm - the packer object
642: . gvec - the global vector
643: . imode - INSERT_VALUES or ADD_VALUES
644: - Vec ... - the individual sequential vectors, NULL for any that are not needed
646: Level: advanced
648: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
650: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
651: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
652: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
654: @*/
655: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
656: {
657: va_list Argp;
658: PetscErrorCode ierr;
659: struct DMCompositeLink *next;
660: DM_Composite *com = (DM_Composite*)dm->data;
661: PetscInt cnt;
662: PetscBool flg;
667: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
668: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
669: if (!com->setup) {
670: DMSetUp(dm);
671: }
673: /* loop over packed objects, handling one at at time */
674: va_start(Argp,gvec);
675: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
676: Vec local;
677: local = va_arg(Argp, Vec);
678: if (local) {
679: PetscScalar *array;
680: Vec global;
682: DMGetGlobalVector(next->dm,&global);
683: VecGetArray(gvec,&array);
684: VecPlaceArray(global,array+next->rstart);
685: DMLocalToGlobalBegin(next->dm,local,imode,global);
686: DMLocalToGlobalEnd(next->dm,local,imode,global);
687: VecRestoreArray(gvec,&array);
688: VecResetArray(global);
689: DMRestoreGlobalVector(next->dm,&global);
690: }
691: }
692: va_end(Argp);
693: return(0);
694: }
696: /*@
697: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
699: Collective on dm
701: Input Parameter:
702: + dm - the packer object
703: . gvec - the global vector
704: . imode - INSERT_VALUES or ADD_VALUES
705: - lvecs - the individual sequential vectors, NULL for any that are not needed
707: Level: advanced
709: Notes:
710: This is a non-variadic alternative to DMCompositeGather().
712: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
713: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
714: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
715: @*/
716: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
717: {
718: PetscErrorCode ierr;
719: struct DMCompositeLink *next;
720: DM_Composite *com = (DM_Composite*)dm->data;
721: PetscInt i;
722: PetscBool flg;
727: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
728: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
729: if (!com->setup) {
730: DMSetUp(dm);
731: }
733: /* loop over packed objects, handling one at at time */
734: for (next=com->next,i=0; next; next=next->next,i++) {
735: if (lvecs[i]) {
736: PetscScalar *array;
737: Vec global;
739: DMGetGlobalVector(next->dm,&global);
740: VecGetArray(gvec,&array);
741: VecPlaceArray(global,array+next->rstart);
742: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
743: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
744: VecRestoreArray(gvec,&array);
745: VecResetArray(global);
746: DMRestoreGlobalVector(next->dm,&global);
747: }
748: }
749: return(0);
750: }
752: /*@
753: DMCompositeAddDM - adds a DM vector to a DMComposite
755: Collective on dm
757: Input Parameter:
758: + dmc - the DMComposite (packer) object
759: - dm - the DM object
761: Level: advanced
763: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
764: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
765: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
767: @*/
768: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
769: {
770: PetscErrorCode ierr;
771: PetscInt n,nlocal;
772: struct DMCompositeLink *mine,*next;
773: Vec global,local;
774: DM_Composite *com = (DM_Composite*)dmc->data;
775: PetscBool flg;
780: PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
781: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
782: next = com->next;
783: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
785: /* create new link */
786: PetscNew(&mine);
787: PetscObjectReference((PetscObject)dm);
788: DMGetGlobalVector(dm,&global);
789: VecGetLocalSize(global,&n);
790: DMRestoreGlobalVector(dm,&global);
791: DMGetLocalVector(dm,&local);
792: VecGetSize(local,&nlocal);
793: DMRestoreLocalVector(dm,&local);
795: mine->n = n;
796: mine->nlocal = nlocal;
797: mine->dm = dm;
798: mine->next = NULL;
799: com->n += n;
800: com->nghost += nlocal;
802: /* add to end of list */
803: if (!next) com->next = mine;
804: else {
805: while (next->next) next = next->next;
806: next->next = mine;
807: }
808: com->nDM++;
809: com->nmine++;
810: return(0);
811: }
813: #include <petscdraw.h>
814: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
815: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
816: {
817: DM dm;
818: PetscErrorCode ierr;
819: struct DMCompositeLink *next;
820: PetscBool isdraw;
821: DM_Composite *com;
824: VecGetDM(gvec, &dm);
825: if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
826: com = (DM_Composite*)dm->data;
827: next = com->next;
829: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
830: if (!isdraw) {
831: /* do I really want to call this? */
832: VecView_MPI(gvec,viewer);
833: } else {
834: PetscInt cnt = 0;
836: /* loop over packed objects, handling one at at time */
837: while (next) {
838: Vec vec;
839: const PetscScalar *array;
840: PetscInt bs;
842: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
843: DMGetGlobalVector(next->dm,&vec);
844: VecGetArrayRead(gvec,&array);
845: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
846: VecRestoreArrayRead(gvec,&array);
847: VecView(vec,viewer);
848: VecResetArray(vec);
849: VecGetBlockSize(vec,&bs);
850: DMRestoreGlobalVector(next->dm,&vec);
851: PetscViewerDrawBaseAdd(viewer,bs);
852: cnt += bs;
853: next = next->next;
854: }
855: PetscViewerDrawBaseAdd(viewer,-cnt);
856: }
857: return(0);
858: }
860: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
861: {
863: DM_Composite *com = (DM_Composite*)dm->data;
867: DMSetUp(dm);
868: VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
869: VecSetDM(*gvec, dm);
870: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
871: return(0);
872: }
874: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
875: {
877: DM_Composite *com = (DM_Composite*)dm->data;
881: if (!com->setup) {
882: DMSetUp(dm);
883: }
884: VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
885: VecSetDM(*lvec, dm);
886: return(0);
887: }
889: /*@C
890: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
892: Collective on DM
894: Input Parameter:
895: . dm - the packer object
897: Output Parameters:
898: . ltogs - the individual mappings for each packed vector. Note that this includes
899: all the ghost points that individual ghosted DMDA's may have.
901: Level: advanced
903: Notes:
904: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
906: Not available from Fortran
908: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
909: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
910: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
912: @*/
913: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
914: {
915: PetscErrorCode ierr;
916: PetscInt i,*idx,n,cnt;
917: struct DMCompositeLink *next;
918: PetscMPIInt rank;
919: DM_Composite *com = (DM_Composite*)dm->data;
920: PetscBool flg;
924: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
925: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
926: DMSetUp(dm);
927: PetscMalloc1(com->nDM,ltogs);
928: next = com->next;
929: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
931: /* loop over packed objects, handling one at at time */
932: cnt = 0;
933: while (next) {
934: ISLocalToGlobalMapping ltog;
935: PetscMPIInt size;
936: const PetscInt *suboff,*indices;
937: Vec global;
939: /* Get sub-DM global indices for each local dof */
940: DMGetLocalToGlobalMapping(next->dm,<og);
941: ISLocalToGlobalMappingGetSize(ltog,&n);
942: ISLocalToGlobalMappingGetIndices(ltog,&indices);
943: PetscMalloc1(n,&idx);
945: /* Get the offsets for the sub-DM global vector */
946: DMGetGlobalVector(next->dm,&global);
947: VecGetOwnershipRanges(global,&suboff);
948: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
950: /* Shift the sub-DM definition of the global space to the composite global space */
951: for (i=0; i<n; i++) {
952: PetscInt subi = indices[i],lo = 0,hi = size,t;
953: /* There's no consensus on what a negative index means,
954: except for skipping when setting the values in vectors and matrices */
955: if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
956: /* Binary search to find which rank owns subi */
957: while (hi-lo > 1) {
958: t = lo + (hi-lo)/2;
959: if (suboff[t] > subi) hi = t;
960: else lo = t;
961: }
962: idx[i] = subi - suboff[lo] + next->grstarts[lo];
963: }
964: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
965: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
966: DMRestoreGlobalVector(next->dm,&global);
967: next = next->next;
968: cnt++;
969: }
970: return(0);
971: }
973: /*@C
974: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
976: Not Collective
978: Input Arguments:
979: . dm - composite DM
981: Output Arguments:
982: . is - array of serial index sets for each each component of the DMComposite
984: Level: intermediate
986: Notes:
987: At present, a composite local vector does not normally exist. This function is used to provide index sets for
988: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
989: scatter to a composite local vector. The user should not typically need to know which is being done.
991: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
993: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
995: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
997: Not available from Fortran
999: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
1000: @*/
1001: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
1002: {
1003: PetscErrorCode ierr;
1004: DM_Composite *com = (DM_Composite*)dm->data;
1005: struct DMCompositeLink *link;
1006: PetscInt cnt,start;
1007: PetscBool flg;
1012: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1013: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1014: PetscMalloc1(com->nmine,is);
1015: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1016: PetscInt bs;
1017: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1018: DMGetBlockSize(link->dm,&bs);
1019: ISSetBlockSize((*is)[cnt],bs);
1020: }
1021: return(0);
1022: }
1024: /*@C
1025: DMCompositeGetGlobalISs - Gets the index sets for each composed object
1027: Collective on dm
1029: Input Parameter:
1030: . dm - the packer object
1032: Output Parameters:
1033: . is - the array of index sets
1035: Level: advanced
1037: Notes:
1038: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1040: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1042: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1043: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1044: indices.
1046: Fortran Notes:
1048: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1050: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1051: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1052: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1054: @*/
1055: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1056: {
1057: PetscErrorCode ierr;
1058: PetscInt cnt = 0;
1059: struct DMCompositeLink *next;
1060: PetscMPIInt rank;
1061: DM_Composite *com = (DM_Composite*)dm->data;
1062: PetscBool flg;
1066: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1067: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1068: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1069: PetscMalloc1(com->nDM,is);
1070: next = com->next;
1071: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1073: /* loop over packed objects, handling one at at time */
1074: while (next) {
1075: PetscDS prob;
1077: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1078: DMGetDS(dm, &prob);
1079: if (prob) {
1080: MatNullSpace space;
1081: Mat pmat;
1082: PetscObject disc;
1083: PetscInt Nf;
1085: PetscDSGetNumFields(prob, &Nf);
1086: if (cnt < Nf) {
1087: PetscDSGetDiscretization(prob, cnt, &disc);
1088: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1089: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1090: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1091: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1092: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1093: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1094: }
1095: }
1096: cnt++;
1097: next = next->next;
1098: }
1099: return(0);
1100: }
1102: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1103: {
1104: PetscInt nDM;
1105: DM *dms;
1106: PetscInt i;
1110: DMCompositeGetNumberDM(dm, &nDM);
1111: if (numFields) *numFields = nDM;
1112: DMCompositeGetGlobalISs(dm, fields);
1113: if (fieldNames) {
1114: PetscMalloc1(nDM, &dms);
1115: PetscMalloc1(nDM, fieldNames);
1116: DMCompositeGetEntriesArray(dm, dms);
1117: for (i=0; i<nDM; i++) {
1118: char buf[256];
1119: const char *splitname;
1121: /* Split naming precedence: object name, prefix, number */
1122: splitname = ((PetscObject) dm)->name;
1123: if (!splitname) {
1124: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1125: if (splitname) {
1126: size_t len;
1127: PetscStrncpy(buf,splitname,sizeof(buf));
1128: buf[sizeof(buf) - 1] = 0;
1129: PetscStrlen(buf,&len);
1130: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1131: splitname = buf;
1132: }
1133: }
1134: if (!splitname) {
1135: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1136: splitname = buf;
1137: }
1138: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1139: }
1140: PetscFree(dms);
1141: }
1142: return(0);
1143: }
1145: /*
1146: This could take over from DMCreateFieldIS(), as it is more general,
1147: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1148: At this point it's probably best to be less intrusive, however.
1149: */
1150: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1151: {
1152: PetscInt nDM;
1153: PetscInt i;
1157: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1158: if (dmlist) {
1159: DMCompositeGetNumberDM(dm, &nDM);
1160: PetscMalloc1(nDM, dmlist);
1161: DMCompositeGetEntriesArray(dm, *dmlist);
1162: for (i=0; i<nDM; i++) {
1163: PetscObjectReference((PetscObject)((*dmlist)[i]));
1164: }
1165: }
1166: return(0);
1167: }
1171: /* -------------------------------------------------------------------------------------*/
1172: /*@C
1173: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1174: Use DMCompositeRestoreLocalVectors() to return them.
1176: Not Collective
1178: Input Parameter:
1179: . dm - the packer object
1181: Output Parameter:
1182: . Vec ... - the individual sequential Vecs
1184: Level: advanced
1186: Not available from Fortran
1188: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1189: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1190: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1192: @*/
1193: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1194: {
1195: va_list Argp;
1196: PetscErrorCode ierr;
1197: struct DMCompositeLink *next;
1198: DM_Composite *com = (DM_Composite*)dm->data;
1199: PetscBool flg;
1203: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1204: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1205: next = com->next;
1206: /* loop over packed objects, handling one at at time */
1207: va_start(Argp,dm);
1208: while (next) {
1209: Vec *vec;
1210: vec = va_arg(Argp, Vec*);
1211: if (vec) {DMGetLocalVector(next->dm,vec);}
1212: next = next->next;
1213: }
1214: va_end(Argp);
1215: return(0);
1216: }
1218: /*@C
1219: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1221: Not Collective
1223: Input Parameter:
1224: . dm - the packer object
1226: Output Parameter:
1227: . Vec ... - the individual sequential Vecs
1229: Level: advanced
1231: Not available from Fortran
1233: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1234: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1235: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1237: @*/
1238: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1239: {
1240: va_list Argp;
1241: PetscErrorCode ierr;
1242: struct DMCompositeLink *next;
1243: DM_Composite *com = (DM_Composite*)dm->data;
1244: PetscBool flg;
1248: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1249: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1250: next = com->next;
1251: /* loop over packed objects, handling one at at time */
1252: va_start(Argp,dm);
1253: while (next) {
1254: Vec *vec;
1255: vec = va_arg(Argp, Vec*);
1256: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1257: next = next->next;
1258: }
1259: va_end(Argp);
1260: return(0);
1261: }
1263: /* -------------------------------------------------------------------------------------*/
1264: /*@C
1265: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1267: Not Collective
1269: Input Parameter:
1270: . dm - the packer object
1272: Output Parameter:
1273: . DM ... - the individual entries (DMs)
1275: Level: advanced
1277: Fortran Notes:
1278: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1280: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1281: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1282: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1283: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1285: @*/
1286: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1287: {
1288: va_list Argp;
1289: struct DMCompositeLink *next;
1290: DM_Composite *com = (DM_Composite*)dm->data;
1291: PetscBool flg;
1292: PetscErrorCode ierr;
1296: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1297: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1298: next = com->next;
1299: /* loop over packed objects, handling one at at time */
1300: va_start(Argp,dm);
1301: while (next) {
1302: DM *dmn;
1303: dmn = va_arg(Argp, DM*);
1304: if (dmn) *dmn = next->dm;
1305: next = next->next;
1306: }
1307: va_end(Argp);
1308: return(0);
1309: }
1311: /*@C
1312: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1314: Not Collective
1316: Input Parameter:
1317: . dm - the packer object
1319: Output Parameter:
1320: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1322: Level: advanced
1324: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1325: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1326: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1327: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1329: @*/
1330: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1331: {
1332: struct DMCompositeLink *next;
1333: DM_Composite *com = (DM_Composite*)dm->data;
1334: PetscInt i;
1335: PetscBool flg;
1336: PetscErrorCode ierr;
1340: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1341: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1342: /* loop over packed objects, handling one at at time */
1343: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1344: return(0);
1345: }
1347: typedef struct {
1348: DM dm;
1349: PetscViewer *subv;
1350: Vec *vecs;
1351: } GLVisViewerCtx;
1353: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1354: {
1355: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1356: PetscInt i,n;
1360: DMCompositeGetNumberDM(ctx->dm,&n);
1361: for (i = 0; i < n; i++) {
1362: PetscViewerDestroy(&ctx->subv[i]);
1363: }
1364: PetscFree2(ctx->subv,ctx->vecs);
1365: DMDestroy(&ctx->dm);
1366: PetscFree(ctx);
1367: return(0);
1368: }
1370: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1371: {
1372: Vec X = (Vec)oX;
1373: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1374: PetscInt i,n,cumf;
1378: DMCompositeGetNumberDM(ctx->dm,&n);
1379: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1380: for (i = 0, cumf = 0; i < n; i++) {
1381: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1382: void *fctx;
1383: PetscInt nfi;
1385: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1386: if (!nfi) continue;
1387: if (g2l) {
1388: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1389: } else {
1390: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1391: }
1392: cumf += nfi;
1393: }
1394: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1395: return(0);
1396: }
1398: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1399: {
1400: DM dm = (DM)odm, *dms;
1401: Vec *Ufds;
1402: GLVisViewerCtx *ctx;
1403: PetscInt i,n,tnf,*sdim;
1404: char **fecs;
1408: PetscNew(&ctx);
1409: PetscObjectReference((PetscObject)dm);
1410: ctx->dm = dm;
1411: DMCompositeGetNumberDM(dm,&n);
1412: PetscMalloc1(n,&dms);
1413: DMCompositeGetEntriesArray(dm,dms);
1414: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1415: for (i = 0, tnf = 0; i < n; i++) {
1416: PetscInt nf;
1418: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1419: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1420: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1421: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1422: tnf += nf;
1423: }
1424: PetscFree(dms);
1425: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1426: for (i = 0, tnf = 0; i < n; i++) {
1427: PetscInt *sd,nf,f;
1428: const char **fec;
1429: Vec *Uf;
1431: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1432: for (f = 0; f < nf; f++) {
1433: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1434: Ufds[tnf+f] = Uf[f];
1435: sdim[tnf+f] = sd[f];
1436: }
1437: tnf += nf;
1438: }
1439: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1440: for (i = 0; i < tnf; i++) {
1441: PetscFree(fecs[i]);
1442: }
1443: PetscFree3(fecs,sdim,Ufds);
1444: return(0);
1445: }
1447: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1448: {
1449: PetscErrorCode ierr;
1450: struct DMCompositeLink *next;
1451: DM_Composite *com = (DM_Composite*)dmi->data;
1452: DM dm;
1456: if (comm == MPI_COMM_NULL) {
1457: PetscObjectGetComm((PetscObject)dmi,&comm);
1458: }
1459: DMSetUp(dmi);
1460: next = com->next;
1461: DMCompositeCreate(comm,fine);
1463: /* loop over packed objects, handling one at at time */
1464: while (next) {
1465: DMRefine(next->dm,comm,&dm);
1466: DMCompositeAddDM(*fine,dm);
1467: PetscObjectDereference((PetscObject)dm);
1468: next = next->next;
1469: }
1470: return(0);
1471: }
1473: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1474: {
1475: PetscErrorCode ierr;
1476: struct DMCompositeLink *next;
1477: DM_Composite *com = (DM_Composite*)dmi->data;
1478: DM dm;
1482: DMSetUp(dmi);
1483: if (comm == MPI_COMM_NULL) {
1484: PetscObjectGetComm((PetscObject)dmi,&comm);
1485: }
1486: next = com->next;
1487: DMCompositeCreate(comm,fine);
1489: /* loop over packed objects, handling one at at time */
1490: while (next) {
1491: DMCoarsen(next->dm,comm,&dm);
1492: DMCompositeAddDM(*fine,dm);
1493: PetscObjectDereference((PetscObject)dm);
1494: next = next->next;
1495: }
1496: return(0);
1497: }
1499: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1500: {
1501: PetscErrorCode ierr;
1502: PetscInt m,n,M,N,nDM,i;
1503: struct DMCompositeLink *nextc;
1504: struct DMCompositeLink *nextf;
1505: Vec gcoarse,gfine,*vecs;
1506: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1507: DM_Composite *comfine = (DM_Composite*)fine->data;
1508: Mat *mats;
1513: DMSetUp(coarse);
1514: DMSetUp(fine);
1515: /* use global vectors only for determining matrix layout */
1516: DMGetGlobalVector(coarse,&gcoarse);
1517: DMGetGlobalVector(fine,&gfine);
1518: VecGetLocalSize(gcoarse,&n);
1519: VecGetLocalSize(gfine,&m);
1520: VecGetSize(gcoarse,&N);
1521: VecGetSize(gfine,&M);
1522: DMRestoreGlobalVector(coarse,&gcoarse);
1523: DMRestoreGlobalVector(fine,&gfine);
1525: nDM = comfine->nDM;
1526: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1527: PetscCalloc1(nDM*nDM,&mats);
1528: if (v) {
1529: PetscCalloc1(nDM,&vecs);
1530: }
1532: /* loop over packed objects, handling one at at time */
1533: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1534: if (!v) {
1535: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1536: } else {
1537: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1538: }
1539: }
1540: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1541: if (v) {
1542: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1543: }
1544: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1545: PetscFree(mats);
1546: if (v) {
1547: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1548: PetscFree(vecs);
1549: }
1550: return(0);
1551: }
1553: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1554: {
1555: DM_Composite *com = (DM_Composite*)dm->data;
1556: ISLocalToGlobalMapping *ltogs;
1557: PetscInt i;
1558: PetscErrorCode ierr;
1561: /* Set the ISLocalToGlobalMapping on the new matrix */
1562: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1563: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1564: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1565: PetscFree(ltogs);
1566: return(0);
1567: }
1570: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1571: {
1572: PetscErrorCode ierr;
1573: PetscInt n,i,cnt;
1574: ISColoringValue *colors;
1575: PetscBool dense = PETSC_FALSE;
1576: ISColoringValue maxcol = 0;
1577: DM_Composite *com = (DM_Composite*)dm->data;
1581: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1582: else if (ctype == IS_COLORING_GLOBAL) {
1583: n = com->n;
1584: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1585: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1587: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1588: if (dense) {
1589: for (i=0; i<n; i++) {
1590: colors[i] = (ISColoringValue)(com->rstart + i);
1591: }
1592: maxcol = com->N;
1593: } else {
1594: struct DMCompositeLink *next = com->next;
1595: PetscMPIInt rank;
1597: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1598: cnt = 0;
1599: while (next) {
1600: ISColoring lcoloring;
1602: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1603: for (i=0; i<lcoloring->N; i++) {
1604: colors[cnt++] = maxcol + lcoloring->colors[i];
1605: }
1606: maxcol += lcoloring->n;
1607: ISColoringDestroy(&lcoloring);
1608: next = next->next;
1609: }
1610: }
1611: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1612: return(0);
1613: }
1615: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1616: {
1617: PetscErrorCode ierr;
1618: struct DMCompositeLink *next;
1619: PetscScalar *garray,*larray;
1620: DM_Composite *com = (DM_Composite*)dm->data;
1626: if (!com->setup) {
1627: DMSetUp(dm);
1628: }
1630: VecGetArray(gvec,&garray);
1631: VecGetArray(lvec,&larray);
1633: /* loop over packed objects, handling one at at time */
1634: next = com->next;
1635: while (next) {
1636: Vec local,global;
1637: PetscInt N;
1639: DMGetGlobalVector(next->dm,&global);
1640: VecGetLocalSize(global,&N);
1641: VecPlaceArray(global,garray);
1642: DMGetLocalVector(next->dm,&local);
1643: VecPlaceArray(local,larray);
1644: DMGlobalToLocalBegin(next->dm,global,mode,local);
1645: DMGlobalToLocalEnd(next->dm,global,mode,local);
1646: VecResetArray(global);
1647: VecResetArray(local);
1648: DMRestoreGlobalVector(next->dm,&global);
1649: DMRestoreLocalVector(next->dm,&local);
1651: larray += next->nlocal;
1652: garray += next->n;
1653: next = next->next;
1654: }
1656: VecRestoreArray(gvec,NULL);
1657: VecRestoreArray(lvec,NULL);
1658: return(0);
1659: }
1661: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1662: {
1667: return(0);
1668: }
1670: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1671: {
1672: PetscErrorCode ierr;
1673: struct DMCompositeLink *next;
1674: PetscScalar *larray,*garray;
1675: DM_Composite *com = (DM_Composite*)dm->data;
1682: if (!com->setup) {
1683: DMSetUp(dm);
1684: }
1686: VecGetArray(lvec,&larray);
1687: VecGetArray(gvec,&garray);
1689: /* loop over packed objects, handling one at at time */
1690: next = com->next;
1691: while (next) {
1692: Vec global,local;
1694: DMGetLocalVector(next->dm,&local);
1695: VecPlaceArray(local,larray);
1696: DMGetGlobalVector(next->dm,&global);
1697: VecPlaceArray(global,garray);
1698: DMLocalToGlobalBegin(next->dm,local,mode,global);
1699: DMLocalToGlobalEnd(next->dm,local,mode,global);
1700: VecResetArray(local);
1701: VecResetArray(global);
1702: DMRestoreGlobalVector(next->dm,&global);
1703: DMRestoreLocalVector(next->dm,&local);
1705: garray += next->n;
1706: larray += next->nlocal;
1707: next = next->next;
1708: }
1710: VecRestoreArray(gvec,NULL);
1711: VecRestoreArray(lvec,NULL);
1713: return(0);
1714: }
1716: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1717: {
1722: return(0);
1723: }
1725: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1726: {
1727: PetscErrorCode ierr;
1728: struct DMCompositeLink *next;
1729: PetscScalar *array1,*array2;
1730: DM_Composite *com = (DM_Composite*)dm->data;
1737: if (!com->setup) {
1738: DMSetUp(dm);
1739: }
1741: VecGetArray(vec1,&array1);
1742: VecGetArray(vec2,&array2);
1744: /* loop over packed objects, handling one at at time */
1745: next = com->next;
1746: while (next) {
1747: Vec local1,local2;
1749: DMGetLocalVector(next->dm,&local1);
1750: VecPlaceArray(local1,array1);
1751: DMGetLocalVector(next->dm,&local2);
1752: VecPlaceArray(local2,array2);
1753: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1754: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1755: VecResetArray(local2);
1756: DMRestoreLocalVector(next->dm,&local2);
1757: VecResetArray(local1);
1758: DMRestoreLocalVector(next->dm,&local1);
1760: array1 += next->nlocal;
1761: array2 += next->nlocal;
1762: next = next->next;
1763: }
1765: VecRestoreArray(vec1,NULL);
1766: VecRestoreArray(vec2,NULL);
1768: return(0);
1769: }
1771: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1772: {
1777: return(0);
1778: }
1780: /*MC
1781: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1783: Level: intermediate
1785: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1786: M*/
1789: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1790: {
1792: DM_Composite *com;
1795: PetscNewLog(p,&com);
1796: p->data = com;
1797: com->n = 0;
1798: com->nghost = 0;
1799: com->next = NULL;
1800: com->nDM = 0;
1802: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1803: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1804: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1805: p->ops->createfieldis = DMCreateFieldIS_Composite;
1806: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1807: p->ops->refine = DMRefine_Composite;
1808: p->ops->coarsen = DMCoarsen_Composite;
1809: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1810: p->ops->creatematrix = DMCreateMatrix_Composite;
1811: p->ops->getcoloring = DMCreateColoring_Composite;
1812: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1813: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1814: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1815: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1816: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1817: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1818: p->ops->destroy = DMDestroy_Composite;
1819: p->ops->view = DMView_Composite;
1820: p->ops->setup = DMSetUp_Composite;
1822: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1823: return(0);
1824: }
1826: /*@
1827: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1828: vectors made up of several subvectors.
1830: Collective
1832: Input Parameter:
1833: . comm - the processors that will share the global vector
1835: Output Parameters:
1836: . packer - the packer object
1838: Level: advanced
1840: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1841: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1842: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1844: @*/
1845: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1846: {
1851: DMCreate(comm,packer);
1852: DMSetType(*packer,DMCOMPOSITE);
1853: return(0);
1854: }