Actual source code: isdiff.c

petsc-master 2020-05-29
Report Typos and Errors

  2:  #include <petsc/private/isimpl.h>
  3:  #include <petsc/private/sectionimpl.h>
  4:  #include <petscbt.h>

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

  9:    Collective on IS

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

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

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

 23:    If is2 is NULL, the result is the same as for an empty IS, i.e., a duplicate of is1.

 25:    Level: intermediate

 27: .seealso: ISDestroy(), ISView(), ISSum(), ISExpand()
 28: @*/
 29: PetscErrorCode  ISDifference(IS is1,IS is2,IS *isout)
 30: {
 32:   PetscInt       i,n1,n2,imin,imax,nout,*iout;
 33:   const PetscInt *i1,*i2;
 34:   PetscBT        mask;
 35:   MPI_Comm       comm;

 40:   if (!is2) {
 41:     ISDuplicate(is1, isout);
 42:     return(0);
 43:   }

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

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

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

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

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

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

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

 98:    Only sequential version (at the moment)

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

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

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

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

112:    Level: intermediate

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


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

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

134:   ISSorted(is1,&f);
135:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
136:   ISSorted(is2,&f);
137:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");

139:   ISGetLocalSize(is1,&n1);
140:   ISGetLocalSize(is2,&n2);
141:   if (!n2) {
142:     ISDuplicate(is1,is3);
143:     return(0);
144:   }
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: }

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

224:    Collective on IS

226:    Input Parameter:
227: +  is1 - first index set
228: -  is2 - index values to be added

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

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

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

240:    Level: intermediate

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


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


259:   if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
260:   if (!is1) {ISDuplicate(is2, isout);return(0);}
261:   if (!is2) {ISDuplicate(is1, isout);return(0);}
262:   ISGetIndices(is1,&i1);
263:   ISGetLocalSize(is1,&n1);
264:   ISGetIndices(is2,&i2);
265:   ISGetLocalSize(is2,&n2);

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

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

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

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

307: /*@
308:    ISIntersect - Computes the intersection of two index sets, by sorting and comparing.

310:    Collective on IS

312:    Input Parameter:
313: +  is1 - first index set
314: -  is2 - second index set

316:    Output Parameters:
317: .  isout - the sorted intersection of is1 and is2

319:    Notes:
320:    Negative values are removed from the lists. This requires O(min(is1,is2))
321:    memory and O(max(is1,is2)log(max(is1,is2))) work

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

325:    Level: intermediate

327: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
328: @*/
329: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
330: {
332:   PetscInt       i,n1,n2,nout,*iout;
333:   const PetscInt *i1,*i2;
334:   IS             is1sorted = NULL, is2sorted = NULL;
335:   PetscBool      sorted, lsorted;
336:   MPI_Comm       comm;

343:   PetscObjectGetComm((PetscObject)is1,&comm);

345:   ISGetLocalSize(is1,&n1);
346:   ISGetLocalSize(is2,&n2);
347:   if (n1 < n2) {
348:     IS       tempis = is1;
349:     PetscInt ntemp = n1;

351:     is1 = is2;
352:     is2 = tempis;
353:     n1  = n2;
354:     n2  = ntemp;
355:   }
356:   ISSorted(is1,&lsorted);
357:   MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
358:   if (!sorted) {
359:     ISDuplicate(is1,&is1sorted);
360:     ISSort(is1sorted);
361:     ISGetIndices(is1sorted,&i1);
362:   } else {
363:     is1sorted = is1;
364:     PetscObjectReference((PetscObject)is1);
365:     ISGetIndices(is1,&i1);
366:   }
367:   ISSorted(is2,&lsorted);
368:   MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
369:   if (!sorted) {
370:     ISDuplicate(is2,&is2sorted);
371:     ISSort(is2sorted);
372:     ISGetIndices(is2sorted,&i2);
373:   } else {
374:     is2sorted = is2;
375:     PetscObjectReference((PetscObject)is2);
376:     ISGetIndices(is2,&i2);
377:   }

379:   PetscMalloc1(n2,&iout);

381:   for (i = 0, nout = 0; i < n2; i++) {
382:     PetscInt key = i2[i];
383:     PetscInt loc;

385:     ISLocate(is1sorted,key,&loc);
386:     if (loc >= 0) {
387:       if (!nout || iout[nout-1] < key) {
388:         iout[nout++] = key;
389:       }
390:     }
391:   }
392:   PetscRealloc(nout*sizeof(PetscInt),&iout);

394:   /* create the new IS containing the sum */
395:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

397:   ISRestoreIndices(is2sorted,&i2);
398:   ISDestroy(&is2sorted);
399:   ISRestoreIndices(is1sorted,&i1);
400:   ISDestroy(&is1sorted);
401:   return(0);
402: }

404: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
405: {

409:   *isect = NULL;
410:   if (is2 && is1) {
411:     char           composeStr[33] = {0};
412:     PetscObjectId  is2id;

414:     PetscObjectGetId((PetscObject)is2,&is2id);
415:     PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%x",is2id);
416:     PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
417:     if (*isect == NULL) {
418:       ISIntersect(is1, is2, isect);
419:       PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
420:     } else {
421:       PetscObjectReference((PetscObject) *isect);
422:     }
423:   }
424:   return(0);
425: }

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


