Actual source code: isltog.c
1: #define PETSCVEC_DLL
3: #include petscvec.h
4: #include src/vec/is/isimpl.h
6: PetscCookie IS_LTOGM_COOKIE = -1;
10: /*@C
11: ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping.
13: Not Collective
15: Input Parameter:
16: . ltog - local to global mapping
18: Output Parameter:
19: . n - the number of entries in the local mapping
21: Level: advanced
23: Concepts: mapping^local to global
25: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
26: @*/
27: PetscErrorCode ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping,PetscInt *n)
28: {
32: *n = mapping->n;
33: return(0);
34: }
38: /*@C
39: ISLocalToGlobalMappingView - View a local to global mapping
41: Not Collective
43: Input Parameters:
44: + ltog - local to global mapping
45: - viewer - viewer
47: Level: advanced
49: Concepts: mapping^local to global
51: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
52: @*/
53: PetscErrorCode ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping,PetscViewer viewer)
54: {
55: PetscInt i;
56: PetscMPIInt rank;
57: PetscTruth iascii;
58: PetscErrorCode ierr;
62: if (!viewer) {
63: PetscViewerASCIIGetStdout(mapping->comm,&viewer);
64: }
67: MPI_Comm_rank(mapping->comm,&rank);
68: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
69: if (iascii) {
70: for (i=0; i<mapping->n; i++) {
71: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d %d\n",rank,i,mapping->indices[i]);
72: }
73: PetscViewerFlush(viewer);
74: } else {
75: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for ISLocalToGlobalMapping",((PetscObject)viewer)->type_name);
76: }
78: return(0);
79: }
83: /*@
84: ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
85: ordering and a global parallel ordering.
87: Not collective
89: Input Parameter:
90: . is - index set containing the global numbers for each local
92: Output Parameter:
93: . mapping - new mapping data structure
95: Level: advanced
97: Concepts: mapping^local to global
99: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
100: @*/
101: PetscErrorCode ISLocalToGlobalMappingCreateIS(IS is,ISLocalToGlobalMapping *mapping)
102: {
104: PetscInt n,*indices;
105: MPI_Comm comm;
111: PetscObjectGetComm((PetscObject)is,&comm);
112: ISGetLocalSize(is,&n);
113: ISGetIndices(is,&indices);
114: ISLocalToGlobalMappingCreate(comm,n,indices,mapping);
115: ISRestoreIndices(is,&indices);
117: return(0);
118: }
123: /*@
124: ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
125: ordering and a global parallel ordering.
127: Not Collective, but communicator may have more than one process
129: Input Parameters:
130: + comm - MPI communicator
131: . n - the number of local elements
132: - indices - the global index for each local element
134: Output Parameter:
135: . mapping - new mapping data structure
137: Level: advanced
139: Concepts: mapping^local to global
141: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreateNC()
142: @*/
143: PetscErrorCode ISLocalToGlobalMappingCreate(MPI_Comm cm,PetscInt n,const PetscInt indices[],ISLocalToGlobalMapping *mapping)
144: {
146: PetscInt *in;
151: PetscMalloc(n*sizeof(PetscInt),&in);
152: PetscMemcpy(in,indices,n*sizeof(PetscInt));
153: ISLocalToGlobalMappingCreateNC(cm,n,in,mapping);
154: return(0);
155: }
159: /*@C
160: ISLocalToGlobalMappingCreateNC - Creates a mapping between a local (0 to n)
161: ordering and a global parallel ordering.
163: Not Collective, but communicator may have more than one process
165: Input Parameters:
166: + comm - MPI communicator
167: . n - the number of local elements
168: - indices - the global index for each local element
170: Output Parameter:
171: . mapping - new mapping data structure
173: Level: developer
175: Notes: Does not copy the indices, just keeps the pointer to the indices. The ISLocalToGlobalMappingDestroy()
176: will free the space so it must be obtained with PetscMalloc() and it must not be freed elsewhere.
178: Concepts: mapping^local to global
180: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate()
181: @*/
182: PetscErrorCode ISLocalToGlobalMappingCreateNC(MPI_Comm cm,PetscInt n,const PetscInt indices[],ISLocalToGlobalMapping *mapping)
183: {
187: if (n) {
189: }
191: *mapping = PETSC_NULL;
192: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
193: VecInitializePackage(PETSC_NULL);
194: #endif
195: if (IS_LTOGM_COOKIE == -1) {
196: PetscLogClassRegister(&IS_LTOGM_COOKIE,"IS Local to global mapping");
197: }
199: PetscHeaderCreate(*mapping,_p_ISLocalToGlobalMapping,int,IS_LTOGM_COOKIE,0,"ISLocalToGlobalMapping",
200: cm,ISLocalToGlobalMappingDestroy,ISLocalToGlobalMappingView);
201: PetscLogObjectMemory(*mapping,sizeof(struct _p_ISLocalToGlobalMapping)+n*sizeof(PetscInt));
203: (*mapping)->n = n;
204: (*mapping)->indices = (PetscInt*)indices;
206: /*
207: Do not create the global to local mapping. This is only created if
208: ISGlobalToLocalMapping() is called
209: */
210: (*mapping)->globals = 0;
211: return(0);
212: }
216: /*@
217: ISLocalToGlobalMappingBlock - Creates a blocked index version of an
218: ISLocalToGlobalMapping that is appropriate for MatSetLocalToGlobalMappingBlock()
219: and VecSetLocalToGlobalMappingBlock().
221: Not Collective, but communicator may have more than one process
223: Input Parameters:
224: + inmap - original point-wise mapping
225: - bs - block size
227: Output Parameter:
228: . outmap - block based mapping; the indices are relative to BLOCKS, not individual vector or matrix entries.
230: Level: advanced
232: Concepts: mapping^local to global
234: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingCreateIS()
235: @*/
236: PetscErrorCode ISLocalToGlobalMappingBlock(ISLocalToGlobalMapping inmap,PetscInt bs,ISLocalToGlobalMapping *outmap)
237: {
239: PetscInt *ii,i,n;
242: if (bs > 1) {
243: n = inmap->n/bs;
244: if (n*bs != inmap->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Pointwise mapping length is not divisible by block size");
245: PetscMalloc(n*sizeof(PetscInt),&ii);
246: for (i=0; i<n; i++) {
247: ii[i] = inmap->indices[bs*i]/bs;
248: }
249: ISLocalToGlobalMappingCreate(inmap->comm,n,ii,outmap);
250: PetscFree(ii);
251: } else {
252: *outmap = inmap;
253: PetscObjectReference((PetscObject)inmap);
254: }
255: return(0);
256: }
257:
260: /*@
261: ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
262: ordering and a global parallel ordering.
264: Note Collective
266: Input Parameters:
267: . mapping - mapping data structure
269: Level: advanced
271: .seealso: ISLocalToGlobalMappingCreate()
272: @*/
273: PetscErrorCode ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping mapping)
274: {
278: if (--mapping->refct > 0) return(0);
279: if (mapping->refct < 0) {
280: SETERRQ(PETSC_ERR_PLIB,"Mapping already destroyed");
281: }
283: PetscFree(mapping->indices);
284: PetscFree(mapping->globals);
285: PetscHeaderDestroy(mapping);
286: return(0);
287: }
288:
291: /*@
292: ISLocalToGlobalMappingApplyIS - Creates from an IS in the local numbering
293: a new index set using the global numbering defined in an ISLocalToGlobalMapping
294: context.
296: Not collective
298: Input Parameters:
299: + mapping - mapping between local and global numbering
300: - is - index set in local numbering
302: Output Parameters:
303: . newis - index set in global numbering
305: Level: advanced
307: Concepts: mapping^local to global
309: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
310: ISLocalToGlobalMappingDestroy(), ISGlobalToLocalMappingApply()
311: @*/
312: PetscErrorCode ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping,IS is,IS *newis)
313: {
315: PetscInt n,i,*idxin,*idxmap,*idxout,Nmax = mapping->n;
322: ISGetLocalSize(is,&n);
323: ISGetIndices(is,&idxin);
324: idxmap = mapping->indices;
325:
326: PetscMalloc(n*sizeof(PetscInt),&idxout);
327: for (i=0; i<n; i++) {
328: if (idxin[i] >= Nmax) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Local index %d too large %d (max) at %d",idxin[i],Nmax-1,i);
329: idxout[i] = idxmap[idxin[i]];
330: }
331: ISRestoreIndices(is,&idxin);
332: ISCreateGeneral(PETSC_COMM_SELF,n,idxout,newis);
333: PetscFree(idxout);
334: return(0);
335: }
337: /*MC
338: ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
339: and converts them to the global numbering.
341: Not collective
343: Input Parameters:
344: + mapping - the local to global mapping context
345: . N - number of integers
346: - in - input indices in local numbering
348: Output Parameter:
349: . out - indices in global numbering
351: Synopsis:
352: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping,int N,int in[],int out[])
354: Notes:
355: The in and out array parameters may be identical.
357: Level: advanced
359: .seealso: ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
360: ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
361: AOPetscToApplication(), ISGlobalToLocalMappingApply()
363: Concepts: mapping^local to global
365: M*/
367: /* -----------------------------------------------------------------------------------------*/
371: /*
372: Creates the global fields in the ISLocalToGlobalMapping structure
373: */
374: static PetscErrorCode ISGlobalToLocalMappingSetUp_Private(ISLocalToGlobalMapping mapping)
375: {
377: PetscInt i,*idx = mapping->indices,n = mapping->n,end,start,*globals;
380: end = 0;
381: start = 100000000;
383: for (i=0; i<n; i++) {
384: if (idx[i] < 0) continue;
385: if (idx[i] < start) start = idx[i];
386: if (idx[i] > end) end = idx[i];
387: }
388: if (start > end) {start = 0; end = -1;}
389: mapping->globalstart = start;
390: mapping->globalend = end;
392: PetscMalloc((end-start+2)*sizeof(PetscInt),&globals);
393: mapping->globals = globals;
394: for (i=0; i<end-start+1; i++) {
395: globals[i] = -1;
396: }
397: for (i=0; i<n; i++) {
398: if (idx[i] < 0) continue;
399: globals[idx[i] - start] = i;
400: }
402: PetscLogObjectMemory(mapping,(end-start+1)*sizeof(PetscInt));
403: return(0);
404: }
408: /*@
409: ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
410: specified with a global numbering.
412: Not collective
414: Input Parameters:
415: + mapping - mapping between local and global numbering
416: . type - IS_GTOLM_MASK - replaces global indices with no local value with -1
417: IS_GTOLM_DROP - drops the indices with no local value from the output list
418: . n - number of global indices to map
419: - idx - global indices to map
421: Output Parameters:
422: + nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
423: - idxout - local index of each global index, one must pass in an array long enough
424: to hold all the indices. You can call ISGlobalToLocalMappingApply() with
425: idxout == PETSC_NULL to determine the required length (returned in nout)
426: and then allocate the required space and call ISGlobalToLocalMappingApply()
427: a second time to set the values.
429: Notes:
430: Either nout or idxout may be PETSC_NULL. idx and idxout may be identical.
432: This is not scalable in memory usage. Each processor requires O(Nglobal) size
433: array to compute these.
435: Level: advanced
437: Concepts: mapping^global to local
439: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
440: ISLocalToGlobalMappingDestroy()
441: @*/
442: PetscErrorCode ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingType type,
443: PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
444: {
445: PetscInt i,*globals,nf = 0,tmp,start,end;
449: if (!mapping->globals) {
450: ISGlobalToLocalMappingSetUp_Private(mapping);
451: }
452: globals = mapping->globals;
453: start = mapping->globalstart;
454: end = mapping->globalend;
456: if (type == IS_GTOLM_MASK) {
457: if (idxout) {
458: for (i=0; i<n; i++) {
459: if (idx[i] < 0) idxout[i] = idx[i];
460: else if (idx[i] < start) idxout[i] = -1;
461: else if (idx[i] > end) idxout[i] = -1;
462: else idxout[i] = globals[idx[i] - start];
463: }
464: }
465: if (nout) *nout = n;
466: } else {
467: if (idxout) {
468: for (i=0; i<n; i++) {
469: if (idx[i] < 0) continue;
470: if (idx[i] < start) continue;
471: if (idx[i] > end) continue;
472: tmp = globals[idx[i] - start];
473: if (tmp < 0) continue;
474: idxout[nf++] = tmp;
475: }
476: } else {
477: for (i=0; i<n; i++) {
478: if (idx[i] < 0) continue;
479: if (idx[i] < start) continue;
480: if (idx[i] > end) continue;
481: tmp = globals[idx[i] - start];
482: if (tmp < 0) continue;
483: nf++;
484: }
485: }
486: if (nout) *nout = nf;
487: }
489: return(0);
490: }
494: /*@C
495: ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and
496: each index shared by more than one processor
498: Collective on ISLocalToGlobalMapping
500: Input Parameters:
501: . mapping - the mapping from local to global indexing
503: Output Parameter:
504: + nproc - number of processors that are connected to this one
505: . proc - neighboring processors
506: . numproc - number of indices for each subdomain (processor)
507: - indices - indices of local nodes shared with neighbor (sorted by global numbering)
509: Level: advanced
511: Concepts: mapping^local to global
513: Fortran Usage:
514: $ ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
515: $ ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
516: PetscInt indices[nproc][numprocmax],ierr)
517: There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
518: indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.
521: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
522: ISLocalToGlobalMappingRestoreInfo()
523: @*/
524: PetscErrorCode ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
525: {
527: PetscMPIInt size,rank,tag1,tag2,tag3,*len,*source,imdex;
528: PetscInt i,n = mapping->n,Ng,ng,max = 0,*lindices = mapping->indices;
529: PetscInt *nprocs,*owner,nsends,*sends,j,*starts,nmax,nrecvs,*recvs,proc;
530: PetscInt cnt,scale,*ownedsenders,*nownedsenders,rstart,nowned;
531: PetscInt node,nownedm,nt,*sends2,nsends2,*starts2,*lens2,*dest,nrecvs2,*starts3,*recvs2,k,*bprocs,*tmp;
532: PetscInt first_procs,first_numprocs,*first_indices;
533: MPI_Request *recv_waits,*send_waits;
534: MPI_Status recv_status,*send_status,*recv_statuses;
535: MPI_Comm comm = mapping->comm;
536: PetscTruth debug = PETSC_FALSE;
539: MPI_Comm_size(comm,&size);
540: MPI_Comm_rank(comm,&rank);
541: if (size == 1) {
542: *nproc = 0;
543: *procs = PETSC_NULL;
544: PetscMalloc(sizeof(PetscInt),numprocs);
545: (*numprocs)[0] = 0;
546: PetscMalloc(sizeof(PetscInt*),indices);
547: (*indices)[0] = PETSC_NULL;
548: return(0);
549: }
551: PetscOptionsHasName(PETSC_NULL,"-islocaltoglobalmappinggetinfo_debug",&debug);
553: /*
554: Notes on ISLocalToGlobalMappingGetInfo
556: globally owned node - the nodes that have been assigned to this processor in global
557: numbering, just for this routine.
559: nontrivial globally owned node - node assigned to this processor that is on a subdomain
560: boundary (i.e. is has more than one local owner)
562: locally owned node - node that exists on this processors subdomain
564: nontrivial locally owned node - node that is not in the interior (i.e. has more than one
565: local subdomain
566: */
567: PetscObjectGetNewTag((PetscObject)mapping,&tag1);
568: PetscObjectGetNewTag((PetscObject)mapping,&tag2);
569: PetscObjectGetNewTag((PetscObject)mapping,&tag3);
571: for (i=0; i<n; i++) {
572: if (lindices[i] > max) max = lindices[i];
573: }
574: MPI_Allreduce(&max,&Ng,1,MPIU_INT,MPI_MAX,comm);
575: Ng++;
576: MPI_Comm_size(comm,&size);
577: MPI_Comm_rank(comm,&rank);
578: scale = Ng/size + 1;
579: ng = scale; if (rank == size-1) ng = Ng - scale*(size-1); ng = PetscMax(1,ng);
580: rstart = scale*rank;
582: /* determine ownership ranges of global indices */
583: PetscMalloc(2*size*sizeof(PetscInt),&nprocs);
584: PetscMemzero(nprocs,2*size*sizeof(PetscInt));
586: /* determine owners of each local node */
587: PetscMalloc(n*sizeof(PetscInt),&owner);
588: for (i=0; i<n; i++) {
589: proc = lindices[i]/scale; /* processor that globally owns this index */
590: nprocs[2*proc+1] = 1; /* processor globally owns at least one of ours */
591: owner[i] = proc;
592: nprocs[2*proc]++; /* count of how many that processor globally owns of ours */
593: }
594: nsends = 0; for (i=0; i<size; i++) nsends += nprocs[2*i+1];
595: PetscInfo1(0,"Number of global owners for my local data %d\n",nsends);
597: /* inform other processors of number of messages and max length*/
598: PetscMaxSum(comm,nprocs,&nmax,&nrecvs);
599: PetscInfo1(0,"Number of local owners for my global data %d\n",nrecvs);
601: /* post receives for owned rows */
602: PetscMalloc((2*nrecvs+1)*(nmax+1)*sizeof(PetscInt),&recvs);
603: PetscMalloc((nrecvs+1)*sizeof(MPI_Request),&recv_waits);
604: for (i=0; i<nrecvs; i++) {
605: MPI_Irecv(recvs+2*nmax*i,2*nmax,MPIU_INT,MPI_ANY_SOURCE,tag1,comm,recv_waits+i);
606: }
608: /* pack messages containing lists of local nodes to owners */
609: PetscMalloc((2*n+1)*sizeof(PetscInt),&sends);
610: PetscMalloc((size+1)*sizeof(PetscInt),&starts);
611: starts[0] = 0;
612: for (i=1; i<size; i++) { starts[i] = starts[i-1] + 2*nprocs[2*i-2];}
613: for (i=0; i<n; i++) {
614: sends[starts[owner[i]]++] = lindices[i];
615: sends[starts[owner[i]]++] = i;
616: }
617: PetscFree(owner);
618: starts[0] = 0;
619: for (i=1; i<size; i++) { starts[i] = starts[i-1] + 2*nprocs[2*i-2];}
621: /* send the messages */
622: PetscMalloc((nsends+1)*sizeof(MPI_Request),&send_waits);
623: PetscMalloc((nsends+1)*sizeof(PetscInt),&dest);
624: cnt = 0;
625: for (i=0; i<size; i++) {
626: if (nprocs[2*i]) {
627: MPI_Isend(sends+starts[i],2*nprocs[2*i],MPIU_INT,i,tag1,comm,send_waits+cnt);
628: dest[cnt] = i;
629: cnt++;
630: }
631: }
632: PetscFree(starts);
634: /* wait on receives */
635: PetscMalloc((nrecvs+1)*sizeof(PetscMPIInt),&source);
636: PetscMalloc((nrecvs+1)*sizeof(PetscMPIInt),&len);
637: cnt = nrecvs;
638: PetscMalloc((ng+1)*sizeof(PetscInt),&nownedsenders);
639: PetscMemzero(nownedsenders,ng*sizeof(PetscInt));
640: while (cnt) {
641: MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
642: /* unpack receives into our local space */
643: MPI_Get_count(&recv_status,MPIU_INT,&len[imdex]);
644: source[imdex] = recv_status.MPI_SOURCE;
645: len[imdex] = len[imdex]/2;
646: /* count how many local owners for each of my global owned indices */
647: for (i=0; i<len[imdex]; i++) nownedsenders[recvs[2*imdex*nmax+2*i]-rstart]++;
648: cnt--;
649: }
650: PetscFree(recv_waits);
652: /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
653: nowned = 0;
654: nownedm = 0;
655: for (i=0; i<ng; i++) {
656: if (nownedsenders[i] > 1) {nownedm += nownedsenders[i]; nowned++;}
657: }
659: /* create single array to contain rank of all local owners of each globally owned index */
660: PetscMalloc((nownedm+1)*sizeof(PetscInt),&ownedsenders);
661: PetscMalloc((ng+1)*sizeof(PetscInt),&starts);
662: starts[0] = 0;
663: for (i=1; i<ng; i++) {
664: if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
665: else starts[i] = starts[i-1];
666: }
668: /* for each nontrival globally owned node list all arriving processors */
669: for (i=0; i<nrecvs; i++) {
670: for (j=0; j<len[i]; j++) {
671: node = recvs[2*i*nmax+2*j]-rstart;
672: if (nownedsenders[node] > 1) {
673: ownedsenders[starts[node]++] = source[i];
674: }
675: }
676: }
678: if (debug) { /* ----------------------------------- */
679: starts[0] = 0;
680: for (i=1; i<ng; i++) {
681: if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
682: else starts[i] = starts[i-1];
683: }
684: for (i=0; i<ng; i++) {
685: if (nownedsenders[i] > 1) {
686: PetscSynchronizedPrintf(comm,"[%d] global node %d local owner processors: ",rank,i+rstart);
687: for (j=0; j<nownedsenders[i]; j++) {
688: PetscSynchronizedPrintf(comm,"%d ",ownedsenders[starts[i]+j]);
689: }
690: PetscSynchronizedPrintf(comm,"\n");
691: }
692: }
693: PetscSynchronizedFlush(comm);
694: }/* ----------------------------------- */
696: /* wait on original sends */
697: if (nsends) {
698: PetscMalloc(nsends*sizeof(MPI_Status),&send_status);
699: MPI_Waitall(nsends,send_waits,send_status);
700: PetscFree(send_status);
701: }
702: PetscFree(send_waits);
703: PetscFree(sends);
704: PetscFree(nprocs);
706: /* pack messages to send back to local owners */
707: starts[0] = 0;
708: for (i=1; i<ng; i++) {
709: if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
710: else starts[i] = starts[i-1];
711: }
712: nsends2 = nrecvs;
713: PetscMalloc((nsends2+1)*sizeof(PetscInt),&nprocs); /* length of each message */
714: for (i=0; i<nrecvs; i++) {
715: nprocs[i] = 1;
716: for (j=0; j<len[i]; j++) {
717: node = recvs[2*i*nmax+2*j]-rstart;
718: if (nownedsenders[node] > 1) {
719: nprocs[i] += 2 + nownedsenders[node];
720: }
721: }
722: }
723: nt = 0; for (i=0; i<nsends2; i++) nt += nprocs[i];
724: PetscMalloc((nt+1)*sizeof(PetscInt),&sends2);
725: PetscMalloc((nsends2+1)*sizeof(PetscInt),&starts2);
726: starts2[0] = 0; for (i=1; i<nsends2; i++) starts2[i] = starts2[i-1] + nprocs[i-1];
727: /*
728: Each message is 1 + nprocs[i] long, and consists of
729: (0) the number of nodes being sent back
730: (1) the local node number,
731: (2) the number of processors sharing it,
732: (3) the processors sharing it
733: */
734: for (i=0; i<nsends2; i++) {
735: cnt = 1;
736: sends2[starts2[i]] = 0;
737: for (j=0; j<len[i]; j++) {
738: node = recvs[2*i*nmax+2*j]-rstart;
739: if (nownedsenders[node] > 1) {
740: sends2[starts2[i]]++;
741: sends2[starts2[i]+cnt++] = recvs[2*i*nmax+2*j+1];
742: sends2[starts2[i]+cnt++] = nownedsenders[node];
743: PetscMemcpy(&sends2[starts2[i]+cnt],&ownedsenders[starts[node]],nownedsenders[node]*sizeof(PetscInt));
744: cnt += nownedsenders[node];
745: }
746: }
747: }
749: /* receive the message lengths */
750: nrecvs2 = nsends;
751: PetscMalloc((nrecvs2+1)*sizeof(PetscInt),&lens2);
752: PetscMalloc((nrecvs2+1)*sizeof(PetscInt),&starts3);
753: PetscMalloc((nrecvs2+1)*sizeof(MPI_Request),&recv_waits);
754: for (i=0; i<nrecvs2; i++) {
755: MPI_Irecv(&lens2[i],1,MPIU_INT,dest[i],tag2,comm,recv_waits+i);
756: }
758: /* send the message lengths */
759: for (i=0; i<nsends2; i++) {
760: MPI_Send(&nprocs[i],1,MPIU_INT,source[i],tag2,comm);
761: }
763: /* wait on receives of lens */
764: if (nrecvs2) {
765: PetscMalloc(nrecvs2*sizeof(MPI_Status),&recv_statuses);
766: MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
767: PetscFree(recv_statuses);
768: }
769: PetscFree(recv_waits);
771: starts3[0] = 0;
772: nt = 0;
773: for (i=0; i<nrecvs2-1; i++) {
774: starts3[i+1] = starts3[i] + lens2[i];
775: nt += lens2[i];
776: }
777: nt += lens2[nrecvs2-1];
779: PetscMalloc((nt+1)*sizeof(PetscInt),&recvs2);
780: PetscMalloc((nrecvs2+1)*sizeof(MPI_Request),&recv_waits);
781: for (i=0; i<nrecvs2; i++) {
782: MPI_Irecv(recvs2+starts3[i],lens2[i],MPIU_INT,dest[i],tag3,comm,recv_waits+i);
783: }
784:
785: /* send the messages */
786: PetscMalloc((nsends2+1)*sizeof(MPI_Request),&send_waits);
787: for (i=0; i<nsends2; i++) {
788: MPI_Isend(sends2+starts2[i],nprocs[i],MPIU_INT,source[i],tag3,comm,send_waits+i);
789: }
791: /* wait on receives */
792: if (nrecvs2) {
793: PetscMalloc(nrecvs2*sizeof(MPI_Status),&recv_statuses);
794: MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
795: PetscFree(recv_statuses);
796: }
797: PetscFree(recv_waits);
798: PetscFree(nprocs);
800: if (debug) { /* ----------------------------------- */
801: cnt = 0;
802: for (i=0; i<nrecvs2; i++) {
803: nt = recvs2[cnt++];
804: for (j=0; j<nt; j++) {
805: PetscSynchronizedPrintf(comm,"[%d] local node %d number of subdomains %d: ",rank,recvs2[cnt],recvs2[cnt+1]);
806: for (k=0; k<recvs2[cnt+1]; k++) {
807: PetscSynchronizedPrintf(comm,"%d ",recvs2[cnt+2+k]);
808: }
809: cnt += 2 + recvs2[cnt+1];
810: PetscSynchronizedPrintf(comm,"\n");
811: }
812: }
813: PetscSynchronizedFlush(comm);
814: } /* ----------------------------------- */
816: /* count number subdomains for each local node */
817: PetscMalloc(size*sizeof(PetscInt),&nprocs);
818: PetscMemzero(nprocs,size*sizeof(PetscInt));
819: cnt = 0;
820: for (i=0; i<nrecvs2; i++) {
821: nt = recvs2[cnt++];
822: for (j=0; j<nt; j++) {
823: for (k=0; k<recvs2[cnt+1]; k++) {
824: nprocs[recvs2[cnt+2+k]]++;
825: }
826: cnt += 2 + recvs2[cnt+1];
827: }
828: }
829: nt = 0; for (i=0; i<size; i++) nt += (nprocs[i] > 0);
830: *nproc = nt;
831: PetscMalloc((nt+1)*sizeof(PetscInt),procs);
832: PetscMalloc((nt+1)*sizeof(PetscInt),numprocs);
833: PetscMalloc((nt+1)*sizeof(PetscInt*),indices);
834: PetscMalloc(size*sizeof(PetscInt),&bprocs);
835: cnt = 0;
836: for (i=0; i<size; i++) {
837: if (nprocs[i] > 0) {
838: bprocs[i] = cnt;
839: (*procs)[cnt] = i;
840: (*numprocs)[cnt] = nprocs[i];
841: PetscMalloc(nprocs[i]*sizeof(PetscInt),&(*indices)[cnt]);
842: cnt++;
843: }
844: }
846: /* make the list of subdomains for each nontrivial local node */
847: PetscMemzero(*numprocs,nt*sizeof(PetscInt));
848: cnt = 0;
849: for (i=0; i<nrecvs2; i++) {
850: nt = recvs2[cnt++];
851: for (j=0; j<nt; j++) {
852: for (k=0; k<recvs2[cnt+1]; k++) {
853: (*indices)[bprocs[recvs2[cnt+2+k]]][(*numprocs)[bprocs[recvs2[cnt+2+k]]]++] = recvs2[cnt];
854: }
855: cnt += 2 + recvs2[cnt+1];
856: }
857: }
858: PetscFree(bprocs);
859: PetscFree(recvs2);
861: /* sort the node indexing by their global numbers */
862: nt = *nproc;
863: for (i=0; i<nt; i++) {
864: PetscMalloc(((*numprocs)[i])*sizeof(PetscInt),&tmp);
865: for (j=0; j<(*numprocs)[i]; j++) {
866: tmp[j] = lindices[(*indices)[i][j]];
867: }
868: PetscSortIntWithArray((*numprocs)[i],tmp,(*indices)[i]);
869: PetscFree(tmp);
870: }
872: if (debug) { /* ----------------------------------- */
873: nt = *nproc;
874: for (i=0; i<nt; i++) {
875: PetscSynchronizedPrintf(comm,"[%d] subdomain %d number of indices %d: ",rank,(*procs)[i],(*numprocs)[i]);
876: for (j=0; j<(*numprocs)[i]; j++) {
877: PetscSynchronizedPrintf(comm,"%d ",(*indices)[i][j]);
878: }
879: PetscSynchronizedPrintf(comm,"\n");
880: }
881: PetscSynchronizedFlush(comm);
882: } /* ----------------------------------- */
884: /* wait on sends */
885: if (nsends2) {
886: PetscMalloc(nsends2*sizeof(MPI_Status),&send_status);
887: MPI_Waitall(nsends2,send_waits,send_status);
888: PetscFree(send_status);
889: }
891: PetscFree(starts3);
892: PetscFree(dest);
893: PetscFree(send_waits);
895: PetscFree(nownedsenders);
896: PetscFree(ownedsenders);
897: PetscFree(starts);
898: PetscFree(starts2);
899: PetscFree(lens2);
901: PetscFree(source);
902: PetscFree(len);
903: PetscFree(recvs);
904: PetscFree(nprocs);
905: PetscFree(sends2);
907: /* put the information about myself as the first entry in the list */
908: first_procs = (*procs)[0];
909: first_numprocs = (*numprocs)[0];
910: first_indices = (*indices)[0];
911: for (i=0; i<*nproc; i++) {
912: if ((*procs)[i] == rank) {
913: (*procs)[0] = (*procs)[i];
914: (*numprocs)[0] = (*numprocs)[i];
915: (*indices)[0] = (*indices)[i];
916: (*procs)[i] = first_procs;
917: (*numprocs)[i] = first_numprocs;
918: (*indices)[i] = first_indices;
919: break;
920: }
921: }
922: return(0);
923: }
927: /*@C
928: ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by ISLocalToGlobalMappingGetInfo()
930: Collective on ISLocalToGlobalMapping
932: Input Parameters:
933: . mapping - the mapping from local to global indexing
935: Output Parameter:
936: + nproc - number of processors that are connected to this one
937: . proc - neighboring processors
938: . numproc - number of indices for each processor
939: - indices - indices of local nodes shared with neighbor (sorted by global numbering)
941: Level: advanced
943: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
944: ISLocalToGlobalMappingGetInfo()
945: @*/
946: PetscErrorCode ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
947: {
949: PetscInt i;
952: PetscFree(*procs);
953: PetscFree(*numprocs);
954: if (*indices) {
955: PetscFree((*indices)[0]);
956: for (i=1; i<*nproc; i++) {
957: PetscFree((*indices)[i]);
958: }
959: PetscFree(*indices);
960: }
961: return(0);
962: }