Actual source code: isdiff.c

petsc-3.8.3 2017-12-09
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) return(0);
141:   ISGetIndices(is1,&i1);
142:   ISGetIndices(is2,&i2);

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

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

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

210:   ISRestoreIndices(is1,&i1);
211:   ISRestoreIndices(is2,&i2);
212:   ISCreateGeneral(comm,n3,iout,PETSC_OWN_POINTER,is3);
213:   return(0);
214: }

216: /*@
217:    ISExpand - Computes the union of two index sets, by concatenating 2 lists and
218:    removing duplicates.

220:    Collective on IS

222:    Input Parameter:
223: +  is1 - first index set
224: -  is2 - index values to be added

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

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

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

236:    Level: intermediate

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

240:    Concepts: index sets^difference
241:    Concepts: IS^difference

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


257:   ISGetIndices(is1,&i1);
258:   ISGetLocalSize(is1,&n1);
259:   ISGetIndices(is2,&i2);
260:   ISGetLocalSize(is2,&n2);

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

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

294:   /* create the new IS containing the sum */
295:   PetscObjectGetComm((PetscObject)is1,&comm);
296:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

298:   PetscBTDestroy(&mask);
299:   return(0);
300: }

302: /*@
303:    ISIntersect - Computes the union of two index sets, by concatenating 2 lists, sorting,
304:    and finding duplicates.

306:    Collective on IS

308:    Input Parameter:
309: +  is1 - first index set
310: -  is2 - second index set

312:    Output Parameters:
313: .  isout - the sorted intersection of is1 and is2

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

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

321:    Level: intermediate

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

325:    Concepts: index sets^intersection
326:    Concepts: IS^intersection

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;
336:   MPI_Comm       comm;


343:   ISGetLocalSize(is1,&n1);
344:   ISGetLocalSize(is2,&n2);
345:   if (n1 < n2) {
346:     IS tempis = is1;
347:     int ntemp = n1;

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

375:   PetscMalloc1(n2,&iout);

377:   for (i = 0, nout = 0; i < n2; i++) {
378:     PetscInt key = i2[i];
379:     PetscInt loc;

381:     ISLocate(is1sorted,key,&loc);
382:     if (loc >= 0) {
383:       if (!nout || iout[nout-1] < key) {
384:         iout[nout++] = key;
385:       }
386:     }
387:   }
388:   PetscRealloc(sizeof (PetscInt) * (size_t) nout,&iout);

390:   /* create the new IS containing the sum */
391:   PetscObjectGetComm((PetscObject)is1,&comm);
392:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

394:   ISRestoreIndices(is2sorted,&i2);
395:   ISDestroy(&is2sorted);
396:   ISRestoreIndices(is1sorted,&i1);
397:   ISDestroy(&is1sorted);
398:   return(0);
399: }

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


405:    Collective on comm.

407:    Input Parameter:
408: +  comm    - communicator of the concatenated IS.
409: .  len     - size of islist array (nonnegative)
410: -  islist  - array of index sets



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

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

419:    Level: intermediate

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

423:    Concepts: index sets^concatenation
424:    Concepts: IS^concatenation

426: @*/
427: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
428: {
430:   PetscInt i,n,N;
431:   const PetscInt *iidx;
432:   PetscInt *idx;

436: #if defined(PETSC_USE_DEBUG)
438: #endif
440:   if (!len) {
441:     ISCreateStride(comm, 0,0,0, isout);
442:     return(0);
443:   }
444:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
445:   N = 0;
446:   for (i = 0; i < len; ++i) {
447:     ISGetLocalSize(islist[i], &n);
448:     N   += n;
449:   }
450:   PetscMalloc1(N, &idx);
451:   N = 0;
452:   for (i = 0; i < len; ++i) {
453:     ISGetLocalSize(islist[i], &n);
454:     ISGetIndices(islist[i], &iidx);
455:     PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
456:     ISRestoreIndices(islist[i], &iidx);
457:     N   += n;
458:   }
459:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
460:   return(0);
461: }

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


469:   Collective on comm.

471:   Input arguments:
472: + comm    -  MPI_Comm
473: . listlen -  IS list length
474: - islist  -  IS list

476:   Output arguments:
477: + xis -  domain IS
478: - yis -  range  IS

480:   Level: advanced

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

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

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

