Actual source code: sorti.c

petsc-master 2019-10-20
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>
  7:  #include <petsc/private/hashseti.h>

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

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

 20: /* Swap one, two or three pairs. Each pair can have its own type */
 21: #define SWAP1(a,b,t1)               do {t1=a;a=b;b=t1;} while(0)
 22: #define SWAP2(a,b,c,d,t1,t2)        do {t1=a;a=b;b=t1; t2=c;c=d;d=t2;} while(0)
 23: #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)

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

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

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

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

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

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

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

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

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

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

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

163:    Not Collective

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

169:    Level: intermediate

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

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

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

186:    Not Collective

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

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

195:    Level: intermediate

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

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

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

217:    Not Collective

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

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

226:    Level: intermediate

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

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

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

243:    Not Collective

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

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

253:    Level: intermediate

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

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

274: /*@

277:    Not Collective

279:    Input Parameters:
280: +  n  - number of values in the array
281: -  X  - array of integers


284:    Output Parameter:
285: .  dups - True if the array has dups, otherwise false

287:    Level: intermediate

289: .seealso: PetscSortRemoveDupsInt()
290: @*/
292: {
294:   PetscInt       i;
295:   PetscHSetI     ht;
296:   PetscBool      missing;

300:   *dups = PETSC_FALSE;
301:   if (n > 1) {
302:     PetscHSetICreate(&ht);
303:     PetscHSetIResize(ht,n);
304:     for (i=0; i<n; i++) {
305:       PetscHSetIQueryAdd(ht,X[i],&missing);
306:       if(!missing) {*dups = PETSC_TRUE; break;}
307:     }
308:     PetscHSetIDestroy(&ht);
309:   }
310:   return(0);
311: }

313: /*@
314:   PetscFindMPIInt - Finds MPI integer in a sorted array of integers

316:    Not Collective

318:    Input Parameters:
319: +  key - the integer to locate
320: .  n   - number of values in the array
321: -  X   - array of integers

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

326:    Level: intermediate

328: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
329: @*/
330: PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc)
331: {
332:   PetscInt lo = 0,hi = n;

336:   if (!n) {*loc = -1; return(0);}
338:   while (hi - lo > 1) {
339:     PetscInt mid = lo + (hi - lo)/2;
340:     if (key < X[mid]) hi = mid;
341:     else               lo = mid;
342:   }
343:   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
344:   return(0);
345: }

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

351:    Not Collective

353:    Input Parameters:
354: +  n  - number of values
355: .  X  - array of integers
356: -  Y  - second array of integers

358:    Level: intermediate

360: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
361: @*/
362: PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
363: {
365:   PetscInt       pivot,t1,t2;

368:   QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr);
369:   return(0);
370: }

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

376:    Not Collective

378:    Input Parameters:
379: +  n  - number of values
380: .  X  - array of integers
381: .  Y  - second array of integers (first array of the pair)
382: -  Z  - third array of integers  (second array of the pair)

384:    Level: intermediate

386: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
387: @*/
388: PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
389: {
391:   PetscInt       pivot,t1,t2,t3;

394:   QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr);
395:   return(0);
396: }

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

401:    Not Collective

403:    Input Parameters:
404: +  n  - number of values
405: -  X  - array of integers

407:    Level: intermediate

409: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
410: @*/
411: PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
412: {
414:   PetscMPIInt    pivot,t1;

417:   QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
418:   return(0);
419: }

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

424:    Not Collective

426:    Input Parameters:
427: +  n  - number of values
428: -  X  - array of integers

430:    Output Parameter:
431: .  n - number of non-redundant values

433:    Level: intermediate

435: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
436: @*/
437: PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
438: {
440:   PetscInt       i,s = 0,N = *n, b = 0;

443:   PetscSortMPIInt(N,X);
444:   for (i=0; i<N-1; i++) {
445:     if (X[b+s+1] != X[b]) {
446:       X[b+1] = X[b+s+1]; b++;
447:     } else s++;
448:   }
449:   *n = N - s;
450:   return(0);
451: }

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

457:    Not Collective

459:    Input Parameters:
460: +  n  - number of values
461: .  X  - array of integers
462: -  Y  - second array of integers

464:    Level: intermediate

466: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
467: @*/
468: PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
469: {
471:   PetscMPIInt    pivot,t1,t2;

474:   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
475:   return(0);
476: }

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

482:    Not Collective

484:    Input Parameters:
485: +  n  - number of values
486: .  X  - array of MPI integers
487: -  Y  - second array of Petsc integers

