Actual source code: sorti.c

petsc-master 2019-09-22
Report Typos and Errors

  2: /*
  3:    This file contains routines for sorting integers. Values are sorted in place.
  4:    One can use src/sys/examples/tests/ex52.c for benchmarking.
  5:  */
  6:  #include <petsc/private/petscimpl.h>

  8: #define MEDIAN3(v,a,b,c)                                                        \
  9:   (v[a]<v[b]                                                                    \
 10:    ? (v[b]<v[c]                                                                 \
 11:       ? (b)                                                                     \
 12:       : (v[a]<v[c] ? (c) : (a)))                                                \
 13:    : (v[c]<v[b]                                                                 \
 14:       ? (b)                                                                     \
 15:       : (v[a]<v[c] ? (a) : (c))))

 17: #define MEDIAN(v,right) MEDIAN3(v,right/4,right/2,right/4*3)

 19: /* Swap one, two or three pairs. Each pair can have its own type */
 20: #define SWAP1(a,b,t1)               do {t1=a;a=b;b=t1;} while(0)
 21: #define SWAP2(a,b,c,d,t1,t2)        do {t1=a;a=b;b=t1; t2=c;c=d;d=t2;} while(0)
 22: #define SWAP3(a,b,c,d,e,f,t1,t2,t3) do {t1=a;a=b;b=t1; t2=c;c=d;d=t2; t3=e;e=f;f=t3;} while(0)

 24: /* Swap a & b, *c & *d. c, d, t2 are pointers to a type of size <siz> */
 25: #define SWAP2Data(a,b,c,d,t1,t2,siz)                                             \
 26:   do {                                                                           \
 28:     t1=a; a=b; b=t1;                                                             \
 29:     PetscMemcpy(t2,c,siz);                                  \
 30:     PetscMemcpy(c,d,siz);                                   \
 31:     PetscMemcpy(d,t2,siz);                                  \
 32:   } while(0)

 34: /*
 35:    Partition X[lo,hi] into two parts: X[lo,l) <= pivot; X[r,hi] > pivot

 37:    Input Parameters:
 38:     + X         - array to partition
 39:     . pivot     - a pivot of X[]
 40:     . t1        - temp variable for X
 41:     - lo,hi     - lower and upper bound of the array

 43:    Output Parameters:
 44:     + l,r       - of type PetscInt

 46:    Notes:
 47:     The TwoWayPartition2/3 variants also partition other arrays along with X.
 48:     These arrays can have different types, so they provide their own temp t2,t3
 49:  */
 50: #define TwoWayPartition1(X,pivot,t1,lo,hi,l,r)                                   \
 51:   do {                                                                           \
 52:     l = lo;                                                                      \
 53:     r = hi;                                                                      \
 54:     while(1) {                                                                   \
 55:       while (X[l] < pivot) l++;                                                  \
 56:       while (X[r] > pivot) r--;                                                  \
 57:       if (l >= r) {r++; break;}                                                  \
 58:       SWAP1(X[l],X[r],t1);                                                       \
 59:       l++;                                                                       \
 60:       r--;                                                                       \
 61:     }                                                                            \
 62:   } while(0)

 64: #define TwoWayPartition2(X,Y,pivot,t1,t2,lo,hi,l,r)                              \
 65:   do {                                                                           \
 66:     l = lo;                                                                      \
 67:     r = hi;                                                                      \
 68:     while(1) {                                                                   \
 69:       while (X[l] < pivot) l++;                                                  \
 70:       while (X[r] > pivot) r--;                                                  \
 71:       if (l >= r) {r++; break;}                                                  \
 72:       SWAP2(X[l],X[r],Y[l],Y[r],t1,t2);                                          \
 73:       l++;                                                                       \
 74:       r--;                                                                       \
 75:     }                                                                            \
 76:   } while(0)

 78: #define TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,lo,hi,l,r)                         \
 79:   do {                                                                           \
 80:     l = lo;                                                                      \
 81:     r = hi;                                                                      \
 82:     while(1) {                                                                   \
 83:       while (X[l] < pivot) l++;                                                  \
 84:       while (X[r] > pivot) r--;                                                  \
 85:       if (l >= r) {r++; break;}                                                  \
 86:       SWAP3(X[l],X[r],Y[l],Y[r],Z[l],Z[r],t1,t2,t3);                             \
 87:       l++;                                                                       \
 88:       r--;                                                                       \
 89:     }                                                                            \
 90:   } while(0)

 92: /* Templates for similar functions used below */
 93: #define QuickSort1(FuncName,X,n,pivot,t1,ierr)                                   \
 94:   do {                                                                           \
 95:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
 96:     if (n < 8) {                                                                 \
 97:       for (i=0; i<n; i++) {                                                      \
 98:         pivot = X[i];                                                            \
 99:         for (j=i+1; j<n; j++) {                                                  \
100:           if (pivot > X[j]) {                                                    \
101:             SWAP1(X[i],X[j],t1);                                                 \
102:             pivot = X[i];                                                        \
103:           }                                                                      \
104:         }                                                                        \
105:       }                                                                          \
106:     } else {                                                                     \
107:       p     = MEDIAN(X,hi);                                                      \
108:       pivot = X[p];                                                              \
109:       TwoWayPartition1(X,pivot,t1,0,hi,l,r);                                     \
110:       FuncName(l,X);                                       \
111:       FuncName(hi-r+1,X+r);                                \
112:     }                                                                            \
113:   } while(0)

