Actual source code: sorti.c

petsc-master 2019-12-13
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: /*
 66:    Partition X[lo,hi] into two parts: X[lo,l) >= pivot; X[r,hi] < pivot

 68:    Input Parameters:
 69:     + X         - array to partition
 70:     . pivot     - a pivot of X[]
 71:     . t1        - temp variable for X
 72:     - lo,hi     - lower and upper bound of the array

 74:    Output Parameters:
 75:     + l,r       - of type PetscInt

 77:    Notes:
 78:     The TwoWayPartition2/3 variants also partition other arrays along with X.
 79:     These arrays can have different types, so they provide their own temp t2,t3
 80:  */
 81: #define TwoWayPartitionReverse1(X,pivot,t1,lo,hi,l,r)                                   \
 82:   do {                                                                           \
 83:     l = lo;                                                                      \
 84:     r = hi;                                                                      \
 85:     while(1) {                                                                   \
 86:       while (X[l] > pivot) l++;                                                  \
 87:       while (X[r] < pivot) r--;                                                  \
 88:       if (l >= r) {r++; break;}                                                  \
 89:       SWAP1(X[l],X[r],t1);                                                       \
 90:       l++;                                                                       \
 91:       r--;                                                                       \
 92:     }                                                                            \
 93:   } while(0)

 95: #define TwoWayPartition2(X,Y,pivot,t1,t2,lo,hi,l,r)                              \
 96:   do {                                                                           \
 97:     l = lo;                                                                      \
 98:     r = hi;                                                                      \
 99:     while(1) {                                                                   \
100:       while (X[l] < pivot) l++;                                                  \
101:       while (X[r] > pivot) r--;                                                  \
102:       if (l >= r) {r++; break;}                                                  \
103:       SWAP2(X[l],X[r],Y[l],Y[r],t1,t2);                                          \
104:       l++;                                                                       \
105:       r--;                                                                       \
106:     }                                                                            \
107:   } while(0)

109: #define TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,lo,hi,l,r)                         \
110:   do {                                                                           \
111:     l = lo;                                                                      \
112:     r = hi;                                                                      \
113:     while(1) {                                                                   \
114:       while (X[l] < pivot) l++;                                                  \
115:       while (X[r] > pivot) r--;                                                  \
116:       if (l >= r) {r++; break;}                                                  \
117:       SWAP3(X[l],X[r],Y[l],Y[r],Z[l],Z[r],t1,t2,t3);                             \
118:       l++;                                                                       \
119:       r--;                                                                       \
120:     }                                                                            \
121:   } while(0)

123: /* Templates for similar functions used below */
124: #define QuickSort1(FuncName,X,n,pivot,t1,ierr)                                   \
125:   do {                                                                           \
126:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
127:     if (n < 8) {                                                                 \
128:       for (i=0; i<n; i++) {                                                      \
129:         pivot = X[i];                                                            \
130:         for (j=i+1; j<n; j++) {                                                  \
131:           if (pivot > X[j]) {                                                    \
132:             SWAP1(X[i],X[j],t1);                                                 \
133:             pivot = X[i];                                                        \
134:           }                                                                      \
135:         }                                                                        \
136:       }                                                                          \
137:     } else {                                                                     \
138:       p     = MEDIAN(X,hi);                                                      \
139:       pivot = X[p];                                                              \
140:       TwoWayPartition1(X,pivot,t1,0,hi,l,r);                                     \
141:       FuncName(l,X);                                       \
142:       FuncName(hi-r+1,X+r);                                \
143:     }                                                                            \
144:   } while(0)

146: /* Templates for similar functions used below */
147: #define QuickSortReverse1(FuncName,X,n,pivot,t1,ierr)                            \
148:   do {                                                                           \
149:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
150:     if (n < 8) {                                                                 \
151:       for (i=0; i<n; i++) {                                                      \
152:         pivot = X[i];                                                            \
153:         for (j=i+1; j<n; j++) {                                                  \
154:           if (pivot < X[j]) {                                                    \
155:             SWAP1(X[i],X[j],t1);                                                 \
156:             pivot = X[i];                                                        \
157:           }                                                                      \
158:         }                                                                        \
159:       }                                                                          \
160:     } else {                                                                     \
161:       p     = MEDIAN(X,hi);                                                      \
162:       pivot = X[p];                                                              \
163:       TwoWayPartitionReverse1(X,pivot,t1,0,hi,l,r);                              \
164:       FuncName(l,X);                                       \
165:       FuncName(hi-r+1,X+r);                                \
166:     }                                                                            \
167:   } while(0)