495: .seealso ISPairToList()
496: @*/
497: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
498: {
500:   PetscInt       ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
501:   const PetscInt *indsi;

504:   PetscMalloc1(listlen, &colors);
505:   PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject*)islist,&ncolors, colors);
506:   len  = 0;
507:   for (i = 0; i < listlen; ++i) {
508:     ISGetLocalSize(islist[i], &leni);
509:     len += leni;
510:   }
511:   PetscMalloc1(len, &xinds);
512:   PetscMalloc1(len, &yinds);
513:   k    = 0;
514:   for (i = 0; i < listlen; ++i) {
515:     ISGetLocalSize(islist[i], &leni);
516:     ISGetIndices(islist[i],&indsi);
517:     for (j = 0; j < leni; ++j) {
518:       xinds[k] = indsi[j];
519:       yinds[k] = colors[i];
520:       ++k;
521:     }
522:   }
523:   PetscFree(colors);
524:   ISCreateGeneral(comm,len,xinds,PETSC_OWN_POINTER,xis);
525:   ISCreateGeneral(comm,len,yinds,PETSC_OWN_POINTER,yis);
526:   return(0);
527: }


530: /*@
531:    ISPairToList   -   convert an IS pair encoding an integer map to a list of ISs.
532:                      Each IS on the output list contains the preimage for each index on the second input IS.
533:                      The ISs on the output list are constructed on the subcommunicators of the input IS pair.
534:                      Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
535:                      exactly the ranks that assign some indices i to j.  This is essentially the inverse of
536:                      ISListToPair().

538:   Collective on indis.

540:   Input arguments:
541: + xis -  domain IS
542: - yis -  range IS

544:   Output arguments:
545: + listlen -  length of islist
546: - islist  -  list of ISs breaking up indis by color

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

552:   Level: advanced

554: .seealso ISListToPair()
555:  @*/
556: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
557: {
559:   IS             indis = xis, coloris = yis;
560:   PetscInt       *inds, *colors, llen, ilen, lstart, lend, lcount,l;
561:   PetscMPIInt    rank, size, llow, lhigh, low, high,color,subsize;
562:   const PetscInt *ccolors, *cinds;
563:   MPI_Comm       comm, subcomm;

571:   PetscObjectGetComm((PetscObject)xis,&comm);
572:   MPI_Comm_rank(comm, &rank);
573:   MPI_Comm_rank(comm, &size);
574:   /* Extract, copy and sort the local indices and colors on the color. */
575:   ISGetLocalSize(coloris, &llen);
576:   ISGetLocalSize(indis,   &ilen);
577:   if (llen != ilen) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", ilen, llen);
578:   ISGetIndices(coloris, &ccolors);
579:   ISGetIndices(indis, &cinds);
580:   PetscMalloc2(ilen,&inds,llen,&colors);
581:   PetscMemcpy(inds,cinds,ilen*sizeof(PetscInt));
582:   PetscMemcpy(colors,ccolors,llen*sizeof(PetscInt));
583:   PetscSortIntWithArray(llen, colors, inds);
584:   /* Determine the global extent of colors. */
585:   llow   = 0; lhigh  = -1;
586:   lstart = 0; lcount = 0;
587:   while (lstart < llen) {
588:     lend = lstart+1;
589:     while (lend < llen && colors[lend] == colors[lstart]) ++lend;
590:     llow  = PetscMin(llow,colors[lstart]);
591:     lhigh = PetscMax(lhigh,colors[lstart]);
592:     ++lcount;
593:   }
594:   MPIU_Allreduce(&llow,&low,1,MPI_INT,MPI_MIN,comm);
595:   MPIU_Allreduce(&lhigh,&high,1,MPI_INT,MPI_MAX,comm);
596:   *listlen = 0;
597:   if (low <= high) {
598:     if (lcount > 0) {
599:       *listlen = lcount;
600:       if (!*islist) {
601:         PetscMalloc1(lcount, islist);
602:       }
603:     }
604:     /*
605:      Traverse all possible global colors, and participate in the subcommunicators
606:      for the locally-supported colors.
607:      */
608:     lcount = 0;
609:     lstart = 0; lend = 0;
610:     for (l = low; l <= high; ++l) {
611:       /*
612:        Find the range of indices with the same color, which is not smaller than l.
613:        Observe that, since colors is sorted, and is a subsequence of [low,high],
614:        as soon as we find a new color, it is >= l.
615:        */
616:       if (lstart < llen) {
617:         /* The start of the next locally-owned color is identified.  Now look for the end. */
618:         if (lstart == lend) {
619:           lend = lstart+1;
620:           while (lend < llen && colors[lend] == colors[lstart]) ++lend;
621:         }
622:         /* Now check whether the identified color segment matches l. */
623:         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);
624:       }
625:       color = (PetscMPIInt)(colors[lstart] == l);
626:       /* Check whether a proper subcommunicator exists. */
627:       MPIU_Allreduce(&color,&subsize,1,MPI_INT,MPI_SUM,comm);

629:       if (subsize == 1) subcomm = PETSC_COMM_SELF;
630:       else if (subsize == size) subcomm = comm;
631:       else {
632:         /* a proper communicator is necessary, so we create it. */
633:         MPI_Comm_split(comm, color, rank, &subcomm);
634:       }
635:       if (colors[lstart] == l) {
636:         /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
637:         ISCreateGeneral(subcomm, lend-lstart,inds+lstart,PETSC_COPY_VALUES,*islist+lcount);
638:         /* Position lstart at the beginning of the next local color. */
639:         lstart = lend;
640:         /* Increment the counter of the local colors split off into an IS. */
641:         ++lcount;
642:       }
643:       if (subsize > 0 && subsize < size) {
644:         /*
645:          Irrespective of color, destroy the split off subcomm:
646:          a subcomm used in the IS creation above is duplicated
647:          into a proper PETSc comm.
648:          */
649:         MPI_Comm_free(&subcomm);
650:       }
651:     } /* for (l = low; l < high; ++l) */
652:   } /* if (low <= high) */
653:   PetscFree2(inds,colors);
654:   return(0);
655: }


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