489:    Level: intermediate

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

493: .seealso: PetscSortMPIIntWithArray()
494: @*/
495: PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
496: {
498:   PetscMPIInt    pivot,t1;
499:   PetscInt       t2;

502:   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
503:   return(0);
504: }

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

510:    Not Collective

512:    Input Parameters:
513: +  n  - number of values
514: .  X  - array of integers
515: -  Y  - second array of scalars

517:    Level: intermediate

519: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
520: @*/
521: PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
522: {
524:   PetscInt       pivot,t1;
525:   PetscScalar    t2;

528:   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
529:   return(0);
530: }

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

537:    Not Collective

539:    Input Parameters:
540: +  n  - number of values
541: .  X  - array of integers
542: .  Y  - second array of data
543: .  size - sizeof elements in the data array in bytes
544: -  t2   - workspace of "size" bytes used when sorting

546:    Level: intermediate

548: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
549: @*/
550: PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
551: {
553:   char           *YY = (char*)Y;
554:   PetscInt       i,j,p,t1,pivot,hi=n-1,l,r;

557:   if (n<8) {
558:     for (i=0; i<n; i++) {
559:       pivot = X[i];
560:       for (j=i+1; j<n; j++) {
561:         if (pivot > X[j]) {
562:           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
563:           pivot = X[i];
564:         }
565:       }
566:     }
567:   } else {
568:     /* Two way partition */
569:     p     = MEDIAN(X,hi);
570:     pivot = X[p];
571:     l     = 0;
572:     r     = hi;
573:     while(1) {
574:       while (X[l] < pivot) l++;
575:       while (X[r] > pivot) r--;
576:       if (l >= r) {r++; break;}
577:       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
578:       l++;
579:       r--;
580:     }
581:     PetscSortIntWithDataArray(l,X,Y,size,t2);
582:     PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);
583:   }
584:   return(0);
585: }

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

590:    Not Collective

592:    Input Parameters:
593: +  an  - number of values in the first array
594: .  aI  - first sorted array of integers
595: .  bn  - number of values in the second array
596: -  bI  - second array of integers

598:    Output Parameters:
599: +  n   - number of values in the merged array
600: -  L   - merged sorted array, this is allocated if an array is not provided

602:    Level: intermediate

604: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
605: @*/
606: PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
607: {
609:   PetscInt       *L_ = *L, ak, bk, k;

611:   if (!L_) {
612:     PetscMalloc1(an+bn, L);
613:     L_   = *L;
614:   }
615:   k = ak = bk = 0;
616:   while (ak < an && bk < bn) {
617:     if (aI[ak] == bI[bk]) {
618:       L_[k] = aI[ak];
619:       ++ak;
620:       ++bk;
621:       ++k;
622:     } else if (aI[ak] < bI[bk]) {
623:       L_[k] = aI[ak];
624:       ++ak;
625:       ++k;
626:     } else {
627:       L_[k] = bI[bk];
628:       ++bk;
629:       ++k;
630:     }
631:   }
632:   if (ak < an) {
633:     PetscArraycpy(L_+k,aI+ak,an-ak);
634:     k   += (an-ak);
635:   }
636:   if (bk < bn) {
637:     PetscArraycpy(L_+k,bI+bk,bn-bk);
638:     k   += (bn-bk);
639:   }
640:   *n = k;
641:   return(0);
642: }

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

649:    Not Collective

651:    Input Parameters:
652: +  an  - number of values in the first array
653: .  aI  - first sorted array of integers
654: .  aJ  - first additional array of integers
655: .  bn  - number of values in the second array
656: .  bI  - second array of integers
657: -  bJ  - second additional of integers

659:    Output Parameters:
660: +  n   - number of values in the merged array (== an + bn)
661: .  L   - merged sorted array
662: -  J   - merged additional array

664:    Notes:
665:     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
666:    Level: intermediate