169: #define QuickSort2(FuncName,X,Y,n,pivot,t1,t2,ierr)                              \
170:   do {                                                                           \
171:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
172:     if (n < 8) {                                                                 \
173:       for (i=0; i<n; i++) {                                                      \
174:         pivot = X[i];                                                            \
175:         for (j=i+1; j<n; j++) {                                                  \
176:           if (pivot > X[j]) {                                                    \
177:             SWAP2(X[i],X[j],Y[i],Y[j],t1,t2);                                    \
178:             pivot = X[i];                                                        \
179:           }                                                                      \
180:         }                                                                        \
181:       }                                                                          \
182:     } else {                                                                     \
183:       p     = MEDIAN(X,hi);                                                      \
184:       pivot = X[p];                                                              \
185:       TwoWayPartition2(X,Y,pivot,t1,t2,0,hi,l,r);                                \
186:       FuncName(l,X,Y);                                     \
187:       FuncName(hi-r+1,X+r,Y+r);                            \
188:     }                                                                            \
189:   } while(0)

191: #define QuickSort3(FuncName,X,Y,Z,n,pivot,t1,t2,t3,ierr)                         \
192:   do {                                                                           \
193:     PetscInt i,j,p,l,r,hi=n-1;                                                   \
194:     if (n < 8) {                                                                 \
195:       for (i=0; i<n; i++) {                                                      \
196:         pivot = X[i];                                                            \
197:         for (j=i+1; j<n; j++) {                                                  \
198:           if (pivot > X[j]) {                                                    \
199:             SWAP3(X[i],X[j],Y[i],Y[j],Z[i],Z[j],t1,t2,t3);                       \
200:             pivot = X[i];                                                        \
201:           }                                                                      \
202:         }                                                                        \
203:       }                                                                          \
204:     } else {                                                                     \
205:       p     = MEDIAN(X,hi);                                                      \
206:       pivot = X[p];                                                              \
207:       TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,0,hi,l,r);                           \
208:       FuncName(l,X,Y,Z);                                   \
209:       FuncName(hi-r+1,X+r,Y+r,Z+r);                        \
210:     }                                                                            \
211:   } while(0)

213: /*@
214:    PetscSortedInt - Determines whether the array is sorted.

216:    Not Collective

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

222:    Output Parameters:
223: .  sorted - flag whether the array is sorted

225:    Level: intermediate

227: .seealso: PetscSortInt(), PetscSortedMPIInt(), PetscSortedReal()
228: @*/
229: PetscErrorCode  PetscSortedInt(PetscInt n,const PetscInt X[],PetscBool *sorted)
230: {
232:   PetscSorted(n,X,*sorted);
233:   return(0);
234: }

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

239:    Not Collective

241:    Input Parameters:
242: +  n  - number of values
243: -  X  - array of integers

245:    Level: intermediate

247: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
248: @*/
249: PetscErrorCode  PetscSortInt(PetscInt n,PetscInt X[])
250: {
252:   PetscInt       pivot,t1;

255:   QuickSort1(PetscSortInt,X,n,pivot,t1,ierr);
256:   return(0);
257: }

259: /*@
260:    PetscSortReverseInt - Sorts an array of integers in place in decreasing order.

262:    Not Collective

264:    Input Parameters:
265: +  n  - number of values
266: -  X  - array of integers

268:    Level: intermediate

270: .seealso: PetscSortInt(), PetscSortIntWithPermutation()
271: @*/
272: PetscErrorCode  PetscSortReverseInt(PetscInt n,PetscInt X[])
273: {
275:   PetscInt       pivot,t1;

278:   QuickSortReverse1(PetscSortReverseInt,X,n,pivot,t1,ierr);
279:   return(0);
280: }

