Actual source code: isdiff.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/isimpl.h>                    /*I "petscis.h"  I*/
  3: #include <petscbt.h>

  7: /*@
  8:    ISDifference - Computes the difference between two index sets.

 10:    Collective on IS

 12:    Input Parameter:
 13: +  is1 - first index, to have items removed from it
 14: -  is2 - index values to be removed

 16:    Output Parameters:
 17: .  isout - is1 - is2

 19:    Notes:
 20:    Negative values are removed from the lists. is2 may have values
 21:    that are not in is1. This requires O(imax-imin) memory and O(imax-imin)
 22:    work, where imin and imax are the bounds on the indices in is1.

 24:    Level: intermediate

 26:    Concepts: index sets^difference
 27:    Concepts: IS^difference

 29: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()

 31: @*/
 32: PetscErrorCode  ISDifference(IS is1,IS is2,IS *isout)
 33: {
 35:   PetscInt       i,n1,n2,imin,imax,nout,*iout;
 36:   const PetscInt *i1,*i2;
 37:   PetscBT        mask;
 38:   MPI_Comm       comm;


 45:   ISGetIndices(is1,&i1);
 46:   ISGetLocalSize(is1,&n1);

 48:   /* Create a bit mask array to contain required values */
 49:   if (n1) {
 50:     imin = PETSC_MAX_INT;
 51:     imax = 0;
 52:     for (i=0; i<n1; i++) {
 53:       if (i1[i] < 0) continue;
 54:       imin = PetscMin(imin,i1[i]);
 55:       imax = PetscMax(imax,i1[i]);
 56:     }
 57:   } else imin = imax = 0;

 59:   PetscBTCreate(imax-imin,&mask);
 60:   /* Put the values from is1 */
 61:   for (i=0; i<n1; i++) {
 62:     if (i1[i] < 0) continue;
 63:     PetscBTSet(mask,i1[i] - imin);
 64:   }
 65:   ISRestoreIndices(is1,&i1);
 66:   /* Remove the values from is2 */
 67:   ISGetIndices(is2,&i2);
 68:   ISGetLocalSize(is2,&n2);
 69:   for (i=0; i<n2; i++) {
 70:     if (i2[i] < imin || i2[i] > imax) continue;
 71:     PetscBTClear(mask,i2[i] - imin);
 72:   }
 73:   ISRestoreIndices(is2,&i2);

 75:   /* Count the number in the difference */
 76:   nout = 0;
 77:   for (i=0; i<imax-imin+1; i++) {
 78:     if (PetscBTLookup(mask,i)) nout++;
 79:   }

 81:   /* create the new IS containing the difference */
 82:   PetscMalloc1(nout,&iout);
 83:   nout = 0;
 84:   for (i=0; i<imax-imin+1; i++) {
 85:     if (PetscBTLookup(mask,i)) iout[nout++] = i + imin;
 86:   }
 87:   PetscObjectGetComm((PetscObject)is1,&comm);
 88:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

 90:   PetscBTDestroy(&mask);
 91:   return(0);
 92: }

 96: /*@
 97:    ISSum - Computes the sum (union) of two index sets.

 99:    Only sequential version (at the moment)

101:    Input Parameter:
102: +  is1 - index set to be extended
103: -  is2 - index values to be added

105:    Output Parameter:
106: .   is3 - the sum; this can not be is1 or is2

108:    Notes:
109:    If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;

111:    Both index sets need to be sorted on input.

113:    Level: intermediate

115: .seealso: ISDestroy(), ISView(), ISDifference(), ISExpand()

117:    Concepts: index sets^union
118:    Concepts: IS^union

120: @*/
121: PetscErrorCode  ISSum(IS is1,IS is2,IS *is3)
122: {
123:   MPI_Comm       comm;
124:   PetscBool      f;
125:   PetscMPIInt    size;
126:   const PetscInt *i1,*i2;
127:   PetscInt       n1,n2,n3, p1,p2, *iout;

133:   PetscObjectGetComm((PetscObject)(is1),&comm);
134:   MPI_Comm_size(comm,&size);
135:   if (size>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for uni-processor IS");

137:   ISSorted(is1,&f);
138:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
139:   ISSorted(is2,&f);
140:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");

142:   ISGetLocalSize(is1,&n1);
143:   ISGetLocalSize(is2,&n2);
144:   if (!n2) return(0);
145:   ISGetIndices(is1,&i1);
146:   ISGetIndices(is2,&i2);

148:   p1 = 0; p2 = 0; n3 = 0;
149:   do {
150:     if (p1==n1) { /* cleanup for is2 */ n3 += n2-p2; break;
151:     } else {
152:       while (p2<n2 && i2[p2]<i1[p1]) {
153:         n3++; p2++;
154:       }
155:       if (p2==n2) {
156:         /* cleanup for is1 */
157:         n3 += n1-p1; break;
158:       } else {
159:         if (i2[p2]==i1[p1]) { n3++; p1++; p2++; }
160:       }
161:     }
162:     if (p2==n2) {
163:       /* cleanup for is1 */
164:       n3 += n1-p1; break;
165:     } else {
166:       while (p1<n1 && i1[p1]<i2[p2]) {
167:         n3++; p1++;
168:       }
169:       if (p1==n1) {
170:         /* clean up for is2 */
171:         n3 += n2-p2; break;
172:       } else {
173:         if (i1[p1]==i2[p2]) { n3++; p1++; p2++; }
174:       }
175:     }
176:   } while (p1<n1 || p2<n2);

178:   if (n3==n1) { /* no new elements to be added */
179:     ISRestoreIndices(is1,&i1);
180:     ISRestoreIndices(is2,&i2);
181:     ISDuplicate(is1,is3);
182:     return(0);
183:   }
184:   PetscMalloc1(n3,&iout);

186:   p1 = 0; p2 = 0; n3 = 0;
187:   do {
188:     if (p1==n1) { /* cleanup for is2 */
189:       while (p2<n2) iout[n3++] = i2[p2++];
190:       break;
191:     } else {
192:       while (p2<n2 && i2[p2]<i1[p1]) iout[n3++] = i2[p2++];
193:       if (p2==n2) { /* cleanup for is1 */
194:         while (p1<n1) iout[n3++] = i1[p1++];
195:         break;
196:       } else {
197:         if (i2[p2]==i1[p1]) { iout[n3++] = i1[p1++]; p2++; }
198:       }
199:     }
200:     if (p2==n2) { /* cleanup for is1 */
201:       while (p1<n1) iout[n3++] = i1[p1++];
202:       break;
203:     } else {
204:       while (p1<n1 && i1[p1]<i2[p2]) iout[n3++] = i1[p1++];
205:       if (p1==n1) { /* clean up for is2 */
206:         while (p2<n2) iout[n3++] = i2[p2++];
207:         break;
208:       } else {
209:         if (i1[p1]==i2[p2]) { iout[n3++] = i1[p1++]; p2++; }
210:       }
211:     }
212:   } while (p1<n1 || p2<n2);

214:   ISRestoreIndices(is1,&i1);
215:   ISRestoreIndices(is2,&i2);
216:   ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
217:   return(0);
218: }

222: /*@
223:    ISExpand - Computes the union of two index sets, by concatenating 2 lists and
224:    removing duplicates.

226:    Collective on IS

228:    Input Parameter:
229: +  is1 - first index set
230: -  is2 - index values to be added

232:    Output Parameters:
233: .  isout - is1 + is2 The index set is2 is appended to is1 removing duplicates

235:    Notes:
236:    Negative values are removed from the lists. This requires O(imax-imin)
237:    memory and O(imax-imin) work, where imin and imax are the bounds on the
238:    indices in is1 and is2.

240:    The IS's do not need to be sorted.

242:    Level: intermediate

244: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum()

246:    Concepts: index sets^difference
247:    Concepts: IS^difference

249: @*/
250: PetscErrorCode ISExpand(IS is1,IS is2,IS *isout)
251: {
253:   PetscInt       i,n1,n2,imin,imax,nout,*iout;
254:   const PetscInt *i1,*i2;
255:   PetscBT        mask;
256:   MPI_Comm       comm;


263:   ISGetIndices(is1,&i1);
264:   ISGetLocalSize(is1,&n1);
265:   ISGetIndices(is2,&i2);
266:   ISGetLocalSize(is2,&n2);

268:   /* Create a bit mask array to contain required values */
269:   if (n1 || n2) {
270:     imin = PETSC_MAX_INT;
271:     imax = 0;
272:     for (i=0; i<n1; i++) {
273:       if (i1[i] < 0) continue;
274:       imin = PetscMin(imin,i1[i]);
275:       imax = PetscMax(imax,i1[i]);
276:     }
277:     for (i=0; i<n2; i++) {
278:       if (i2[i] < 0) continue;
279:       imin = PetscMin(imin,i2[i]);
280:       imax = PetscMax(imax,i2[i]);
281:     }
282:   } else imin = imax = 0;

284:   PetscMalloc1((n1+n2),&iout);
285:   nout = 0;
286:   PetscBTCreate(imax-imin,&mask);
287:   /* Put the values from is1 */
288:   for (i=0; i<n1; i++) {
289:     if (i1[i] < 0) continue;
290:     if (!PetscBTLookupSet(mask,i1[i] - imin)) iout[nout++] = i1[i];
291:   }
292:   ISRestoreIndices(is1,&i1);
293:   /* Put the values from is2 */
294:   for (i=0; i<n2; i++) {
295:     if (i2[i] < 0) continue;
296:     if (!PetscBTLookupSet(mask,i2[i] - imin)) iout[nout++] = i2[i];
297:   }
298:   ISRestoreIndices(is2,&i2);

300:   /* create the new IS containing the sum */
301:   PetscObjectGetComm((PetscObject)is1,&comm);
302:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

304:   PetscBTDestroy(&mask);
305:   return(0);
306: }

310: /*@
311:    ISConcatenate - Forms a new IS by locally concatenating the indices from an IS list without reordering.


314:    Collective on comm.

316:    Input Parameter:
317: +  comm    - communicator of the concatenated IS.
318: .  len     - size of islist array (nonnegative)
319: -  islist  - array of index sets



323:    Output Parameters:
324: .  isout   - The concatenated index set; empty, if len == 0.

326:    Notes: The semantics of calling this on comm imply that the comms of the members if islist also contain this rank.

328:    Level: intermediate

330: .seealso: ISDifference(), ISSum(), ISExpand()

332:    Concepts: index sets^concatenation
333:    Concepts: IS^concatenation

335: @*/
336: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
337: {
339:   PetscInt i,n,N;
340:   const PetscInt *iidx;
341:   PetscInt *idx;

345: #if defined(PETSC_USE_DEBUG)
347: #endif
349:   if (!len) {
350:     ISCreateStride(comm, 0,0,0, isout);
351:     return(0);
352:   }
353:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
354:   N = 0;
355:   for (i = 0; i < len; ++i) {
356:     ISGetLocalSize(islist[i], &n);
357:     N   += n;
358:   }
359:   PetscMalloc(sizeof(PetscInt)*N, &idx);
360:   N = 0;
361:   for (i = 0; i < len; ++i) {
362:     ISGetLocalSize(islist[i], &n);
363:     ISGetIndices(islist[i], &iidx);
364:     PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
365:     N   += n;
366:   }
367:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
368:   return(0);
369: }

371: /*@
372:    ISListToPair     -    convert an IS list to a pair of ISs of equal length defining an equivalent integer multimap.
373:                         Each IS on the input list is assigned an integer j so that all of the indices of that IS are
374:                         mapped to j.


377:   Collective on comm.

379:   Input arguments:
380: + comm    -  MPI_Comm
381: . listlen -  IS list length
382: - islist  -  IS list

384:   Output arguments:
385: + xis -  domain IS
386: - yis -  range  IS

388:   Level: advanced

390:   Notes:
391:   The global integers assigned to the ISs of the local input list might not correspond to the
392:   local numbers of the ISs on that list, but the two *orderings* are the same: the global
393:   integers assigned to the ISs on the local list form a strictly increasing sequence.

395:   The ISs on the input list can belong to subcommunicators of comm, and the subcommunicators
396:   on the input IS list are assumed to be in a "deadlock-free" order.

398:   Local lists of PetscObjects (or their subcommes) on a comm are "deadlock-free" if subcomm1
399:   preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
400:   Equivalently, the local numbers of the subcomms on each local list are drawn from some global
401:   numbering. This is ensured, for example, by ISPairToList().

403: .seealso ISPairToList()
404: @*/
405: #undef  __FUNCT__
407: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
408: {
410:   PetscInt       ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
411:   const PetscInt *indsi;

414:   PetscMalloc1(listlen, &colors);
415:   PetscObjectsGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
416:   len  = 0;
417:   for (i = 0; i < listlen; ++i) {
418:     ISGetLocalSize(islist[i], &leni);
419:     len += leni;
420:   }
421:   PetscMalloc1(len, &xinds);
422:   PetscMalloc1(len, &yinds);
423:   k    = 0;
424:   for (i = 0; i < listlen; ++i) {
425:     ISGetLocalSize(islist[i], &leni);
426:     ISGetIndices(islist[i],&indsi);
427:     for (j = 0; j < leni; ++j) {
428:       xinds[k] = indsi[j];
429:       yinds[k] = colors[i];
430:       ++k;
431:     }
432:   }
433:   PetscFree(colors);
434:   ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
435:   ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
436:   return(0);
437: }


440: /*@
441:    ISPairToList   -   convert an IS pair encoding an integer map to a list of ISs.
442:                      Each IS on the output list contains the preimage for each index on the second input IS.
443:                      The ISs on the output list are constructed on the subcommunicators of the input IS pair.
444:                      Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
445:                      exactly the ranks that assign some indices i to j.  This is essentially the inverse of
446:                      ISListToPair().

448:   Collective on indis.

450:   Input arguments:
451: + xis -  domain IS
452: - yis -  range IS

454:   Output arguments:
455: + listlen -  length of islist
456: - islist  -  list of ISs breaking up indis by color

458:   Note:
459: + xis and yis must be of the same length and have congruent communicators.
460: - The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).

462:   Level: advanced

464: .seealso ISListToPair()
465:  @*/
466: #undef  __FUNCT__
468: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
469: {
471:   IS             indis = xis, coloris = yis;
472:   PetscInt       *inds, *colors, llen, ilen, lstart, lend, lcount,l;
473:   PetscMPIInt    rank, size, llow, lhigh, low, high,color,subsize;
474:   const PetscInt *ccolors, *cinds;
475:   MPI_Comm       comm, subcomm;

483:   PetscObjectGetComm((PetscObject)xis,&comm);
484:   MPI_Comm_rank(comm, &rank);
485:   MPI_Comm_rank(comm, &size);
486:   /* Extract, copy and sort the local indices and colors on the color. */
487:   ISGetLocalSize(coloris, &llen);
488:   ISGetLocalSize(indis,   &ilen);
489:   if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
490:   ISGetIndices(coloris, &ccolors);
491:   ISGetIndices(indis, &cinds);
492:   PetscMalloc2(ilen,&inds,llen,&colors);
493:   PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
494:   PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
495:   PetscSortIntWithArray(llen, colors, inds);
496:   /* Determine the global extent of colors. */
497:   llow   = 0; lhigh  = -1;
498:   lstart = 0; lcount = 0;
499:   while (lstart < llen) {
500:     lend = lstart+1;
501:     while (lend < llen && colors[lend] == colors[lstart]) ++lend;
502:     llow  = PetscMin(llow,colors[lstart]);
503:     lhigh = PetscMax(lhigh,colors[lstart]);
504:     ++lcount;
505:   }
506:   MPI_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
507:   MPI_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
508:   *listlen = 0;
509:   if (low <= high) {
510:     if (lcount > 0) {
511:       *listlen = lcount;
512:       if (!*islist) {
513:         PetscMalloc(sizeof(IS)*lcount, islist);
514:       }
515:     }
516:     /*
517:      Traverse all possible global colors, and participate in the subcommunicators
518:      for the locally-supported colors.
519:      */
520:     lcount = 0;
521:     lstart = 0; lend = 0;
522:     for (l = low; l <= high; ++l) {
523:       /*
524:        Find the range of indices with the same color, which is not smaller than l.
525:        Observe that, since colors is sorted, and is a subsequence of [low,high],
526:        as soon as we find a new color, it is >= l.
527:        */
528:       if (lstart < llen) {
529:         /* The start of the next locally-owned color is identified.  Now look for the end. */
530:         if (lstart == lend) {
531:           lend = lstart+1;
532:           while (lend < llen && colors[lend] == colors[lstart]) ++lend;
533:         }
534:         /* Now check whether the identified color segment matches l. */
535:         if (colors[lstart] < l) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Locally owned color %D at location %D is < than the next global color %D", colors[lstart], lcount, l);
536:       }
537:       color = (PetscMPIInt)(colors[lstart] == l);
538:       /* Check whether a proper subcommunicator exists. */
539:       MPI_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);

541:       if (subsize == 1) subcomm = PETSC_COMM_SELF;
542:       else if (subsize == size) subcomm = comm;
543:       else {
544:         /* a proper communicator is necessary, so we create it. */
545:         MPI_Comm_split(comm, color, rank, &subcomm);
546:       }
547:       if (colors[lstart] == l) {
548:         /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
549:         ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
550:         /* Position lstart at the beginning of the next local color. */
551:         lstart = lend;
552:         /* Increment the counter of the local colors split off into an IS. */
553:         ++lcount;
554:       }
555:       if (subsize > 0 && subsize < size) {
556:         /*
557:          Irrespective of color, destroy the split off subcomm:
558:          a subcomm used in the IS creation above is duplicated
559:          into a proper PETSc comm.
560:          */
561:         MPI_Comm_free(&subcomm);
562:       }
563:     } /* for (l = low; l < high; ++l) */
564:   } /* if (low <= high) */
565:   PetscFree2(inds,colors);
566:   return(0);
567: }