115: #define QuickSort2(FuncName,X,Y,n,pivot,t1,t2,ierr)                              \
116:   do {                                                                           \
117:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
118:     if (n < 8) {                                                                 \
119:       for (i=0; i<n; i++) {                                                      \
120:         pivot = X[i];                                                            \
121:         for (j=i+1; j<n; j++) {                                                  \
122:           if (pivot > X[j]) {                                                    \
123:             SWAP2(X[i],X[j],Y[i],Y[j],t1,t2);                                    \
124:             pivot = X[i];                                                        \
125:           }                                                                      \
126:         }                                                                        \
127:       }                                                                          \
128:     } else {                                                                     \
129:       p     = MEDIAN(X,hi);                                                      \
130:       pivot = X[p];                                                              \
131:       TwoWayPartition2(X,Y,pivot,t1,t2,0,hi,l,r);                                \
132:       FuncName(l,X,Y);                                     \
133:       FuncName(hi-r+1,X+r,Y+r);                            \
134:     }                                                                            \
135:   } while(0)

137: #define QuickSort3(FuncName,X,Y,Z,n,pivot,t1,t2,t3,ierr)                         \
138:   do {                                                                           \
139:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
140:     if (n < 8) {                                                                 \
141:       for (i=0; i<n; i++) {                                                      \
142:         pivot = X[i];                                                            \
143:         for (j=i+1; j<n; j++) {                                                  \
144:           if (pivot > X[j]) {                                                    \
145:             SWAP3(X[i],X[j],Y[i],Y[j],Z[i],Z[j],t1,t2,t3);                       \
146:             pivot = X[i];                                                        \
147:           }                                                                      \
148:         }                                                                        \
149:       }                                                                          \
150:     } else {                                                                     \
151:       p     = MEDIAN(X,hi);                                                      \
152:       pivot = X[p];                                                              \
153:       TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,0,hi,l,r);                           \
154:       FuncName(l,X,Y,Z);                                   \
155:       FuncName(hi-r+1,X+r,Y+r,Z+r);                        \
156:     }                                                                            \
157:   } while(0)

159: /*@
160:    PetscSortInt - Sorts an array of integers in place in increasing order.

162:    Not Collective

164:    Input Parameters:
165: +  n  - number of values
166: -  X  - array of integers

168:    Level: intermediate

170: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
171: @*/
172: PetscErrorCode  PetscSortInt(PetscInt n,PetscInt X[])
173: {
175:   PetscInt       pivot,t1;

178:   QuickSort1(PetscSortInt,X,n,pivot,t1,ierr);
179:   return(0);
180: }

182: /*@
183:    PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array

185:    Not Collective

187:    Input Parameters:
188: +  n  - number of values
189: -  X  - sorted array of integers

191:    Output Parameter:
192: .  n - number of non-redundant values

194:    Level: intermediate

196: .seealso: PetscSortInt()
197: @*/
198: PetscErrorCode  PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[])
199: {
200:   PetscInt i,s = 0,N = *n, b = 0;

203:   for (i=0; i<N-1; i++) {
204:     if (PetscUnlikely(X[b+s+1] < X[b])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Input array is not sorted");
205:     if (X[b+s+1] != X[b]) {
206:       X[b+1] = X[b+s+1]; b++;
207:     } else s++;
208:   }
209:   *n = N - s;
210:   return(0);
211: }

213: /*@
214:    PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries

216:    Not Collective

218:    Input Parameters:
219: +  n  - number of values
220: -  X  - array of integers

222:    Output Parameter:
223: .  n - number of non-redundant values

225:    Level: intermediate

227: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt()
228: @*/
229: PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[])
230: {

234:   PetscSortInt(*n,X);
235:   PetscSortedRemoveDupsInt(n,X);
236:   return(0);
237: }

