Actual source code: pack.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  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,&ltog);
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,&ltogs);
1371:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1372:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[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: }