Actual source code: isdiff.c

petsc-master 2019-06-15
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


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

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


 41:   ISGetIndices(is1,&i1);
 42:   ISGetLocalSize(is1,&n1);

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

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

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

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

 86:   PetscBTDestroy(&mask);
 87:   return(0);
 88: }

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

 93:    Only sequential version (at the moment)

 95:    Input Parameter:
 96: +  is1 - index set to be extended
 97: -  is2 - index values to be added

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

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

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

107:    Level: intermediate

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


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

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

129:   ISSorted(is1,&f);
130:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 1 is not sorted");
131:   ISSorted(is2,&f);
132:   if (!f) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Arg 2 is not sorted");

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 Parameter:
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()


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


254:   if (!is1 && !is2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
255:   if (!is1) {ISDuplicate(is2, isout);return(0);}
256:   if (!is2) {ISDuplicate(is1, isout);return(0);}
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 intersection of two index sets, by sorting and comparing.

305:    Collective on IS

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

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

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

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

320:    Level: intermediate

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


325: @*/
326: PetscErrorCode ISIntersect(IS is1,IS is2,IS *isout)
327: {
329:   PetscInt       i,n1,n2,nout,*iout;
330:   const PetscInt *i1,*i2;
331:   IS             is1sorted = NULL, is2sorted = NULL;
332:   PetscBool      sorted;
333:   MPI_Comm       comm;


340:   ISGetLocalSize(is1,&n1);
341:   ISGetLocalSize(is2,&n2);
342:   if (n1 < n2) {
343:     IS tempis = is1;
344:     int ntemp = n1;

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

372:   PetscMalloc1(n2,&iout);

374:   for (i = 0, nout = 0; i < n2; i++) {
375:     PetscInt key = i2[i];
376:     PetscInt loc;

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

387:   /* create the new IS containing the sum */
388:   PetscObjectGetComm((PetscObject)is1,&comm);
389:   ISCreateGeneral(comm,nout,iout,PETSC_OWN_POINTER,isout);

391:   ISRestoreIndices(is2sorted,&i2);
392:   ISDestroy(&is2sorted);
393:   ISRestoreIndices(is1sorted,&i1);
394:   ISDestroy(&is1sorted);
395:   return(0);
396: }

398: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
399: {

403:   *isect = NULL;
404:   if (is2 && is1) {
405:     char           composeStr[33] = {0};
406:     PetscObjectId  is2id;

408:     PetscObjectGetId((PetscObject)is2,&is2id);
409:     PetscSNPrintf(composeStr,32,"ISIntersect_Caching_%x",is2id);
410:     PetscObjectQuery((PetscObject) is1, composeStr, (PetscObject *) isect);
411:     if (*isect == NULL) {
412:       ISIntersect(is1, is2, isect);
413:       PetscObjectCompose((PetscObject) is1, composeStr, (PetscObject) *isect);
414:     } else {
415:       PetscObjectReference((PetscObject) *isect);
416:     }
417:   }
418:   return(0);
419: }

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


425:    Collective.

427:    Input Parameter:
428: +  comm    - communicator of the concatenated IS.
429: .  len     - size of islist array (nonnegative)
430: -  islist  - array of index sets

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

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

438:    Level: intermediate

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


443: @*/
444: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
445: {
447:   PetscInt i,n,N;
448:   const PetscInt *iidx;
449:   PetscInt *idx;

453: #if defined(PETSC_USE_DEBUG)
455: #endif
457:   if (!len) {
458:     ISCreateStride(comm, 0,0,0, isout);
459:     return(0);
460:   }
461:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %D", len);
462:   N = 0;
463:   for (i = 0; i < len; ++i) {
464:     if (islist[i]) {
465:       ISGetLocalSize(islist[i], &n);
466:       N   += n;
467:     }
468:   }
469:   PetscMalloc1(N, &idx);
470:   N = 0;
471:   for (i = 0; i < len; ++i) {
472:     if (islist[i]) {
473:       ISGetLocalSize(islist[i], &n);
474:       ISGetIndices(islist[i], &iidx);
475:       PetscMemcpy(idx+N,iidx, sizeof(PetscInt)*n);
476:       ISRestoreIndices(islist[i], &iidx);
477:       N   += n;
478:     }
479:   }
480:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
481:   return(0);
482: }

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


490:   Collective.

492:   Input arguments:
493: + comm    -  MPI_Comm
494: . listlen -  IS list length
495: - islist  -  IS list

497:   Output arguments:
498: + xis -  domain IS
499: - yis -  range  IS

501:   Level: advanced

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

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

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

516: .seealso ISPairToList()
517: @*/
518: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
519: {
521:   PetscInt       ncolors, *colors,i, leni,len,*xinds, *yinds,k,j;
522:   const PetscInt *indsi;

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


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

559:   Collective on indis.

561:   Input arguments:
562: + xis -  domain IS
563: - yis -  range IS

565:   Output arguments:
566: + listlen -  length of islist
567: - islist  -  list of ISs breaking up indis by color

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

573:   Level: advanced

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

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

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


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

684:   Not collective.

686:   Input arguments:
687: + a    -  IS to embed
688: . b    -  IS to embed into
689: - drop -  flag indicating whether to drop a's indices that are not in b.

691:   Output arguments:
692: . c    -  local embedding indices

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

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

701:   Level: advanced

703: .seealso ISLocalToGlobalMapping
704:  @*/
705: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
706: {
707:   PetscErrorCode             ierr;
708:   ISLocalToGlobalMapping     ltog;
709:   ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
710:   PetscInt                   alen, clen, *cindices, *cindices2;
711:   const PetscInt             *aindices;

717:   ISLocalToGlobalMappingCreateIS(b, &ltog);
718:   ISGetLocalSize(a, &alen);
719:   ISGetIndices(a, &aindices);
720:   PetscMalloc1(alen, &cindices);
721:   if (!drop) gtoltype = IS_GTOLM_MASK;
722:   ISGlobalToLocalMappingApply(ltog,gtoltype,alen,aindices,&clen,cindices);
723:   ISLocalToGlobalMappingDestroy(&ltog);
724:   if (clen != alen) {
725:     cindices2 = cindices;
726:     PetscMalloc1(clen, &cindices);
727:     PetscMemcpy(cindices,cindices2,clen*sizeof(PetscInt));
728:     PetscFree(cindices2);
729:   }
730:   ISCreateGeneral(PETSC_COMM_SELF,clen,cindices,PETSC_OWN_POINTER,c);
731:   return(0);
732: }


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

738:   Not collective.

740:   Input arguments:
741: + f      -  IS to sort
742: - always -  build the permutation even when f's indices are nondecreasin.

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


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

753:   Level: advanced

755: .seealso ISLocalToGlobalMapping, ISSort(), PetscIntSortWithPermutation()
756:  @*/
757: PetscErrorCode ISSortPermutation(IS f,PetscBool always,IS *h)
758: {
759:   PetscErrorCode  ierr;
760:   const PetscInt  *findices;
761:   PetscInt        fsize,*hindices,i;
762:   PetscBool       isincreasing;

767:   ISGetLocalSize(f,&fsize);
768:   ISGetIndices(f,&findices);
769:   *h = NULL;
770:   if (!always) {
771:     isincreasing = PETSC_TRUE;
772:     for (i = 1; i < fsize; ++i) {
773:       if (findices[i] <= findices[i-1]) {
774:         isincreasing = PETSC_FALSE;
775:         break;
776:       }
777:     }
778:     if (isincreasing) {
779:       ISRestoreIndices(f,&findices);
780:       return(0);
781:     }
782:   }
783:   PetscMalloc1(fsize,&hindices);
784:   for (i = 0; i < fsize; ++i) hindices[i] = i;
785:   PetscSortIntWithPermutation(fsize,findices,hindices);
786:   ISRestoreIndices(f,&findices);
787:   ISCreateGeneral(PETSC_COMM_SELF,fsize,hindices,PETSC_OWN_POINTER,h);
788:   return(0);
789: }