668: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
669: @*/
670: PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
671: {
673:   PetscInt       n_, *L_, *J_, ak, bk, k;

678:   n_ = an + bn;
679:   *n = n_;
680:   if (!*L) {
681:     PetscMalloc1(n_, L);
682:   }
683:   L_ = *L;
684:   if (!*J) {
685:     PetscMalloc1(n_, J);
686:   }
687:   J_   = *J;
688:   k = ak = bk = 0;
689:   while (ak < an && bk < bn) {
690:     if (aI[ak] <= bI[bk]) {
691:       L_[k] = aI[ak];
692:       J_[k] = aJ[ak];
693:       ++ak;
694:       ++k;
695:     } else {
696:       L_[k] = bI[bk];
697:       J_[k] = bJ[bk];
698:       ++bk;
699:       ++k;
700:     }
701:   }
702:   if (ak < an) {
703:     PetscArraycpy(L_+k,aI+ak,an-ak);
704:     PetscArraycpy(J_+k,aJ+ak,an-ak);
705:     k   += (an-ak);
706:   }
707:   if (bk < bn) {
708:     PetscArraycpy(L_+k,bI+bk,bn-bk);
709:     PetscArraycpy(J_+k,bJ+bk,bn-bk);
710:   }
711:   return(0);
712: }

714: /*@
715:    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.

717:    Not Collective

719:    Input Parameters:
720: +  an  - number of values in the first array
721: .  aI  - first sorted array of integers
722: .  bn  - number of values in the second array
723: -  bI  - second array of integers

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

729:    Level: intermediate

731: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
732: @*/
733: PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
734: {
736:   PetscInt       ai,bi,k;

739:   if (!*L) {PetscMalloc1((an+bn),L);}
740:   for (ai=0,bi=0,k=0; ai<an || bi<bn; ) {
741:     PetscInt t = -1;
742:     for ( ; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
743:     for ( ; bi<bn && bI[bi] == t; bi++);
744:     for ( ; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
745:     for ( ; ai<an && aI[ai] == t; ai++);
746:   }
747:   *n = k;
748:   return(0);
749: }

751: /*@C
752:    PetscProcessTree - Prepares tree data to be displayed graphically

754:    Not Collective

756:    Input Parameters:
757: +  n  - number of values
758: .  mask - indicates those entries in the tree, location 0 is always masked
759: -  parentid - indicates the parent of each entry

761:    Output Parameters:
762: +  Nlevels - the number of levels
763: .  Level - for each node tells its level
764: .  Levelcnts - the number of nodes on each level
765: .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
766: -  Column - for each id tells its column index

768:    Level: developer

770:    Notes:
771:     This code is not currently used

773: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
774: @*/
775: PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
776: {
777:   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
779:   PetscBool      done = PETSC_FALSE;

782:   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
783:   for (i=0; i<n; i++) {
784:     if (mask[i]) continue;
785:     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
786:     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
787:   }

789:   for (i=0; i<n; i++) {
790:     if (!mask[i]) nmask++;
791:   }

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

796:   level[0] = 1;
797:   while (!done) {
798:     done = PETSC_TRUE;
799:     for (i=0; i<n; i++) {
800:       if (mask[i]) continue;
801:       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
802:       else if (!level[i]) done = PETSC_FALSE;
803:     }
804:   }
805:   for (i=0; i<n; i++) {
806:     level[i]--;
807:     nlevels = PetscMax(nlevels,level[i]);
808:   }

810:   /* count the number of nodes on each level and its max */
811:   PetscCalloc1(nlevels,&levelcnt);
812:   for (i=0; i<n; i++) {
813:     if (mask[i]) continue;
814:     levelcnt[level[i]-1]++;
815:   }
816:   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);

818:   /* for each level sort the ids by the parent id */
819:   PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
820:   PetscMalloc1(nmask,&idbylevel);
821:   for (j=1; j<=nlevels;j++) {
822:     cnt = 0;
823:     for (i=0; i<n; i++) {
824:       if (mask[i]) continue;
825:       if (level[i] != j) continue;
826:       workid[cnt]         = i;
827:       workparentid[cnt++] = parentid[i];
828:     }
829:     /*  PetscIntView(cnt,workparentid,0);
830:     PetscIntView(cnt,workid,0);
831:     PetscSortIntWithArray(cnt,workparentid,workid);
832:     PetscIntView(cnt,workparentid,0);
833:     PetscIntView(cnt,workid,0);*/
834:     PetscArraycpy(idbylevel+tcnt,workid,cnt);
835:     tcnt += cnt;
836:   }
837:   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
838:   PetscFree2(workid,workparentid);

840:   /* for each node list its column */
841:   PetscMalloc1(n,&column);
842:   cnt = 0;
843:   for (j=0; j<nlevels; j++) {
844:     for (i=0; i<levelcnt[j]; i++) {
845:       column[idbylevel[cnt++]] = i;
846:     }
847:   }

849:   *Nlevels   = nlevels;
850:   *Level     = level;
851:   *Levelcnt  = levelcnt;
852:   *Idbylevel = idbylevel;
853:   *Column    = column;
854:   return(0);
855: }