282: /*@
283:    PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array

285:    Not Collective

287:    Input Parameters:
288: +  n  - number of values
289: -  X  - sorted array of integers

291:    Output Parameter:
292: .  n - number of non-redundant values

294:    Level: intermediate

296: .seealso: PetscSortInt()
297: @*/
298: PetscErrorCode  PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[])
299: {
300:   PetscInt i,s = 0,N = *n, b = 0;

304:   for (i=0; i<N-1; i++) {
305:     if (X[b+s+1] != X[b]) {
306:       X[b+1] = X[b+s+1]; b++;
307:     } else s++;
308:   }
309:   *n = N - s;
310:   return(0);
311: }

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

316:    Not Collective

318:    Input Parameters:
319: +  n  - number of values
320: -  X  - array of integers

322:    Output Parameter:
323: .  n - number of non-redundant values

325:    Level: intermediate

327: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt()
328: @*/
329: PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[])
330: {

334:   PetscSortInt(*n,X);
335:   PetscSortedRemoveDupsInt(n,X);
336:   return(0);
337: }

339: /*@
340:   PetscFindInt - Finds integer in a sorted array of integers

342:    Not Collective

344:    Input Parameters:
345: +  key - the integer to locate
346: .  n   - number of values in the array
347: -  X  - array of integers

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

352:    Level: intermediate

354: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
355: @*/
356: PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt X[], PetscInt *loc)
357: {
358:   PetscInt lo = 0,hi = n;

362:   if (!n) {*loc = -1; return(0);}
365:   while (hi - lo > 1) {
366:     PetscInt mid = lo + (hi - lo)/2;
367:     if (key < X[mid]) hi = mid;
368:     else               lo = mid;
369:   }
370:   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
371:   return(0);
372: }

374: /*@

377:    Not Collective

379:    Input Parameters:
380: +  n  - number of values in the array
381: -  X  - array of integers


384:    Output Parameter:
385: .  dups - True if the array has dups, otherwise false

387:    Level: intermediate

389: .seealso: PetscSortRemoveDupsInt()
390: @*/
392: {
394:   PetscInt       i;
395:   PetscHSetI     ht;
396:   PetscBool      missing;

400:   *dups = PETSC_FALSE;
401:   if (n > 1) {
402:     PetscHSetICreate(&ht);
403:     PetscHSetIResize(ht,n);
404:     for (i=0; i<n; i++) {
405:       PetscHSetIQueryAdd(ht,X[i],&missing);
406:       if(!missing) {*dups = PETSC_TRUE; break;}
407:     }
408:     PetscHSetIDestroy(&ht);
409:   }
410:   return(0);
411: }

413: /*@
414:   PetscFindMPIInt - Finds MPI integer in a sorted array of integers

416:    Not Collective

418:    Input Parameters:
419: +  key - the integer to locate
420: .  n   - number of values in the array
421: -  X   - array of integers

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

426:    Level: intermediate

428: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
429: @*/
430: PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc)
431: {
432:   PetscInt lo = 0,hi = n;

436:   if (!n) {*loc = -1; return(0);}
439:   while (hi - lo > 1) {
440:     PetscInt mid = lo + (hi - lo)/2;
441:     if (key < X[mid]) hi = mid;
442:     else               lo = mid;
443:   }
444:   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
445:   return(0);
446: }

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

452:    Not Collective

454:    Input Parameters:
455: +  n  - number of values
456: .  X  - array of integers
457: -  Y  - second array of integers

459:    Level: intermediate

461: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
462: @*/
463: PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
464: {
466:   PetscInt       pivot,t1,t2;

469:   QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr);
470:   return(0);
471: }

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

477:    Not Collective

479:    Input Parameters:
480: +  n  - number of values
481: .  X  - array of integers
482: .  Y  - second array of integers (first array of the pair)
483: -  Z  - third array of integers  (second array of the pair)

485:    Level: intermediate

487: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
488: @*/
489: PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
490: {
492:   PetscInt       pivot,t1,t2,t3;

495:   QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr);
496:   return(0);
497: }