570: /*@
571:    ISEmbed   -   embed IS a into IS b by finding the locations in b that have the same indices as in a.
572:                  If c is the IS of these locations, we have a = b*c, regarded as a composition of the 
573:                  corresponding ISLocalToGlobalMaps.

575:   Not collective.

577:   Input arguments:
578: + a    -  IS to embed
579: . b    -  IS to embed into
580: - drop -  flag indicating whether to drop a's indices that are not in b.

582:   Output arguments:
583: . c    -  local embedding indices

585:   Note:
586:   If some of a's global indices are not among b's indices the embedding is impossible.  The local indices of a
587:   corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop).  In the former
588:   case the size of c is that same as that of a, in the latter case c's size may be smaller.

590:   The resulting IS is sequential, since the index substition it encodes is purely local.

592:   Level: advanced

594: .seealso ISLocalToGlobalMapping
595:  @*/
596: #undef  __FUNCT__
598: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
599: {
600:   PetscErrorCode             ierr;
601:   ISLocalToGlobalMapping     ltog;
602:   ISGlobalToLocalMappingType gtoltype = IS_GTOLM_DROP;
603:   PetscInt                   alen, clen, *cindices, *cindices2;
604:   const PetscInt             *aindices;

610:   ISLocalToGlobalMappingCreateIS(b, &ltog);
611:   ISGetLocalSize(a, &alen);
612:   ISGetIndices(a, &aindices);
613:   PetscMalloc1(alen, &cindices);
614:   if (!drop) gtoltype = IS_GTOLM_MASK;
615:   ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
616:   if (clen != alen) {
617:     cindices2 = cindices;
618:     PetscMalloc1(clen, &cindices);
619:     PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
620:     PetscFree(cindices2);
621:   }
622:   ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c);
623:   return(0);
624: }