xref: /petsc/src/sys/utils/sorti.c (revision 3307d110e72ee4e6d2468971620073eb5ff93529)
1 
2 /*
3    This file contains routines for sorting integers. Values are sorted in place.
4    One can use src/sys/tests/ex52.c for benchmarking.
5  */
6 #include <petsc/private/petscimpl.h>                /*I  "petscsys.h"  I*/
7 #include <petsc/private/hashseti.h>
8 
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))))
17 
18 #define MEDIAN(v,right) MEDIAN3(v,right/4,right/2,right/4*3)
19 
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)
24 
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 {                                                                           \
28     t1=a; a=b; b=t1;                                                             \
29     PetscCall(PetscMemcpy(t2,c,siz));                                              \
30     PetscCall(PetscMemcpy(c,d,siz));                                               \
31     PetscCall(PetscMemcpy(d,t2,siz));                                              \
32   } while (0)
33 
34 /*
35    Partition X[lo,hi] into two parts: X[lo,l) <= pivot; X[r,hi] > pivot
36 
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
42 
43    Output Parameters:
44     + l,r       - of type PetscInt
45 
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)
63 
64 /*
65    Partition X[lo,hi] into two parts: X[lo,l) >= pivot; X[r,hi] < pivot
66 
67    Input Parameters:
68     + X         - array to partition
69     . pivot     - a pivot of X[]
70     . t1        - temp variable for X
71     - lo,hi     - lower and upper bound of the array
72 
73    Output Parameters:
74     + l,r       - of type PetscInt
75 
76    Notes:
77     The TwoWayPartition2/3 variants also partition other arrays along with X.
78     These arrays can have different types, so they provide their own temp t2,t3
79  */
80 #define TwoWayPartitionReverse1(X,pivot,t1,lo,hi,l,r)                            \
81   do {                                                                           \
82     l = lo;                                                                      \
83     r = hi;                                                                      \
84     while (1) {                                                                  \
85       while (X[l] > pivot) l++;                                                  \
86       while (X[r] < pivot) r--;                                                  \
87       if (l >= r) {r++; break;}                                                  \
88       SWAP1(X[l],X[r],t1);                                                       \
89       l++;                                                                       \
90       r--;                                                                       \
91     }                                                                            \
92   } while (0)
93 
94 #define TwoWayPartition2(X,Y,pivot,t1,t2,lo,hi,l,r)                              \
95   do {                                                                           \
96     l = lo;                                                                      \
97     r = hi;                                                                      \
98     while (1) {                                                                  \
99       while (X[l] < pivot) l++;                                                  \
100       while (X[r] > pivot) r--;                                                  \
101       if (l >= r) {r++; break;}                                                  \
102       SWAP2(X[l],X[r],Y[l],Y[r],t1,t2);                                          \
103       l++;                                                                       \
104       r--;                                                                       \
105     }                                                                            \
106   } while (0)
107 
108 #define TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,lo,hi,l,r)                         \
109   do {                                                                           \
110     l = lo;                                                                      \
111     r = hi;                                                                      \
112     while (1) {                                                                  \
113       while (X[l] < pivot) l++;                                                  \
114       while (X[r] > pivot) r--;                                                  \
115       if (l >= r) {r++; break;}                                                  \
116       SWAP3(X[l],X[r],Y[l],Y[r],Z[l],Z[r],t1,t2,t3);                             \
117       l++;                                                                       \
118       r--;                                                                       \
119     }                                                                            \
120   } while (0)
121 
122 /* Templates for similar functions used below */
123 #define QuickSort1(FuncName,X,n,pivot,t1)                                        \
124   do {                                                                           \
125     PetscCount i,j,p,l,r,hi=n-1;                                                 \
126     if (n < 8) {                                                                 \
127       for (i=0; i<n; i++) {                                                      \
128         pivot = X[i];                                                            \
129         for (j=i+1; j<n; j++) {                                                  \
130           if (pivot > X[j]) {                                                    \
131             SWAP1(X[i],X[j],t1);                                                 \
132             pivot = X[i];                                                        \
133           }                                                                      \
134         }                                                                        \
135       }                                                                          \
136     } else {                                                                     \
137       p     = MEDIAN(X,hi);                                                      \
138       pivot = X[p];                                                              \
139       TwoWayPartition1(X,pivot,t1,0,hi,l,r);                                     \
140       PetscCall(FuncName(l,X));                                                    \
141       PetscCall(FuncName(hi-r+1,X+r));                                             \
142     }                                                                            \
143   } while (0)
144 
145 /* Templates for similar functions used below */
146 #define QuickSortReverse1(FuncName,X,n,pivot,t1)                                 \
147   do {                                                                           \
148     PetscCount i,j,p,l,r,hi=n-1;                                                 \
149     if (n < 8) {                                                                 \
150       for (i=0; i<n; i++) {                                                      \
151         pivot = X[i];                                                            \
152         for (j=i+1; j<n; j++) {                                                  \
153           if (pivot < X[j]) {                                                    \
154             SWAP1(X[i],X[j],t1);                                                 \
155             pivot = X[i];                                                        \
156           }                                                                      \
157         }                                                                        \
158       }                                                                          \
159     } else {                                                                     \
160       p     = MEDIAN(X,hi);                                                      \
161       pivot = X[p];                                                              \
162       TwoWayPartitionReverse1(X,pivot,t1,0,hi,l,r);                              \
163       PetscCall(FuncName(l,X));                                                    \
164       PetscCall(FuncName(hi-r+1,X+r));                                             \
165     }                                                                            \
166   } while (0)
167 
168 #define QuickSort2(FuncName,X,Y,n,pivot,t1,t2)                                   \
169   do {                                                                           \
170     PetscCount i,j,p,l,r,hi=n-1;                                                 \
171     if (n < 8) {                                                                 \
172       for (i=0; i<n; i++) {                                                      \
173         pivot = X[i];                                                            \
174         for (j=i+1; j<n; j++) {                                                  \
175           if (pivot > X[j]) {                                                    \
176             SWAP2(X[i],X[j],Y[i],Y[j],t1,t2);                                    \
177             pivot = X[i];                                                        \
178           }                                                                      \
179         }                                                                        \
180       }                                                                          \
181     } else {                                                                     \
182       p     = MEDIAN(X,hi);                                                      \
183       pivot = X[p];                                                              \
184       TwoWayPartition2(X,Y,pivot,t1,t2,0,hi,l,r);                                \
185       PetscCall(FuncName(l,X,Y));                                                  \
186       PetscCall(FuncName(hi-r+1,X+r,Y+r));                                         \
187     }                                                                            \
188   } while (0)
189 
190 #define QuickSort3(FuncName,X,Y,Z,n,pivot,t1,t2,t3)                              \
191   do {                                                                           \
192     PetscCount i,j,p,l,r,hi=n-1;                                                 \
193     if (n < 8) {                                                                 \
194       for (i=0; i<n; i++) {                                                      \
195         pivot = X[i];                                                            \
196         for (j=i+1; j<n; j++) {                                                  \
197           if (pivot > X[j]) {                                                    \
198             SWAP3(X[i],X[j],Y[i],Y[j],Z[i],Z[j],t1,t2,t3);                       \
199             pivot = X[i];                                                        \
200           }                                                                      \
201         }                                                                        \
202       }                                                                          \
203     } else {                                                                     \
204       p     = MEDIAN(X,hi);                                                      \
205       pivot = X[p];                                                              \
206       TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,0,hi,l,r);                           \
207       PetscCall(FuncName(l,X,Y,Z));                                                \
208       PetscCall(FuncName(hi-r+1,X+r,Y+r,Z+r));                                     \
209     }                                                                            \
210   } while (0)
211 
212 /*@
213    PetscSortedInt - Determines whether the array is sorted.
214 
215    Not Collective
216 
217    Input Parameters:
218 +  n  - number of values
219 -  X  - array of integers
220 
221    Output Parameters:
222 .  sorted - flag whether the array is sorted
223 
224    Level: intermediate
225 
226 .seealso: PetscSortInt(), PetscSortedMPIInt(), PetscSortedReal()
227 @*/
228 PetscErrorCode  PetscSortedInt(PetscInt n,const PetscInt X[],PetscBool *sorted)
229 {
230   PetscFunctionBegin;
231   if (n) PetscValidIntPointer(X,2);
232   PetscValidBoolPointer(sorted,3);
233   PetscSorted(n,X,*sorted);
234   PetscFunctionReturn(0);
235 }
236 
237 /*@
238    PetscSortInt - Sorts an array of integers in place in increasing order.
239 
240    Not Collective
241 
242    Input Parameters:
243 +  n  - number of values
244 -  X  - array of integers
245 
246    Notes:
247    This function serves as an alternative to PetscIntSortSemiOrdered(), and may perform faster especially if the array
248    is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
249    code to see which routine is fastest.
250 
251    Level: intermediate
252 
253 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation()
254 @*/
255 PetscErrorCode  PetscSortInt(PetscInt n,PetscInt X[])
256 {
257   PetscInt pivot,t1;
258 
259   PetscFunctionBegin;
260   if (n) PetscValidIntPointer(X,2);
261   QuickSort1(PetscSortInt,X,n,pivot,t1);
262   PetscFunctionReturn(0);
263 }
264 
265 /*@
266    PetscSortReverseInt - Sorts an array of integers in place in decreasing order.
267 
268    Not Collective
269 
270    Input Parameters:
271 +  n  - number of values
272 -  X  - array of integers
273 
274    Level: intermediate
275 
276 .seealso: PetscIntSortSemiOrdered(), PetscSortInt(), PetscSortIntWithPermutation()
277 @*/
278 PetscErrorCode  PetscSortReverseInt(PetscInt n,PetscInt X[])
279 {
280   PetscInt pivot,t1;
281 
282   PetscFunctionBegin;
283   if (n) PetscValidIntPointer(X,2);
284   QuickSortReverse1(PetscSortReverseInt,X,n,pivot,t1);
285   PetscFunctionReturn(0);
286 }
287 
288 /*@
289    PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array
290 
291    Not Collective
292 
293    Input Parameters:
294 +  n  - number of values
295 -  X  - sorted array of integers
296 
297    Output Parameter:
298 .  n - number of non-redundant values
299 
300    Level: intermediate
301 
302 .seealso: PetscSortInt()
303 @*/
304 PetscErrorCode  PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[])
305 {
306   PetscInt i,s = 0,N = *n, b = 0;
307 
308   PetscFunctionBegin;
309   PetscValidIntPointer(n,1);
310   PetscCheckSorted(*n,X);
311   for (i=0; i<N-1; i++) {
312     if (X[b+s+1] != X[b]) {
313       X[b+1] = X[b+s+1]; b++;
314     } else s++;
315   }
316   *n = N - s;
317   PetscFunctionReturn(0);
318 }
319 
320 /*@
321    PetscSortedCheckDupsInt - Checks if a sorted integer array has duplicates
322 
323    Not Collective
324 
325    Input Parameters:
326 +  n  - number of values
327 -  X  - sorted array of integers
328 
329    Output Parameter:
330 .  dups - True if the array has dups, otherwise false
331 
332    Level: intermediate
333 
334 .seealso: PetscSortInt(), PetscCheckDupsInt(), PetscSortRemoveDupsInt(), PetscSortedRemoveDupsInt()
335 @*/
336 PetscErrorCode  PetscSortedCheckDupsInt(PetscInt n,const PetscInt X[],PetscBool *flg)
337 {
338   PetscInt i;
339 
340   PetscFunctionBegin;
341   PetscCheckSorted(n,X);
342   *flg = PETSC_FALSE;
343   for (i=0; i<n-1; i++) {
344     if (X[i+1] == X[i]) {
345       *flg = PETSC_TRUE;
346       break;
347     }
348   }
349   PetscFunctionReturn(0);
350 }
351 
352 /*@
353    PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
354 
355    Not Collective
356 
357    Input Parameters:
358 +  n  - number of values
359 -  X  - array of integers
360 
361    Output Parameter:
362 .  n - number of non-redundant values
363 
364    Level: intermediate
365 
366 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt()
367 @*/
368 PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[])
369 {
370   PetscFunctionBegin;
371   PetscValidIntPointer(n,1);
372   PetscCall(PetscSortInt(*n,X));
373   PetscCall(PetscSortedRemoveDupsInt(n,X));
374   PetscFunctionReturn(0);
375 }
376 
377 /*@
378   PetscFindInt - Finds integer in a sorted array of integers
379 
380    Not Collective
381 
382    Input Parameters:
383 +  key - the integer to locate
384 .  n   - number of values in the array
385 -  X  - array of integers
386 
387    Output Parameter:
388 .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
389 
390    Level: intermediate
391 
392 .seealso: PetscIntSortSemiOrdered(), PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
393 @*/
394 PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt X[], PetscInt *loc)
395 {
396   PetscInt lo = 0,hi = n;
397 
398   PetscFunctionBegin;
399   PetscValidIntPointer(loc,4);
400   if (!n) {*loc = -1; PetscFunctionReturn(0);}
401   PetscValidIntPointer(X,3);
402   PetscCheckSorted(n,X);
403   while (hi - lo > 1) {
404     PetscInt mid = lo + (hi - lo)/2;
405     if (key < X[mid]) hi = mid;
406     else               lo = mid;
407   }
408   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
409   PetscFunctionReturn(0);
410 }
411 
412 /*@
413   PetscCheckDupsInt - Checks if an integer array has duplicates
414 
415    Not Collective
416 
417    Input Parameters:
418 +  n  - number of values in the array
419 -  X  - array of integers
420 
421    Output Parameter:
422 .  dups - True if the array has dups, otherwise false
423 
424    Level: intermediate
425 
426 .seealso: PetscSortRemoveDupsInt(), PetscSortedCheckDupsInt()
427 @*/
428 PetscErrorCode PetscCheckDupsInt(PetscInt n,const PetscInt X[],PetscBool *dups)
429 {
430   PetscInt   i;
431   PetscHSetI ht;
432   PetscBool  missing;
433 
434   PetscFunctionBegin;
435   if (n) PetscValidIntPointer(X,2);
436   PetscValidBoolPointer(dups,3);
437   *dups = PETSC_FALSE;
438   if (n > 1) {
439     PetscCall(PetscHSetICreate(&ht));
440     PetscCall(PetscHSetIResize(ht,n));
441     for (i=0; i<n; i++) {
442       PetscCall(PetscHSetIQueryAdd(ht,X[i],&missing));
443       if (!missing) {*dups = PETSC_TRUE; break;}
444     }
445     PetscCall(PetscHSetIDestroy(&ht));
446   }
447   PetscFunctionReturn(0);
448 }
449 
450 /*@
451   PetscFindMPIInt - Finds MPI integer in a sorted array of integers
452 
453    Not Collective
454 
455    Input Parameters:
456 +  key - the integer to locate
457 .  n   - number of values in the array
458 -  X   - array of integers
459 
460    Output Parameter:
461 .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
462 
463    Level: intermediate
464 
465 .seealso: PetscMPIIntSortSemiOrdered(), PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
466 @*/
467 PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc)
468 {
469   PetscInt lo = 0,hi = n;
470 
471   PetscFunctionBegin;
472   PetscValidIntPointer(loc,4);
473   if (!n) {*loc = -1; PetscFunctionReturn(0);}
474   PetscValidIntPointer(X,3);
475   PetscCheckSorted(n,X);
476   while (hi - lo > 1) {
477     PetscInt mid = lo + (hi - lo)/2;
478     if (key < X[mid]) hi = mid;
479     else               lo = mid;
480   }
481   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
482   PetscFunctionReturn(0);
483 }
484 
485 /*@
486    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
487        changes a second array to match the sorted first array.
488 
489    Not Collective
490 
491    Input Parameters:
492 +  n  - number of values
493 .  X  - array of integers
494 -  Y  - second array of integers
495 
496    Level: intermediate
497 
498 .seealso: PetscIntSortSemiOrderedWithArray(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithCountArray()
499 @*/
500 PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
501 {
502   PetscInt pivot,t1,t2;
503 
504   PetscFunctionBegin;
505   QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2);
506   PetscFunctionReturn(0);
507 }
508 
509 /*@
510    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
511        changes a pair of integer arrays to match the sorted first array.
512 
513    Not Collective
514 
515    Input Parameters:
516 +  n  - number of values
517 .  X  - array of integers
518 .  Y  - second array of integers (first array of the pair)
519 -  Z  - third array of integers  (second array of the pair)
520 
521    Level: intermediate
522 
523 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortIntWithArray(), PetscIntSortSemiOrdered(), PetscSortIntWithIntCountArrayPair()
524 @*/
525 PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
526 {
527   PetscInt pivot,t1,t2,t3;
528 
529   PetscFunctionBegin;
530   QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3);
531   PetscFunctionReturn(0);
532 }
533 
534 /*@
535    PetscSortIntWithCountArray - Sorts an array of integers in place in increasing order;
536        changes a second array to match the sorted first array.
537 
538    Not Collective
539 
540    Input Parameters:
541 +  n  - number of values
542 .  X  - array of integers
543 -  Y  - second array of PetscCounts (signed integers)
544 
545    Level: intermediate
546 
547 .seealso: PetscIntSortSemiOrderedWithArray(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
548 @*/
549 PetscErrorCode  PetscSortIntWithCountArray(PetscCount n,PetscInt X[],PetscCount Y[])
550 {
551   PetscInt   pivot,t1;
552   PetscCount t2;
553 
554   PetscFunctionBegin;
555   QuickSort2(PetscSortIntWithCountArray,X,Y,n,pivot,t1,t2);
556   PetscFunctionReturn(0);
557 }
558 
559 /*@
560    PetscSortIntWithIntCountArrayPair - Sorts an array of integers in place in increasing order;
561        changes an integer array and a PetscCount array to match the sorted first array.
562 
563    Not Collective
564 
565    Input Parameters:
566 +  n  - number of values
567 .  X  - array of integers
568 .  Y  - second array of integers (first array of the pair)
569 -  Z  - third array of PetscCounts  (second array of the pair)
570 
571    Level: intermediate
572 
573    Notes:
574     Usually X, Y are matrix row/column indices, and Z is a permutation array and therefore Z's type is PetscCount to allow 2B+ nonzeros even with 32-bit PetscInt.
575 
576 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray(), PetscIntSortSemiOrdered(), PetscSortIntWithArrayPair()
577 @*/
578 PetscErrorCode  PetscSortIntWithIntCountArrayPair(PetscCount n,PetscInt X[],PetscInt Y[],PetscCount Z[])
579 {
580   PetscInt   pivot,t1,t2; /* pivot is take from X[], so its type is still PetscInt */
581   PetscCount t3; /* temp for Z[] */
582 
583   PetscFunctionBegin;
584   QuickSort3(PetscSortIntWithIntCountArrayPair,X,Y,Z,n,pivot,t1,t2,t3);
585   PetscFunctionReturn(0);
586 }
587 
588 /*@
589    PetscSortedMPIInt - Determines whether the array is sorted.
590 
591    Not Collective
592 
593    Input Parameters:
594 +  n  - number of values
595 -  X  - array of integers
596 
597    Output Parameters:
598 .  sorted - flag whether the array is sorted
599 
600    Level: intermediate
601 
602 .seealso: PetscMPIIntSortSemiOrdered(), PetscSortMPIInt(), PetscSortedInt(), PetscSortedReal()
603 @*/
604 PetscErrorCode  PetscSortedMPIInt(PetscInt n,const PetscMPIInt X[],PetscBool *sorted)
605 {
606   PetscFunctionBegin;
607   PetscSorted(n,X,*sorted);
608   PetscFunctionReturn(0);
609 }
610 
611 /*@
612    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
613 
614    Not Collective
615 
616    Input Parameters:
617 +  n  - number of values
618 -  X  - array of integers
619 
620    Level: intermediate
621 
622    Notes:
623    This function serves as an alternative to PetscMPIIntSortSemiOrdered(), and may perform faster especially if the array
624    is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
625    code to see which routine is fastest.
626 
627 .seealso: PetscMPIIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation()
628 @*/
629 PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
630 {
631   PetscMPIInt pivot,t1;
632 
633   PetscFunctionBegin;
634   QuickSort1(PetscSortMPIInt,X,n,pivot,t1);
635   PetscFunctionReturn(0);
636 }
637 
638 /*@
639    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
640 
641    Not Collective
642 
643    Input Parameters:
644 +  n  - number of values
645 -  X  - array of integers
646 
647    Output Parameter:
648 .  n - number of non-redundant values
649 
650    Level: intermediate
651 
652 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
653 @*/
654 PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
655 {
656   PetscInt s = 0,N = *n,b = 0;
657 
658   PetscFunctionBegin;
659   PetscCall(PetscSortMPIInt(N,X));
660   for (PetscInt i=0; i<N-1; i++) {
661     if (X[b+s+1] != X[b]) {
662       X[b+1] = X[b+s+1]; b++;
663     } else s++;
664   }
665   *n = N - s;
666   PetscFunctionReturn(0);
667 }
668 
669 /*@
670    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
671        changes a second array to match the sorted first array.
672 
673    Not Collective
674 
675    Input Parameters:
676 +  n  - number of values
677 .  X  - array of integers
678 -  Y  - second array of integers
679 
680    Level: intermediate
681 
682 .seealso: PetscMPIIntSortSemiOrderedWithArray(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
683 @*/
684 PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
685 {
686   PetscMPIInt pivot,t1,t2;
687 
688   PetscFunctionBegin;
689   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2);
690   PetscFunctionReturn(0);
691 }
692 
693 /*@
694    PetscSortMPIIntWithIntArray - Sorts an array of MPI integers in place in increasing order;
695        changes a second array of Petsc intergers to match the sorted first array.
696 
697    Not Collective
698 
699    Input Parameters:
700 +  n  - number of values
701 .  X  - array of MPI integers
702 -  Y  - second array of Petsc integers
703 
704    Level: intermediate
705 
706    Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays.
707 
708 .seealso: PetscSortMPIIntWithArray(), PetscIntSortSemiOrderedWithArray(), PetscTimSortWithArray()
709 @*/
710 PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
711 {
712   PetscMPIInt pivot,t1;
713   PetscInt    t2;
714 
715   PetscFunctionBegin;
716   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2);
717   PetscFunctionReturn(0);
718 }
719 
720 /*@
721    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
722        changes a second SCALAR array to match the sorted first INTEGER array.
723 
724    Not Collective
725 
726    Input Parameters:
727 +  n  - number of values
728 .  X  - array of integers
729 -  Y  - second array of scalars
730 
731    Level: intermediate
732 
733 .seealso: PetscTimSortWithArray(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithArray()
734 @*/
735 PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
736 {
737   PetscInt    pivot,t1;
738   PetscScalar t2;
739 
740   PetscFunctionBegin;
741   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2);
742   PetscFunctionReturn(0);
743 }
744 
745 /*@C
746    PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
747        changes a second array to match the sorted first INTEGER array.  Unlike other sort routines, the user must
748        provide workspace (the size of an element in the data array) to use when sorting.
749 
750    Not Collective
751 
752    Input Parameters:
753 +  n  - number of values
754 .  X  - array of integers
755 .  Y  - second array of data
756 .  size - sizeof elements in the data array in bytes
757 -  t2   - workspace of "size" bytes used when sorting
758 
759    Level: intermediate
760 
761 .seealso: PetscTimSortWithArray(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithArray()
762 @*/
763 PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
764 {
765   char     *YY       = (char*)Y;
766   PetscInt  t1,pivot,hi = n-1;
767 
768   PetscFunctionBegin;
769   if (n<8) {
770     for (PetscInt i=0; i<n; i++) {
771       pivot = X[i];
772       for (PetscInt j=i+1; j<n; j++) {
773         if (pivot > X[j]) {
774           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
775           pivot = X[i];
776         }
777       }
778     }
779   } else {
780     /* Two way partition */
781     PetscInt l = 0,r = hi;
782 
783     pivot = X[MEDIAN(X,hi)];
784     while (1) {
785       while (X[l] < pivot) l++;
786       while (X[r] > pivot) r--;
787       if (l >= r) {r++; break;}
788       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
789       l++;
790       r--;
791     }
792     PetscCall(PetscSortIntWithDataArray(l,X,Y,size,t2));
793     PetscCall(PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2));
794   }
795   PetscFunctionReturn(0);
796 }
797 
798 /*@
799    PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.
800 
801    Not Collective
802 
803    Input Parameters:
804 +  an  - number of values in the first array
805 .  aI  - first sorted array of integers
806 .  bn  - number of values in the second array
807 -  bI  - second array of integers
808 
809    Output Parameters:
810 +  n   - number of values in the merged array
811 -  L   - merged sorted array, this is allocated if an array is not provided
812 
813    Level: intermediate
814 
815 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithArray()
816 @*/
817 PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
818 {
819   PetscInt       *L_ = *L, ak, bk, k;
820 
821   PetscFunctionBegin;
822   if (!L_) {
823     PetscCall(PetscMalloc1(an+bn, L));
824     L_   = *L;
825   }
826   k = ak = bk = 0;
827   while (ak < an && bk < bn) {
828     if (aI[ak] == bI[bk]) {
829       L_[k] = aI[ak];
830       ++ak;
831       ++bk;
832       ++k;
833     } else if (aI[ak] < bI[bk]) {
834       L_[k] = aI[ak];
835       ++ak;
836       ++k;
837     } else {
838       L_[k] = bI[bk];
839       ++bk;
840       ++k;
841     }
842   }
843   if (ak < an) {
844     PetscCall(PetscArraycpy(L_+k,aI+ak,an-ak));
845     k   += (an-ak);
846   }
847   if (bk < bn) {
848     PetscCall(PetscArraycpy(L_+k,bI+bk,bn-bk));
849     k   += (bn-bk);
850   }
851   *n = k;
852   PetscFunctionReturn(0);
853 }
854 
855 /*@
856    PetscMergeIntArrayPair -     Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
857                                 The additional arrays are the same length as sorted arrays and are merged
858                                 in the order determined by the merging of the sorted pair.
859 
860    Not Collective
861 
862    Input Parameters:
863 +  an  - number of values in the first array
864 .  aI  - first sorted array of integers
865 .  aJ  - first additional array of integers
866 .  bn  - number of values in the second array
867 .  bI  - second array of integers
868 -  bJ  - second additional of integers
869 
870    Output Parameters:
871 +  n   - number of values in the merged array (== an + bn)
872 .  L   - merged sorted array
873 -  J   - merged additional array
874 
875    Notes:
876     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
877    Level: intermediate
878 
879 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithArray()
880 @*/
881 PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
882 {
883   PetscInt       n_, *L_, *J_, ak, bk, k;
884 
885   PetscFunctionBegin;
886   PetscValidPointer(L,8);
887   PetscValidPointer(J,9);
888   n_ = an + bn;
889   *n = n_;
890   if (!*L) {
891     PetscCall(PetscMalloc1(n_, L));
892   }
893   L_ = *L;
894   if (!*J) {
895     PetscCall(PetscMalloc1(n_, J));
896   }
897   J_   = *J;
898   k = ak = bk = 0;
899   while (ak < an && bk < bn) {
900     if (aI[ak] <= bI[bk]) {
901       L_[k] = aI[ak];
902       J_[k] = aJ[ak];
903       ++ak;
904       ++k;
905     } else {
906       L_[k] = bI[bk];
907       J_[k] = bJ[bk];
908       ++bk;
909       ++k;
910     }
911   }
912   if (ak < an) {
913     PetscCall(PetscArraycpy(L_+k,aI+ak,an-ak));
914     PetscCall(PetscArraycpy(J_+k,aJ+ak,an-ak));
915     k   += (an-ak);
916   }
917   if (bk < bn) {
918     PetscCall(PetscArraycpy(L_+k,bI+bk,bn-bk));
919     PetscCall(PetscArraycpy(J_+k,bJ+bk,bn-bk));
920   }
921   PetscFunctionReturn(0);
922 }
923 
924 /*@
925    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.
926 
927    Not Collective
928 
929    Input Parameters:
930 +  an  - number of values in the first array
931 .  aI  - first sorted array of integers
932 .  bn  - number of values in the second array
933 -  bI  - second array of integers
934 
935    Output Parameters:
936 +  n   - number of values in the merged array (<= an + bn)
937 -  L   - merged sorted array, allocated if address of NULL pointer is passed
938 
939    Level: intermediate
940 
941 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortIntWithArray()
942 @*/
943 PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
944 {
945   PetscInt       ai,bi,k;
946 
947   PetscFunctionBegin;
948   if (!*L) PetscCall(PetscMalloc1((an+bn),L));
949   for (ai=0,bi=0,k=0; ai<an || bi<bn;) {
950     PetscInt t = -1;
951     for (; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
952     for (; bi<bn && bI[bi] == t; bi++);
953     for (; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
954     for (; ai<an && aI[ai] == t; ai++);
955   }
956   *n = k;
957   PetscFunctionReturn(0);
958 }
959 
960 /*@C
961    PetscProcessTree - Prepares tree data to be displayed graphically
962 
963    Not Collective
964 
965    Input Parameters:
966 +  n  - number of values
967 .  mask - indicates those entries in the tree, location 0 is always masked
968 -  parentid - indicates the parent of each entry
969 
970    Output Parameters:
971 +  Nlevels - the number of levels
972 .  Level - for each node tells its level
973 .  Levelcnts - the number of nodes on each level
974 .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
975 -  Column - for each id tells its column index
976 
977    Level: developer
978 
979    Notes:
980     This code is not currently used
981 
982 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
983 @*/
984 PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
985 {
986   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
987   PetscBool      done = PETSC_FALSE;
988 
989   PetscFunctionBegin;
990   PetscCheck(mask[0],PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
991   for (i=0; i<n; i++) {
992     if (mask[i]) continue;
993     PetscCheck(parentid[i]  != i,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
994     PetscCheckFalse(parentid[i] && mask[parentid[i]],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
995   }
996 
997   for (i=0; i<n; i++) {
998     if (!mask[i]) nmask++;
999   }
1000 
1001   /* determine the level in the tree of each node */
1002   PetscCall(PetscCalloc1(n,&level));
1003 
1004   level[0] = 1;
1005   while (!done) {
1006     done = PETSC_TRUE;
1007     for (i=0; i<n; i++) {
1008       if (mask[i]) continue;
1009       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
1010       else if (!level[i]) done = PETSC_FALSE;
1011     }
1012   }
1013   for (i=0; i<n; i++) {
1014     level[i]--;
1015     nlevels = PetscMax(nlevels,level[i]);
1016   }
1017 
1018   /* count the number of nodes on each level and its max */
1019   PetscCall(PetscCalloc1(nlevels,&levelcnt));
1020   for (i=0; i<n; i++) {
1021     if (mask[i]) continue;
1022     levelcnt[level[i]-1]++;
1023   }
1024   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);
1025 
1026   /* for each level sort the ids by the parent id */
1027   PetscCall(PetscMalloc2(levelmax,&workid,levelmax,&workparentid));
1028   PetscCall(PetscMalloc1(nmask,&idbylevel));
1029   for (j=1; j<=nlevels;j++) {
1030     cnt = 0;
1031     for (i=0; i<n; i++) {
1032       if (mask[i]) continue;
1033       if (level[i] != j) continue;
1034       workid[cnt]         = i;
1035       workparentid[cnt++] = parentid[i];
1036     }
1037     /*  PetscIntView(cnt,workparentid,0);
1038     PetscIntView(cnt,workid,0);
1039     PetscCall(PetscSortIntWithArray(cnt,workparentid,workid));
1040     PetscIntView(cnt,workparentid,0);
1041     PetscIntView(cnt,workid,0);*/
1042     PetscCall(PetscArraycpy(idbylevel+tcnt,workid,cnt));
1043     tcnt += cnt;
1044   }
1045   PetscCheck(tcnt == nmask,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
1046   PetscCall(PetscFree2(workid,workparentid));
1047 
1048   /* for each node list its column */
1049   PetscCall(PetscMalloc1(n,&column));
1050   cnt = 0;
1051   for (j=0; j<nlevels; j++) {
1052     for (i=0; i<levelcnt[j]; i++) {
1053       column[idbylevel[cnt++]] = i;
1054     }
1055   }
1056 
1057   *Nlevels   = nlevels;
1058   *Level     = level;
1059   *Levelcnt  = levelcnt;
1060   *Idbylevel = idbylevel;
1061   *Column    = column;
1062   PetscFunctionReturn(0);
1063 }
1064 
1065 /*@
1066   PetscParallelSortedInt - Check whether an integer array, distributed over a communicator, is globally sorted.
1067 
1068   Collective
1069 
1070   Input Parameters:
1071 + comm - the MPI communicator
1072 . n - the local number of integers
1073 - keys - the local array of integers
1074 
1075   Output Parameters:
1076 . is_sorted - whether the array is globally sorted
1077 
1078   Level: developer
1079 
1080 .seealso: PetscParallelSortInt()
1081 @*/
1082 PetscErrorCode PetscParallelSortedInt(MPI_Comm comm, PetscInt n, const PetscInt keys[], PetscBool *is_sorted)
1083 {
1084   PetscBool      sorted;
1085   PetscInt       i, min, max, prevmax;
1086   PetscMPIInt    rank;
1087 
1088   PetscFunctionBegin;
1089   sorted = PETSC_TRUE;
1090   min    = PETSC_MAX_INT;
1091   max    = PETSC_MIN_INT;
1092   if (n) {
1093     min = keys[0];
1094     max = keys[0];
1095   }
1096   for (i = 1; i < n; i++) {
1097     if (keys[i] < keys[i - 1]) break;
1098     min = PetscMin(min,keys[i]);
1099     max = PetscMax(max,keys[i]);
1100   }
1101   if (i < n) sorted = PETSC_FALSE;
1102   prevmax = PETSC_MIN_INT;
1103   PetscCallMPI(MPI_Exscan(&max, &prevmax, 1, MPIU_INT, MPI_MAX, comm));
1104   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1105   if (rank == 0) prevmax = PETSC_MIN_INT;
1106   if (prevmax > min) sorted = PETSC_FALSE;
1107   PetscCallMPI(MPI_Allreduce(&sorted, is_sorted, 1, MPIU_BOOL, MPI_LAND, comm));
1108   PetscFunctionReturn(0);
1109 }
1110