663:   Not collective.

665:   Input arguments:
666: + a    -  IS to embed
667: . b    -  IS to embed into
668: - drop -  flag indicating whether to drop a's indices that are not in b.

670:   Output arguments:
671: . c    -  local embedding indices

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

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

680:   Level: advanced

682: .seealso ISLocalToGlobalMapping
683:  @*/
684: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
685: {
686:   PetscErrorCode             ierr;
687:   ISLocalToGlobalMapping     ltog;
688:   ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
689:   PetscInt                   alen, clen, *cindices, *cindices2;
690:   const PetscInt             *aindices;

696:   ISLocalToGlobalMappingCreateIS(b, &ltog);
697:   ISGetLocalSize(a, &alen);
698:   ISGetIndices(a, &aindices);
699:   PetscMalloc1(alen, &cindices);
700:   if (!drop) gtoltype = IS_GTOLM_MASK;
701:   ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
702:   ISLocalToGlobalMappingDestroy(&ltog);
703:   if (clen != alen) {
704:     cindices2 = cindices;
705:     PetscMalloc1(clen, &cindices);
706:     PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
707:     PetscFree(cindices2);
708:   }
709:   ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
710:   return(0);
711: }


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

717:   Not collective.

719:   Input arguments:
720: + f      -  IS to sort
721: - always -  build the permutation even when f's indices are nondecreasin.

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


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

732:   Level: advanced

734: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
735:  @*/
736: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
737: {
738:   PetscErrorCode  ierr;
739:   const PetscInt  *findices;
740:   PetscInt        fsize,*hindices,i;
741:   PetscBool       isincreasing;

746:   ISGetLocalSize(f,&fsize);
747:   ISGetIndices(f,&findices);
748:   *h = NULL;
749:   if (!always) {
750:     isincreasing = PETSC_TRUE;
751:     for (i = 1; i < fsize; ++i) {
752:       if (findices[i] <= findices[i-1]) {
753:         isincreasing = PETSC_FALSE;
754:         break;
755:       }
756:     }
757:     if (isincreasing) {
758:       ISRestoreIndices(f,&findices);
759:       return(0);
760:     }
761:   }
762:   PetscMalloc1(fsize,&hindices);
763:   for (i = 0; i < fsize; ++i) hindices[i] = i;
764:   PetscSortIntWithPermutation(fsize,findices,hindices);
765:   ISRestoreIndices(f,&findices);
766:   ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
767:   return(0);
768: }