499: /*@
500:    PetscSortedMPIInt - Determines whether the array is sorted.

502:    Not Collective

504:    Input Parameters:
505: +  n  - number of values
506: -  X  - array of integers

508:    Output Parameters:
509: .  sorted - flag whether the array is sorted

511:    Level: intermediate

513: .seealso: PetscSortMPIInt(), PetscSortedInt(), PetscSortedReal()
514: @*/
515: PetscErrorCode  PetscSortedMPIInt(PetscInt n,const PetscMPIInt X[],PetscBool *sorted)
516: {
518:   PetscSorted(n,X,*sorted);
519:   return(0);
520: }

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

525:    Not Collective

527:    Input Parameters:
528: +  n  - number of values
529: -  X  - array of integers

531:    Level: intermediate

533: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
534: @*/
535: PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
536: {
538:   PetscMPIInt    pivot,t1;

541:   QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
542:   return(0);
543: }

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

548:    Not Collective

550:    Input Parameters:
551: +  n  - number of values
552: -  X  - array of integers

554:    Output Parameter:
555: .  n - number of non-redundant values

557:    Level: intermediate

559: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
560: @*/
561: PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
562: {
564:   PetscInt       i,s = 0,N = *n, b = 0;

567:   PetscSortMPIInt(N,X);
568:   for (i=0; i<N-1; i++) {
569:     if (X[b+s+1] != X[b]) {
570:       X[b+1] = X[b+s+1]; b++;
571:     } else s++;
572:   }
573:   *n = N - s;
574:   return(0);
575: }

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

581:    Not Collective

583:    Input Parameters:
584: +  n  - number of values
585: .  X  - array of integers
586: -  Y  - second array of integers

588:    Level: intermediate

590: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
591: @*/
592: PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
593: {
595:   PetscMPIInt    pivot,t1,t2;

598:   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
599:   return(0);
600: }

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

606:    Not Collective

608:    Input Parameters:
609: +  n  - number of values
610: .  X  - array of MPI integers
611: -  Y  - second array of Petsc integers

613:    Level: intermediate

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

617: .seealso: PetscSortMPIIntWithArray()
618: @*/
619: PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
620: {
622:   PetscMPIInt    pivot,t1;
623:   PetscInt       t2;

626:   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
627:   return(0);
628: }

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

634:    Not Collective

636:    Input Parameters:
637: +  n  - number of values
638: .  X  - array of integers
639: -  Y  - second array of scalars

641:    Level: intermediate

643: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
644: @*/
645: PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
646: {
648:   PetscInt       pivot,t1;
649:   PetscScalar    t2;

652:   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
653:   return(0);
654: }

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

661:    Not Collective

663:    Input Parameters:
664: +  n  - number of values
665: .  X  - array of integers
666: .  Y  - second array of data
667: .  size - sizeof elements in the data array in bytes
668: -  t2   - workspace of "size" bytes used when sorting

670:    Level: intermediate

672: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
673: @*/
674: PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
675: {
677:   char           *YY = (char*)Y;
678:   PetscInt       i,j,p,t1,pivot,hi=n-1,l,r;

681:   if (n<8) {
682:     for (i=0; i<n; i++) {
683:       pivot = X[i];
684:       for (j=i+1; j<n; j++) {
685:         if (pivot > X[j]) {
686:           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
687:           pivot = X[i];
688:         }
689:       }
690:     }
691:   } else {
692:     /* Two way partition */
693:     p     = MEDIAN(X,hi);
694:     pivot = X[p];
695:     l     = 0;
696:     r     = hi;
697:     while(1) {
698:       while (X[l] < pivot) l++;
699:       while (X[r] > pivot) r--;
700:       if (l >= r) {r++; break;}
701:       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
702:       l++;
703:       r--;
704:     }
705:     PetscSortIntWithDataArray(l,X,Y,size,t2);
706:     PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);
707:   }
708:   return(0);
709: }

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

714:    Not Collective

716:    Input Parameters:
717: +  an  - number of values in the first array
718: .  aI  - first sorted array of integers
719: .  bn  - number of values in the second array
720: -  bI  - second array of integers

722:    Output Parameters:
723: +  n   - number of values in the merged array
724: -  L   - merged sorted array, this is allocated if an array is not provided

726:    Level: intermediate

