Actual source code: isdiff.c


  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 Parameters:
 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: {
 31:   PetscInt       i,n1,n2,imin,imax,nout,*iout;
 32:   const PetscInt *i1,*i2;
 33:   PetscBT        mask;
 34:   MPI_Comm       comm;

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

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

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

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

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

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

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

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

 96:    Only sequential version (at the moment)

 98:    Input Parameters:
 99: +  is1 - index set to be extended
100: -  is2 - index values to be added

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

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

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

110:    Level: intermediate

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

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

125:   PetscObjectGetComm((PetscObject)(is1),&comm);
126:   MPI_Comm_size(comm,&size);

129:   ISSorted(is1,&f);
131:   ISSorted(is2,&f);

134:   ISGetLocalSize(is1,&n1);
135:   ISGetLocalSize(is2,&n2);
136:   if (!n2) {
137:     ISDuplicate(is1,is3);
138:     return 0;
139:   }
140:   ISGetIndices(is1,&i1);
141:   ISGetIndices(is2,&i2);

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

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

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

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

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

219:    Collective on IS

221:    Input Parameters:
222: +  is1 - first index set
223: -  is2 - index values to be added

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

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

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

235:    Level: intermediate

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

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


252:   if (!is1) {ISDuplicate(is2, isout));PetscFunctionReturn(0;}
253:   if (!is2) {ISDuplicate(is1, isout));PetscFunctionReturn(0;}
254:   ISGetIndices(is1,&i1);
255:   ISGetLocalSize(is1,&n1);
256:   ISGetIndices(is2,&i2);
257:   ISGetLocalSize(is2,&n2);

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

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

291:   /* create the new IS containing the sum */
292:   PetscObjectGetComm((PetscObject)is1,&comm);
293:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

295:   PetscBTDestroy(&mask);
296:   return 0;
297: }

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

302:    Collective on IS

304:    Input Parameters:
305: +  is1 - first index set
306: -  is2 - second index set

308:    Output Parameters:
309: .  isout - the sorted intersection of is1 and is2

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

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

317:    Level: intermediate

319: .seealso: ISDestroy(), ISView(), ISDifference(), ISSum(), ISExpand()
320: @*/
321: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
322: {
323:   PetscInt       i,n1,n2,nout,*iout;
324:   const PetscInt *i1,*i2;
325:   IS             is1sorted = NULL, is2sorted = NULL;
326:   PetscBool      sorted, lsorted;
327:   MPI_Comm       comm;

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

335:   ISGetLocalSize(is1,&n1);
336:   ISGetLocalSize(is2,&n2);
337:   if (n1 < n2) {
338:     IS       tempis = is1;
339:     PetscInt ntemp = n1;

341:     is1 = is2;
342:     is2 = tempis;
343:     n1  = n2;
344:     n2  = ntemp;
345:   }
346:   ISSorted(is1,&lsorted);
347:   MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
348:   if (!sorted) {
349:     ISDuplicate(is1,&is1sorted);
350:     ISSort(is1sorted);
351:     ISGetIndices(is1sorted,&i1);
352:   } else {
353:     is1sorted = is1;
354:     PetscObjectReference((PetscObject)is1);
355:     ISGetIndices(is1,&i1);
356:   }
357:   ISSorted(is2,&lsorted);
358:   MPIU_Allreduce(&lsorted,&sorted,1,MPIU_BOOL,MPI_LAND,comm);
359:   if (!sorted) {
360:     ISDuplicate(is2,&is2sorted);
361:     ISSort(is2sorted);
362:     ISGetIndices(is2sorted,&i2);
363:   } else {
364:     is2sorted = is2;
365:     PetscObjectReference((PetscObject)is2);
366:     ISGetIndices(is2,&i2);
367:   }

369:   PetscMalloc1(n2,&iout);

371:   for (i = 0, nout = 0; i < n2; i++) {
372:     PetscInt key = i2[i];
373:     PetscInt loc;

375:     ISLocate(is1sorted,key,&loc);
376:     if (loc >= 0) {
377:       if (!nout || iout[nout-1] < key) {
378:         iout[nout++] = key;
379:       }
380:     }
381:   }
382:   PetscRealloc(nout*sizeof(PetscInt),&iout);

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

387:   ISRestoreIndices(is2sorted,&i2);
388:   ISDestroy(&is2sorted);
389:   ISRestoreIndices(is1sorted,&i1);
390:   ISDestroy(&is1sorted);
391:   return 0;
392: }

394: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
395: {
396:   *isect = NULL;
397:   if (is2 && is1) {
398:     char           composeStr[33] = {0};
399:     PetscObjectId  is2id;

401:     PetscObjectGetId((PetscObject)is2,&is2id);
402:     PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%" PetscInt64_FMT,is2id);
403:     PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
404:     if (*isect == NULL) {
405:       ISIntersect(is1, is2, isect);
406:       PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
407:     } else {
408:       PetscObjectReference((PetscObject) *isect);
409:     }
410:   }
411:   return 0;
412: }

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

417:    Collective.

419:    Input Parameters:
420: +  comm    - communicator of the concatenated IS.
421: .  len     - size of islist array (nonnegative)
422: -  islist  - array of index sets

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

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

430:    Level: intermediate

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

434: @*/
435: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
436: {
437:   PetscInt i,n,N;
438:   const PetscInt *iidx;
439:   PetscInt *idx;

442:   if (PetscDefined(USE_DEBUG)) {
444:   }
446:   if (!len) {
447:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout);
448:     return 0;
449:   }
451:   N = 0;
452:   for (i = 0; i < len; ++i) {
453:     if (islist[i]) {
454:       ISGetLocalSize(islist[i], &n);
455:       N   += n;
456:     }
457:   }
458:   PetscMalloc1(N, &idx);
459:   N = 0;
460:   for (i = 0; i < len; ++i) {
461:     if (islist[i]) {
462:       ISGetLocalSize(islist[i], &n);
463:       ISGetIndices(islist[i], &iidx);
464:       PetscArraycpy(idx+N, iidx, n);
465:       ISRestoreIndices(islist[i], &iidx);
466:       N   += n;
467:     }
468:   }
469:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
470:   return 0;
471: }

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

478:   Collective.

480:   Input arguments:
481: + comm    -  MPI_Comm
482: . listlen -  IS list length
483: - islist  -  IS list

485:   Output arguments:
486: + xis -  domain IS
487: - yis -  range  IS

489:   Level: advanced

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

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

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

504: .seealso ISPairToList()
505: @*/
506: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
507: {
508:   PetscInt       ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
509:   const PetscInt *indsi;

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

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

544:   Collective on indis.

546:   Input arguments:
547: + xis -  domain IS
548: - yis -  range IS

550:   Output arguments:
551: + listlen -  length of islist
552: - islist  -  list of ISs breaking up indis by color

554:   Note:
555:     xis and yis must be of the same length and have congruent communicators.

557:     The resulting ISs have subcommunicators in a "deadlock-free" order (see ISListToPair()).

559:   Level: advanced

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

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

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

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

667:   Not collective.

669:   Input arguments:
670: + a    -  IS to embed
671: . b    -  IS to embed into
672: - drop -  flag indicating whether to drop a's indices that are not in b.

674:   Output arguments:
675: . c    -  local embedding indices

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

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

684:   Level: advanced

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

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

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

718:   Not collective.

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

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

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()
735:  @*/
736: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
737: {
738:   const PetscInt  *findices;
739:   PetscInt        fsize,*hindices,i;
740:   PetscBool       isincreasing;

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