Actual source code: isdiff.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

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

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

  8:    Collective on IS

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

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

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

 22:    Level: intermediate

 24:    Concepts: index sets^difference
 25:    Concepts: IS^difference

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

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


 43:   ISGetIndices(is1,&i1);
 44:   ISGetLocalSize(is1,&n1);

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

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

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

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

 88:   PetscBTDestroy(&mask);
 89:   return(0);
 90: }

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

 95:    Only sequential version (at the moment)

 97:    Input Parameter:
 98: +  is1 - index set to be extended
 99: -  is2 - index values to be added

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

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

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

109:    Level: intermediate

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

113:    Concepts: index sets^union
114:    Concepts: IS^union

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

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

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

138:   ISGetLocalSize(is1,&n1);
139:   ISGetLocalSize(is2,&n2);
140:   if (!n2) {
141:     ISDuplicate(is1,is3);
142:     return(0);
143:   }
144:   ISGetIndices(is1,&i1);
145:   ISGetIndices(is2,&i2);

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

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

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

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

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

223:    Collective on IS

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

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

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

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

239:    Level: intermediate

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

243:    Concepts: index sets^difference
244:    Concepts: IS^difference

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


260:   if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
261:   if (!is1) {ISDuplicate(is2, isout);return(0);}
262:   if (!is2) {ISDuplicate(is1, isout);return(0);}
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: }

308: /*@
309:    ISIntersect - Computes the union of two index sets, by concatenating 2 lists, sorting,
310:    and finding duplicates.

312:    Collective on IS

314:    Input Parameter:
315: +  is1 - first index set
316: -  is2 - second index set

318:    Output Parameters:
319: .  isout - the sorted intersection of is1 and is2

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

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

327:    Level: intermediate

329: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()

331:    Concepts: index sets^intersection
332:    Concepts: IS^intersection

334: @*/
335: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
336: {
338:   PetscInt       i,n1,n2,nout,*iout;
339:   const PetscInt *i1,*i2;
340:   IS             is1sorted = NULL, is2sorted = NULL;
341:   PetscBool      sorted;
342:   MPI_Comm       comm;


349:   ISGetLocalSize(is1,&n1);
350:   ISGetLocalSize(is2,&n2);
351:   if (n1 < n2) {
352:     IS tempis = is1;
353:     int ntemp = n1;

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

381:   PetscMalloc1(n2,&iout);

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

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

396:   /* create the new IS containing the sum */
397:   PetscObjectGetComm((PetscObject)is1,&comm);
398:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

400:   ISRestoreIndices(is2sorted,&i2);
401:   ISDestroy(&is2sorted);
402:   ISRestoreIndices(is1sorted,&i1);
403:   ISDestroy(&is1sorted);
404:   return(0);
405: }

407: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
408: {

412:   *isect = NULL;
413:   if (is2 && is1) {
414:     char           composeStr[33] = {0};
415:     PetscObjectId  is2id;

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

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


434:    Collective on comm.

436:    Input Parameter:
437: +  comm    - communicator of the concatenated IS.
438: .  len     - size of islist array (nonnegative)
439: -  islist  - array of index sets

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

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

447:    Level: intermediate

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

451:    Concepts: index sets^concatenation
452:    Concepts: IS^concatenation

454: @*/
455: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
456: {
458:   PetscInt i,n,N;
459:   const PetscInt *iidx;
460:   PetscInt *idx;

464: #if defined(PETSC_USE_DEBUG)
466: #endif
468:   if (!len) {
469:     ISCreateStride(comm, 0,0,0, isout);
470:     return(0);
471:   }
472:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
473:   N = 0;
474:   for (i = 0; i < len; ++i) {
475:     if (islist[i]) {
476:       ISGetLocalSize(islist[i], &n);
477:       N   += n;
478:     }
479:   }
480:   PetscMalloc1(N, &idx);
481:   N = 0;
482:   for (i = 0; i < len; ++i) {
483:     if (islist[i]) {
484:       ISGetLocalSize(islist[i], &n);
485:       ISGetIndices(islist[i], &iidx);
486:       PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
487:       ISRestoreIndices(islist[i], &iidx);
488:       N   += n;
489:     }
490:   }
491:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
492:   return(0);
493: }

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


501:   Collective on comm.

503:   Input arguments:
504: + comm    -  MPI_Comm
505: . listlen -  IS list length
506: - islist  -  IS list

508:   Output arguments:
509: + xis -  domain IS
510: - yis -  range  IS

512:   Level: advanced

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

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

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

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

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


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

570:   Collective on indis.

572:   Input arguments:
573: + xis -  domain IS
574: - yis -  range IS

576:   Output arguments:
577: + listlen -  length of islist
578: - islist  -  list of ISs breaking up indis by color

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

584:   Level: advanced

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

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

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


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

695:   Not collective.

697:   Input arguments:
698: + a    -  IS to embed
699: . b    -  IS to embed into
700: - drop -  flag indicating whether to drop a's indices that are not in b.

702:   Output arguments:
703: . c    -  local embedding indices

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

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

712:   Level: advanced

714: .seealso ISLocalToGlobalMapping
715:  @*/
716: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
717: {
718:   PetscErrorCode             ierr;
719:   ISLocalToGlobalMapping     ltog;
720:   ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
721:   PetscInt                   alen, clen, *cindices, *cindices2;
722:   const PetscInt             *aindices;

728:   ISLocalToGlobalMappingCreateIS(b, &ltog);
729:   ISGetLocalSize(a, &alen);
730:   ISGetIndices(a, &aindices);
731:   PetscMalloc1(alen, &cindices);
732:   if (!drop) gtoltype = IS_GTOLM_MASK;
733:   ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
734:   ISLocalToGlobalMappingDestroy(&ltog);
735:   if (clen != alen) {
736:     cindices2 = cindices;
737:     PetscMalloc1(clen, &cindices);
738:     PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
739:     PetscFree(cindices2);
740:   }
741:   ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
742:   return(0);
743: }


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

749:   Not collective.

751:   Input arguments:
752: + f      -  IS to sort
753: - always -  build the permutation even when f's indices are nondecreasin.

755:   Output argument:
756: . h    -  permutation or NULL, if f is nondecreasing and always == PETSC_TRUE.


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

764:   Level: advanced

766: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
767:  @*/
768: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
769: {
770:   PetscErrorCode  ierr;
771:   const PetscInt  *findices;
772:   PetscInt        fsize,*hindices,i;
773:   PetscBool       isincreasing;

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