728: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
729: @*/
730: PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
731: {
733:   PetscInt       *L_ = *L, ak, bk, k;

735:   if (!L_) {
736:     PetscMalloc1(an+bn, L);
737:     L_   = *L;
738:   }
739:   k = ak = bk = 0;
740:   while (ak < an && bk < bn) {
741:     if (aI[ak] == bI[bk]) {
742:       L_[k] = aI[ak];
743:       ++ak;
744:       ++bk;
745:       ++k;
746:     } else if (aI[ak] < bI[bk]) {
747:       L_[k] = aI[ak];
748:       ++ak;
749:       ++k;
750:     } else {
751:       L_[k] = bI[bk];
752:       ++bk;
753:       ++k;
754:     }
755:   }
756:   if (ak < an) {
757:     PetscArraycpy(L_+k,aI+ak,an-ak);
758:     k   += (an-ak);
759:   }
760:   if (bk < bn) {
761:     PetscArraycpy(L_+k,bI+bk,bn-bk);
762:     k   += (bn-bk);
763:   }
764:   *n = k;
765:   return(0);
766: }

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

773:    Not Collective

775:    Input Parameters:
776: +  an  - number of values in the first array
777: .  aI  - first sorted array of integers
778: .  aJ  - first additional array of integers
779: .  bn  - number of values in the second array
780: .  bI  - second array of integers
781: -  bJ  - second additional of integers

783:    Output Parameters:
784: +  n   - number of values in the merged array (== an + bn)
785: .  L   - merged sorted array
786: -  J   - merged additional array

788:    Notes:
789:     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
790:    Level: intermediate

792: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
793: @*/
794: PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
795: {
797:   PetscInt       n_, *L_, *J_, ak, bk, k;

802:   n_ = an + bn;
803:   *n = n_;
804:   if (!*L) {
805:     PetscMalloc1(n_, L);
806:   }
807:   L_ = *L;
808:   if (!*J) {
809:     PetscMalloc1(n_, J);
810:   }
811:   J_   = *J;
812:   k = ak = bk = 0;
813:   while (ak < an && bk < bn) {
814:     if (aI[ak] <= bI[bk]) {
815:       L_[k] = aI[ak];
816:       J_[k] = aJ[ak];
817:       ++ak;
818:       ++k;
819:     } else {
820:       L_[k] = bI[bk];
821:       J_[k] = bJ[bk];
822:       ++bk;
823:       ++k;
824:     }
825:   }
826:   if (ak < an) {
827:     PetscArraycpy(L_+k,aI+ak,an-ak);
828:     PetscArraycpy(J_+k,aJ+ak,an-ak);
829:     k   += (an-ak);
830:   }
831:   if (bk < bn) {
832:     PetscArraycpy(L_+k,bI+bk,bn-bk);
833:     PetscArraycpy(J_+k,bJ+bk,bn-bk);
834:   }
835:   return(0);
836: }

