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