239: /*@
240:   PetscFindInt - Finds integer in a sorted array of integers

242:    Not Collective

244:    Input Parameters:
245: +  key - the integer to locate
246: .  n   - number of values in the array
247: -  X  - array of integers

249:    Output Parameter:
250: .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go

252:    Level: intermediate

254: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
255: @*/
256: PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt X[], PetscInt *loc)
257: {
258:   PetscInt lo = 0,hi = n;

262:   if (!n) {*loc = -1; return(0);}
264:   while (hi - lo > 1) {
265:     PetscInt mid = lo + (hi - lo)/2;
266:     if (key < X[mid]) hi = mid;
267:     else               lo = mid;
268:   }
269:   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
270:   return(0);
271: }

273: /*@
274:   PetscFindMPIInt - Finds MPI integer in a sorted array of integers

276:    Not Collective

278:    Input Parameters:
279: +  key - the integer to locate
280: .  n   - number of values in the array
281: -  X   - array of integers

283:    Output Parameter:
284: .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go

286:    Level: intermediate

288: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
289: @*/
290: PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc)
291: {
292:   PetscInt lo = 0,hi = n;

296:   if (!n) {*loc = -1; return(0);}
298:   while (hi - lo > 1) {
299:     PetscInt mid = lo + (hi - lo)/2;
300:     if (key < X[mid]) hi = mid;
301:     else               lo = mid;
302:   }
303:   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
304:   return(0);
305: }

307: /*@
308:    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
309:        changes a second array to match the sorted first array.

311:    Not Collective

313:    Input Parameters:
314: +  n  - number of values
315: .  X  - array of integers
316: -  Y  - second array of integers

318:    Level: intermediate

320: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
321: @*/
322: PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
323: {
325:   PetscInt       pivot,t1,t2;

328:   QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr);
329:   return(0);
330: }

332: /*@
333:    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
334:        changes a pair of integer arrays to match the sorted first array.

336:    Not Collective

338:    Input Parameters:
339: +  n  - number of values
340: .  X  - array of integers
341: .  Y  - second array of integers (first array of the pair)
342: -  Z  - third array of integers  (second array of the pair)

344:    Level: intermediate

346: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
347: @*/
348: PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
349: {
351:   PetscInt       pivot,t1,t2,t3;

354:   QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr);
355:   return(0);
356: }

358: /*@
359:    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.

361:    Not Collective

363:    Input Parameters:
364: +  n  - number of values
365: -  X  - array of integers

367:    Level: intermediate

369: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
370: @*/
371: PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
372: {
374:   PetscMPIInt    pivot,t1;

377:   QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
378:   return(0);
379: }

381: /*@
382:    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries

384:    Not Collective

386:    Input Parameters:
387: +  n  - number of values
388: -  X  - array of integers

390:    Output Parameter:
391: .  n - number of non-redundant values

393:    Level: intermediate

395: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
396: @*/
397: PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
398: {
400:   PetscInt       i,s = 0,N = *n, b = 0;

403:   PetscSortMPIInt(N,X);
404:   for (i=0; i<N-1; i++) {
405:     if (X[b+s+1] != X[b]) {
406:       X[b+1] = X[b+s+1]; b++;
407:     } else s++;
408:   }
409:   *n = N - s;
410:   return(0);
411: }

413: /*@
414:    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
415:        changes a second array to match the sorted first array.

417:    Not Collective

419:    Input Parameters:
420: +  n  - number of values
421: .  X  - array of integers
422: -  Y  - second array of integers

424:    Level: intermediate

426: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
427: @*/
428: PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
429: {
431:   PetscMPIInt    pivot,t1,t2;

434:   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
435:   return(0);
436: }

438: /*@
439:    PetscSortMPIIntWithIntArray - Sorts an array of MPI integers in place in increasing order;
440:        changes a second array of Petsc intergers to match the sorted first array.

442:    Not Collective

444:    Input Parameters:
445: +  n  - number of values
446: .  X  - array of MPI integers
447: -  Y  - second array of Petsc integers

449:    Level: intermediate

451:    Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays.

453: .seealso: PetscSortMPIIntWithArray()
454: @*/
455: PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
456: {
458:   PetscMPIInt    pivot,t1;
459:   PetscInt       t2;

462:   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
463:   return(0);
464: }

466: /*@
467:    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
468:        changes a second SCALAR array to match the sorted first INTEGER array.

470:    Not Collective

472:    Input Parameters:
473: +  n  - number of values
474: .  X  - array of integers
475: -  Y  - second array of scalars

477:    Level: intermediate

479: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
480: @*/
481: PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
482: {
484:   PetscInt       pivot,t1;
485:   PetscScalar    t2;

488:   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
489:   return(0);
490: }