838: /*@
839:    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.

841:    Not Collective

843:    Input Parameters:
844: +  an  - number of values in the first array
845: .  aI  - first sorted array of integers
846: .  bn  - number of values in the second array
847: -  bI  - second array of integers

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

853:    Level: intermediate

855: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
856: @*/
857: PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
858: {
860:   PetscInt       ai,bi,k;

863:   if (!*L) {PetscMalloc1((an+bn),L);}
864:   for (ai=0,bi=0,k=0; ai<an || bi<bn; ) {
865:     PetscInt t = -1;
866:     for ( ; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
867:     for ( ; bi<bn && bI[bi] == t; bi++);
868:     for ( ; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
869:     for ( ; ai<an && aI[ai] == t; ai++);
870:   }
871:   *n = k;
872:   return(0);
873: }

875: /*@C
876:    PetscProcessTree - Prepares tree data to be displayed graphically

878:    Not Collective

880:    Input Parameters:
881: +  n  - number of values
882: .  mask - indicates those entries in the tree, location 0 is always masked
883: -  parentid - indicates the parent of each entry

885:    Output Parameters:
886: +  Nlevels - the number of levels
887: .  Level - for each node tells its level
888: .  Levelcnts - the number of nodes on each level
889: .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
890: -  Column - for each id tells its column index

892:    Level: developer

894:    Notes:
895:     This code is not currently used

897: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
898: @*/
899: PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
900: {
901:   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
903:   PetscBool      done = PETSC_FALSE;

906:   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
907:   for (i=0; i<n; i++) {
908:     if (mask[i]) continue;
909:     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
910:     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
911:   }

913:   for (i=0; i<n; i++) {
914:     if (!mask[i]) nmask++;
915:   }

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

920:   level[0] = 1;
921:   while (!done) {
922:     done = PETSC_TRUE;
923:     for (i=0; i<n; i++) {
924:       if (mask[i]) continue;
925:       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
926:       else if (!level[i]) done = PETSC_FALSE;
927:     }
928:   }
929:   for (i=0; i<n; i++) {
930:     level[i]--;
931:     nlevels = PetscMax(nlevels,level[i]);
932:   }

934:   /* count the number of nodes on each level and its max */
935:   PetscCalloc1(nlevels,&levelcnt);
936:   for (i=0; i<n; i++) {
937:     if (mask[i]) continue;
938:     levelcnt[level[i]-1]++;
939:   }
940:   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);

942:   /* for each level sort the ids by the parent id */
943:   PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
944:   PetscMalloc1(nmask,&idbylevel);
945:   for (j=1; j<=nlevels;j++) {
946:     cnt = 0;
947:     for (i=0; i<n; i++) {
948:       if (mask[i]) continue;
949:       if (level[i] != j) continue;
950:       workid[cnt]         = i;
951:       workparentid[cnt++] = parentid[i];
952:     }
953:     /*  PetscIntView(cnt,workparentid,0);
954:     PetscIntView(cnt,workid,0);
955:     PetscSortIntWithArray(cnt,workparentid,workid);
956:     PetscIntView(cnt,workparentid,0);
957:     PetscIntView(cnt,workid,0);*/
958:     PetscArraycpy(idbylevel+tcnt,workid,cnt);
959:     tcnt += cnt;
960:   }
961:   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
962:   PetscFree2(workid,workparentid);

964:   /* for each node list its column */
965:   PetscMalloc1(n,&column);
966:   cnt = 0;
967:   for (j=0; j<nlevels; j++) {
968:     for (i=0; i<levelcnt[j]; i++) {
969:       column[idbylevel[cnt++]] = i;
970:     }
971:   }

973:   *Nlevels   = nlevels;
974:   *Level     = level;
975:   *Levelcnt  = levelcnt;
976:   *Idbylevel = idbylevel;
977:   *Column    = column;
978:   return(0);
979: }

981: /*@
982:   PetscParallelSortedInt - Check whether an integer array, distributed over a communicator, is globally sorted.

984:   Collective

986:   Input Parameters:
987: + comm - the MPI communicator
988: . n - the local number of integers
989: - keys - the local array of integers

991:   Output Parameters:
992: . is_sorted - whether the array is globally sorted

994:   Level: developer

996: .seealso: PetscParallelSortInt()
997: @*/
998: PetscErrorCode PetscParallelSortedInt(MPI_Comm comm, PetscInt n, const PetscInt keys[], PetscBool *is_sorted)
999: {
1000:   PetscBool      sorted;
1001:   PetscInt       i, min, max, prevmax;
1002:   PetscMPIInt    rank;

1006:   sorted = PETSC_TRUE;
1007:   min    = PETSC_MAX_INT;
1008:   max    = PETSC_MIN_INT;
1009:   if (n) {
1010:     min = keys[0];
1011:     max = keys[0];
1012:   }
1013:   for (i = 1; i < n; i++) {
1014:     if (keys[i] < keys[i - 1]) break;
1015:     min = PetscMin(min,keys[i]);
1016:     max = PetscMax(max,keys[i]);
1017:   }
1018:   if (i < n) sorted = PETSC_FALSE;
1019:   prevmax = PETSC_MIN_INT;
1020:   MPI_Exscan(&max, &prevmax, 1, MPIU_INT, MPI_MAX, comm);
1021:   MPI_Comm_rank(comm, &rank);
1022:   if (!rank) prevmax = PETSC_MIN_INT;
1023:   if (prevmax > min) sorted = PETSC_FALSE;
1024:   MPI_Allreduce(&sorted, is_sorted, 1, MPIU_BOOL, MPI_LAND, comm);
1025:   return(0);
1026: }