431:    Collective.

433:    Input Parameter:
434: +  comm    - communicator of the concatenated IS.
435: .  len     - size of islist array (nonnegative)
436: -  islist  - array of index sets

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

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

444:    Level: intermediate

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


449: @*/
450: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
451: {
453:   PetscInt i,n,N;
454:   const PetscInt *iidx;
455:   PetscInt *idx;

459:   if (PetscDefined(USE_DEBUG)) {
461:   }
463:   if (!len) {
464:     ISCreateStride(comm, 0,0,0, isout);
465:     return(0);
466:   }
467:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
468:   N = 0;
469:   for (i = 0; i < len; ++i) {
470:     if (islist[i]) {
471:       ISGetLocalSize(islist[i], &n);
472:       N   += n;
473:     }
474:   }
475:   PetscMalloc1(N, &idx);
476:   N = 0;
477:   for (i = 0; i < len; ++i) {
478:     if (islist[i]) {
479:       ISGetLocalSize(islist[i], &n);
480:       ISGetIndices(islist[i], &iidx);
481:       PetscArraycpy(idx+N,iidx, n);
482:       ISRestoreIndices(islist[i], &iidx);
483:       N   += n;
484:     }
485:   }
486:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
487:   return(0);
488: }

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


496:   Collective.

498:   Input arguments:
499: + comm    -  MPI_Comm
500: . listlen -  IS list length
501: - islist  -  IS list

503:   Output arguments:
504: + xis -  domain IS
505: - yis -  range  IS

507:   Level: advanced

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

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

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

522: .seealso ISPairToList()
523: @*/
524: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
525: {
527:   PetscInt       ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
528:   const PetscInt *indsi;

531:   PetscMalloc1(listlen, &colors);
532:   PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
533:   len  = 0;
534:   for (i = 0; i < listlen; ++i) {
535:     ISGetLocalSize(islist[i], &leni);
536:     len += leni;
537:   }
538:   PetscMalloc1(len, &xinds);
539:   PetscMalloc1(len, &yinds);
540:   k    = 0;
541:   for (i = 0; i < listlen; ++i) {
542:     ISGetLocalSize(islist[i], &leni);
543:     ISGetIndices(islist[i],&indsi);
544:     for (j = 0; j < leni; ++j) {
545:       xinds[k] = indsi[j];
546:       yinds[k] = colors[i];
547:       ++k;
548:     }
549:   }
550:   PetscFree(colors);
551:   ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
552:   ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
553:   return(0);
554: }


557: /*@
558:    ISPairToList   -   convert an IS pair encoding an integer map to a list of ISs.
559:                      Each IS on the output list contains the preimage for each index on the second input IS.
560:                      The ISs on the output list are constructed on the subcommunicators of the input IS pair.
561:                      Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
562:                      exactly the ranks that assign some indices i to j.  This is essentially the inverse of
563:                      ISListToPair().

565:   Collective on indis.

567:   Input arguments:
568: + xis -  domain IS
569: - yis -  range IS

571:   Output arguments:
572: + listlen -  length of islist
573: - islist  -  list of ISs breaking up indis by color

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

579:   Level: advanced

581: .seealso ISListToPair()
582:  @*/
583: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
584: {
586:   IS             indis = xis, coloris = yis;
587:   PetscInt       *inds, *colors, llen, ilen, lstart, lend, lcount,l;
588:   PetscMPIInt    rank, size, llow, lhigh, low, high,color,subsize;
589:   const PetscInt *ccolors, *cinds;
590:   MPI_Comm       comm, subcomm;

598:   PetscObjectGetComm((PetscObject)xis,&comm);
599:   MPI_Comm_rank(comm, &rank);
600:   MPI_Comm_rank(comm, &size);
601:   /* Extract, copy and sort the local indices and colors on the color. */
602:   ISGetLocalSize(coloris, &llen);
603:   ISGetLocalSize(indis,   &ilen);
604:   if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
605:   ISGetIndices(coloris, &ccolors);
606:   ISGetIndices(indis, &cinds);
607:   PetscMalloc2(ilen,&inds,llen,&colors);
608:   PetscArraycpy(inds,cinds,ilen);
609:   PetscArraycpy(colors,ccolors,llen);
610:   PetscSortIntWithArray(llen, colors, inds);
611:   /* Determine the global extent of colors. */
612:   llow   = 0; lhigh  = -1;
613:   lstart = 0; lcount = 0;
614:   while (lstart < llen) {
615:     lend = lstart+1;
616:     while (lend < llen && colors[lend] == colors[lstart]) ++lend;
617:     llow  = PetscMin(llow,colors[lstart]);
618:     lhigh = PetscMax(lhigh,colors[lstart]);
619:     ++lcount;
620:   }
621:   MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
622:   MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
623:   *listlen = 0;
624:   if (low <= high) {
625:     if (lcount > 0) {
626:       *listlen = lcount;
627:       if (!*islist) {
628:         PetscMalloc1(lcount, islist);
629:       }
630:     }
631:     /*
632:      Traverse all possible global colors, and participate in the subcommunicators
633:      for the locally-supported colors.
634:      */
635:     lcount = 0;
636:     lstart = 0; lend = 0;
637:     for (l = low; l <= high; ++l) {
638:       /*
639:        Find the range of indices with the same color, which is not smaller than l.
640:        Observe that, since colors is sorted, and is a subsequence of [low,high],
641:        as soon as we find a new color, it is >= l.
642:        */
643:       if (lstart < llen) {
644:         /* The start of the next locally-owned color is identified.  Now look for the end. */
645:         if (lstart == lend) {
646:           lend = lstart+1;
647:           while (lend < llen && colors[lend] == colors[lstart]) ++lend;
648:         }
649:         /* Now check whether the identified color segment matches l. */
650:         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);
651:       }
652:       color = (PetscMPIInt)(colors[lstart] == l);
653:       /* Check whether a proper subcommunicator exists. */
654:       MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);

656:       if (subsize == 1) subcomm = PETSC_COMM_SELF;
657:       else if (subsize == size) subcomm = comm;
658:       else {
659:         /* a proper communicator is necessary, so we create it. */
660:         MPI_Comm_split(comm, color, rank, &subcomm);
661:       }
662:       if (colors[lstart] == l) {
663:         /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
664:         ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
665:         /* Position lstart at the beginning of the next local color. */
666:         lstart = lend;
667:         /* Increment the counter of the local colors split off into an IS. */
668:         ++lcount;
669:       }
670:       if (subsize > 0 && subsize < size) {
671:         /*
672:          Irrespective of color, destroy the split off subcomm:
673:          a subcomm used in the IS creation above is duplicated
674:          into a proper PETSc comm.
675:          */
676:         MPI_Comm_free(&subcomm);
677:       }
678:     } /* for (l = low; l < high; ++l) */
679:   } /* if (low <= high) */
680:   PetscFree2(inds,colors);
681:   return(0);
682: }


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