492: /*@C
493:    PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
494:        changes a second array to match the sorted first INTEGER array.  Unlike other sort routines, the user must
495:        provide workspace (the size of an element in the data array) to use when sorting.

497:    Not Collective

499:    Input Parameters:
500: +  n  - number of values
501: .  X  - array of integers
502: .  Y  - second array of data
503: .  size - sizeof elements in the data array in bytes
504: -  t2   - workspace of "size" bytes used when sorting

506:    Level: intermediate

508: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
509: @*/
510: PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
511: {
513:   char           *YY = (char*)Y;
514:   PetscInt       i,j,p,t1,pivot,hi=n-1,l,r;

517:   if (n<8) {
518:     for (i=0; i<n; i++) {
519:       pivot = X[i];
520:       for (j=i+1; j<n; j++) {
521:         if (pivot > X[j]) {
522:           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
523:           pivot = X[i];
524:         }
525:       }
526:     }
527:   } else {
528:     /* Two way partition */
529:     p     = MEDIAN(X,hi);
530:     pivot = X[p];
531:     l     = 0;
532:     r     = hi;
533:     while(1) {
534:       while (X[l] < pivot) l++;
535:       while (X[r] > pivot) r--;
536:       if (l >= r) {r++; break;}
537:       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
538:       l++;
539:       r--;
540:     }
541:     PetscSortIntWithDataArray(l,X,Y,size,t2);
542:     PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);
543:   }
544:   return(0);
545: }

547: /*@
548:    PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.

550:    Not Collective

552:    Input Parameters:
553: +  an  - number of values in the first array
554: .  aI  - first sorted array of integers
555: .  bn  - number of values in the second array
556: -  bI  - second array of integers

558:    Output Parameters:
559: +  n   - number of values in the merged array
560: -  L   - merged sorted array, this is allocated if an array is not provided

562:    Level: intermediate

564: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
565: @*/
566: PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
567: {
569:   PetscInt       *L_ = *L, ak, bk, k;

571:   if (!L_) {
572:     PetscMalloc1(an+bn, L);
573:     L_   = *L;
574:   }
575:   k = ak = bk = 0;
576:   while (ak < an && bk < bn) {
577:     if (aI[ak] == bI[bk]) {
578:       L_[k] = aI[ak];
579:       ++ak;
580:       ++bk;
581:       ++k;
582:     } else if (aI[ak] < bI[bk]) {
583:       L_[k] = aI[ak];
584:       ++ak;
585:       ++k;
586:     } else {
587:       L_[k] = bI[bk];
588:       ++bk;
589:       ++k;
590:     }
591:   }
592:   if (ak < an) {
593:     PetscArraycpy(L_+k,aI+ak,an-ak);
594:     k   += (an-ak);
595:   }
596:   if (bk < bn) {
597:     PetscArraycpy(L_+k,bI+bk,bn-bk);
598:     k   += (bn-bk);
599:   }
600:   *n = k;
601:   return(0);
602: }

604: /*@
605:    PetscMergeIntArrayPair -     Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
606:                                 The additional arrays are the same length as sorted arrays and are merged
607:                                 in the order determined by the merging of the sorted pair.

609:    Not Collective

611:    Input Parameters:
612: +  an  - number of values in the first array
613: .  aI  - first sorted array of integers
614: .  aJ  - first additional array of integers
615: .  bn  - number of values in the second array
616: .  bI  - second array of integers
617: -  bJ  - second additional of integers

619:    Output Parameters:
620: +  n   - number of values in the merged array (== an + bn)
621: .  L   - merged sorted array
622: -  J   - merged additional array

624:    Notes:
625:     if L or J point to non-null arrays then this routine will assume they are of the approproate size and use them, otherwise this routine will allocate space for them
626:    Level: intermediate

628: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
629: @*/
630: PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
631: {
633:   PetscInt       n_, *L_, *J_, ak, bk, k;

638:   n_ = an + bn;
639:   *n = n_;
640:   if (!*L) {
641:     PetscMalloc1(n_, L);
642:   }
643:   L_ = *L;
644:   if (!*J) {
645:     PetscMalloc1(n_, J);
646:   }
647:   J_   = *J;
648:   k = ak = bk = 0;
649:   while (ak < an && bk < bn) {
650:     if (aI[ak] <= bI[bk]) {
651:       L_[k] = aI[ak];
652:       J_[k] = aJ[ak];
653:       ++ak;
654:       ++k;
655:     } else {
656:       L_[k] = bI[bk];
657:       J_[k] = bJ[bk];
658:       ++bk;
659:       ++k;
660:     }
661:   }
662:   if (ak < an) {
663:     PetscArraycpy(L_+k,aI+ak,an-ak);
664:     PetscArraycpy(J_+k,aJ+ak,an-ak);
665:     k   += (an-ak);
666:   }
667:   if (bk < bn) {
668:     PetscArraycpy(L_+k,bI+bk,bn-bk);
669:     PetscArraycpy(J_+k,bJ+bk,bn-bk);
670:   }
671:   return(0);
672: }