690:   Not collective.

692:   Input arguments:
693: + a    -  IS to embed
694: . b    -  IS to embed into
695: - drop -  flag indicating whether to drop a's indices that are not in b.

697:   Output arguments:
698: . c    -  local embedding indices

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

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

707:   Level: advanced

709: .seealso ISLocalToGlobalMapping
710:  @*/
711: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
712: {
713:   PetscErrorCode             ierr;
714:   ISLocalToGlobalMapping     ltog;
715:   ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
716:   PetscInt                   alen, clen, *cindices, *cindices2;
717:   const PetscInt             *aindices;

723:   ISLocalToGlobalMappingCreateIS(b, &ltog);
724:   ISGetLocalSize(a, &alen);
725:   ISGetIndices(a, &aindices);
726:   PetscMalloc1(alen, &cindices);
727:   if (!drop) gtoltype = IS_GTOLM_MASK;
728:   ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
729:   ISLocalToGlobalMappingDestroy(&ltog);
730:   if (clen != alen) {
731:     cindices2 = cindices;
732:     PetscMalloc1(clen, &cindices);
733:     PetscArraycpy(cindices,cindices2,clen);
734:     PetscFree(cindices2);
735:   }
736:   ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
737:   return(0);
738: }


741: /*@
742:   ISSortPermutation  -  calculate the permutation of the indices into a nondecreasing order.

744:   Not collective.

746:   Input arguments:
747: + f      -  IS to sort
748: - always -  build the permutation even when f's indices are nondecreasing.

750:   Output argument:
751: . h    -  permutation or NULL, if f is nondecreasing and always == PETSC_FALSE.


754:   Note: Indices in f are unchanged. f[h[i]] is the i-th smallest f index.
755:         If always == PETSC_FALSE, an extra check is peformed to see whether
756:         the f indices are nondecreasing. h is built on PETSC_COMM_SELF, since
757:         the permutation has a local meaning only.

759:   Level: advanced

761: .seealso ISLocalToGlobalMapping, ISSort()
762:  @*/
763: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
764: {
765:   PetscErrorCode  ierr;
766:   const PetscInt  *findices;
767:   PetscInt        fsize,*hindices,i;
768:   PetscBool       isincreasing;

773:   ISGetLocalSize(f,&fsize);
774:   ISGetIndices(f,&findices);
775:   *h = NULL;
776:   if (!always) {
777:     isincreasing = PETSC_TRUE;
778:     for (i = 1; i < fsize; ++i) {
779:       if (findices[i] <= findices[i-1]) {
780:         isincreasing = PETSC_FALSE;
781:         break;
782:       }
783:     }
784:     if (isincreasing) {
785:       ISRestoreIndices(f,&findices);
786:       return(0);
787:     }
788:   }
789:   PetscMalloc1(fsize,&hindices);
790:   for (i = 0; i < fsize; ++i) hindices[i] = i;
791:   PetscSortIntWithPermutation(fsize,findices,hindices);
792:   ISRestoreIndices(f,&findices);
793:   ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
794:   ISSetInfo(*h,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,PETSC_TRUE);
795:   return(0);
796: }