674: /*@
675:    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.

677:    Not Collective

679:    Input Parameters:
680: +  an  - number of values in the first array
681: .  aI  - first sorted array of integers
682: .  bn  - number of values in the second array
683: -  bI  - second array of integers

685:    Output Parameters:
686: +  n   - number of values in the merged array (<= an + bn)
687: -  L   - merged sorted array, allocated if address of NULL pointer is passed

689:    Level: intermediate

691: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
692: @*/
693: PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
694: {
696:   PetscInt       ai,bi,k;

699:   if (!*L) {PetscMalloc1((an+bn),L);}
700:   for (ai=0,bi=0,k=0; ai<an || bi<bn; ) {
701:     PetscInt t = -1;
702:     for ( ; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
703:     for ( ; bi<bn && bI[bi] == t; bi++);
704:     for ( ; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
705:     for ( ; ai<an && aI[ai] == t; ai++);
706:   }
707:   *n = k;
708:   return(0);
709: }

711: /*@C
712:    PetscProcessTree - Prepares tree data to be displayed graphically

714:    Not Collective

716:    Input Parameters:
717: +  n  - number of values
718: .  mask - indicates those entries in the tree, location 0 is always masked
719: -  parentid - indicates the parent of each entry

721:    Output Parameters:
722: +  Nlevels - the number of levels
723: .  Level - for each node tells its level
724: .  Levelcnts - the number of nodes on each level
725: .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
726: -  Column - for each id tells its column index

728:    Level: developer

730:    Notes:
731:     This code is not currently used

733: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
734: @*/
735: PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
736: {
737:   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
739:   PetscBool      done = PETSC_FALSE;

742:   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
743:   for (i=0; i<n; i++) {
744:     if (mask[i]) continue;
745:     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
746:     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
747:   }

749:   for (i=0; i<n; i++) {
750:     if (!mask[i]) nmask++;
751:   }

753:   /* determine the level in the tree of each node */
754:   PetscCalloc1(n,&level);

756:   level[0] = 1;
757:   while (!done) {
758:     done = PETSC_TRUE;
759:     for (i=0; i<n; i++) {
760:       if (mask[i]) continue;
761:       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
762:       else if (!level[i]) done = PETSC_FALSE;
763:     }
764:   }
765:   for (i=0; i<n; i++) {
766:     level[i]--;
767:     nlevels = PetscMax(nlevels,level[i]);
768:   }

770:   /* count the number of nodes on each level and its max */
771:   PetscCalloc1(nlevels,&levelcnt);
772:   for (i=0; i<n; i++) {
773:     if (mask[i]) continue;
774:     levelcnt[level[i]-1]++;
775:   }
776:   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);

778:   /* for each level sort the ids by the parent id */
779:   PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
780:   PetscMalloc1(nmask,&idbylevel);
781:   for (j=1; j<=nlevels;j++) {
782:     cnt = 0;
783:     for (i=0; i<n; i++) {
784:       if (mask[i]) continue;
785:       if (level[i] != j) continue;
786:       workid[cnt]         = i;
787:       workparentid[cnt++] = parentid[i];
788:     }
789:     /*  PetscIntView(cnt,workparentid,0);
790:     PetscIntView(cnt,workid,0);
791:     PetscSortIntWithArray(cnt,workparentid,workid);
792:     PetscIntView(cnt,workparentid,0);
793:     PetscIntView(cnt,workid,0);*/
794:     PetscArraycpy(idbylevel+tcnt,workid,cnt);
795:     tcnt += cnt;
796:   }
797:   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
798:   PetscFree2(workid,workparentid);

800:   /* for each node list its column */
801:   PetscMalloc1(n,&column);
802:   cnt = 0;
803:   for (j=0; j<nlevels; j++) {
804:     for (i=0; i<levelcnt[j]; i++) {
805:       column[idbylevel[cnt++]] = i;
806:     }
807:   }

809:   *Nlevels   = nlevels;
810:   *Level     = level;
811:   *Levelcnt  = levelcnt;
812:   *Idbylevel = idbylevel;
813:   *Column    = column;
814:   return(0);
815: }