xref: /petsc/src/sys/utils/sorti.c (revision 4e278199b78715991f5c71ebbd945c1489263e6c)
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     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       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     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       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     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       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     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       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__ recomended 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(), PetscSortIntPermutation(), PetscSortInt()
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(), PetscSortIntPermutation(), PetscSortIntWithArray(), PetscIntSortSemiOrdered()
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    PetscSortedMPIInt - Determines whether the array is sorted.
505 
506    Not Collective
507 
508    Input Parameters:
509 +  n  - number of values
510 -  X  - array of integers
511 
512    Output Parameters:
513 .  sorted - flag whether the array is sorted
514 
515    Level: intermediate
516 
517 .seealso: PetscMPIIntSortSemiOrdered(), PetscSortMPIInt(), PetscSortedInt(), PetscSortedReal()
518 @*/
519 PetscErrorCode  PetscSortedMPIInt(PetscInt n,const PetscMPIInt X[],PetscBool *sorted)
520 {
521   PetscFunctionBegin;
522   PetscSorted(n,X,*sorted);
523   PetscFunctionReturn(0);
524 }
525 
526 /*@
527    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
528 
529    Not Collective
530 
531    Input Parameters:
532 +  n  - number of values
533 -  X  - array of integers
534 
535    Level: intermediate
536 
537    Notes:
538    This function serves as an alternative to PetscMPIIntSortSemiOrdered(), and may perform faster especially if the array
539    is completely random. There are exceptions to this and so it is __highly__ recomended that the user benchmark their
540    code to see which routine is fastest.
541 
542 .seealso: PetscMPIIntSortSemiOrdered(), PetscSortReal(), PetscSortIntWithPermutation()
543 @*/
544 PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
545 {
546   PetscErrorCode ierr;
547   PetscMPIInt    pivot,t1;
548 
549   PetscFunctionBegin;
550   QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
551   PetscFunctionReturn(0);
552 }
553 
554 /*@
555    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
556 
557    Not Collective
558 
559    Input Parameters:
560 +  n  - number of values
561 -  X  - array of integers
562 
563    Output Parameter:
564 .  n - number of non-redundant values
565 
566    Level: intermediate
567 
568 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
569 @*/
570 PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
571 {
572   PetscErrorCode ierr;
573   PetscInt       i,s = 0,N = *n, b = 0;
574 
575   PetscFunctionBegin;
576   ierr = PetscSortMPIInt(N,X);CHKERRQ(ierr);
577   for (i=0; i<N-1; i++) {
578     if (X[b+s+1] != X[b]) {
579       X[b+1] = X[b+s+1]; b++;
580     } else s++;
581   }
582   *n = N - s;
583   PetscFunctionReturn(0);
584 }
585 
586 /*@
587    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
588        changes a second array to match the sorted first array.
589 
590    Not Collective
591 
592    Input Parameters:
593 +  n  - number of values
594 .  X  - array of integers
595 -  Y  - second array of integers
596 
597    Level: intermediate
598 
599 .seealso: PetscMPIIntSortSemiOrderedWithArray(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
600 @*/
601 PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
602 {
603   PetscErrorCode ierr;
604   PetscMPIInt    pivot,t1,t2;
605 
606   PetscFunctionBegin;
607   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
608   PetscFunctionReturn(0);
609 }
610 
611 /*@
612    PetscSortMPIIntWithIntArray - Sorts an array of MPI integers in place in increasing order;
613        changes a second array of Petsc intergers to match the sorted first array.
614 
615    Not Collective
616 
617    Input Parameters:
618 +  n  - number of values
619 .  X  - array of MPI integers
620 -  Y  - second array of Petsc integers
621 
622    Level: intermediate
623 
624    Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays.
625 
626 .seealso: PetscSortMPIIntWithArray(), PetscIntSortSemiOrderedWithArray(), PetscTimSortWithArray()
627 @*/
628 PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
629 {
630   PetscErrorCode ierr;
631   PetscMPIInt    pivot,t1;
632   PetscInt       t2;
633 
634   PetscFunctionBegin;
635   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
636   PetscFunctionReturn(0);
637 }
638 
639 /*@
640    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
641        changes a second SCALAR array to match the sorted first INTEGER array.
642 
643    Not Collective
644 
645    Input Parameters:
646 +  n  - number of values
647 .  X  - array of integers
648 -  Y  - second array of scalars
649 
650    Level: intermediate
651 
652 .seealso: PetscTimSortWithArray(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
653 @*/
654 PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
655 {
656   PetscErrorCode ierr;
657   PetscInt       pivot,t1;
658   PetscScalar    t2;
659 
660   PetscFunctionBegin;
661   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
662   PetscFunctionReturn(0);
663 }
664 
665 /*@C
666    PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
667        changes a second array to match the sorted first INTEGER array.  Unlike other sort routines, the user must
668        provide workspace (the size of an element in the data array) to use when sorting.
669 
670    Not Collective
671 
672    Input Parameters:
673 +  n  - number of values
674 .  X  - array of integers
675 .  Y  - second array of data
676 .  size - sizeof elements in the data array in bytes
677 -  t2   - workspace of "size" bytes used when sorting
678 
679    Level: intermediate
680 
681 .seealso: PetscTimSortWithArray(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
682 @*/
683 PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
684 {
685   PetscErrorCode ierr;
686   char           *YY = (char*)Y;
687   PetscInt       i,j,p,t1,pivot,hi=n-1,l,r;
688 
689   PetscFunctionBegin;
690   if (n<8) {
691     for (i=0; i<n; i++) {
692       pivot = X[i];
693       for (j=i+1; j<n; j++) {
694         if (pivot > X[j]) {
695           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
696           pivot = X[i];
697         }
698       }
699     }
700   } else {
701     /* Two way partition */
702     p     = MEDIAN(X,hi);
703     pivot = X[p];
704     l     = 0;
705     r     = hi;
706     while (1) {
707       while (X[l] < pivot) l++;
708       while (X[r] > pivot) r--;
709       if (l >= r) {r++; break;}
710       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
711       l++;
712       r--;
713     }
714     ierr = PetscSortIntWithDataArray(l,X,Y,size,t2);CHKERRQ(ierr);
715     ierr = PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);CHKERRQ(ierr);
716   }
717   PetscFunctionReturn(0);
718 }
719 
720 /*@
721    PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.
722 
723    Not Collective
724 
725    Input Parameters:
726 +  an  - number of values in the first array
727 .  aI  - first sorted array of integers
728 .  bn  - number of values in the second array
729 -  bI  - second array of integers
730 
731    Output Parameters:
732 +  n   - number of values in the merged array
733 -  L   - merged sorted array, this is allocated if an array is not provided
734 
735    Level: intermediate
736 
737 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
738 @*/
739 PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
740 {
741   PetscErrorCode ierr;
742   PetscInt       *L_ = *L, ak, bk, k;
743 
744   if (!L_) {
745     ierr = PetscMalloc1(an+bn, L);CHKERRQ(ierr);
746     L_   = *L;
747   }
748   k = ak = bk = 0;
749   while (ak < an && bk < bn) {
750     if (aI[ak] == bI[bk]) {
751       L_[k] = aI[ak];
752       ++ak;
753       ++bk;
754       ++k;
755     } else if (aI[ak] < bI[bk]) {
756       L_[k] = aI[ak];
757       ++ak;
758       ++k;
759     } else {
760       L_[k] = bI[bk];
761       ++bk;
762       ++k;
763     }
764   }
765   if (ak < an) {
766     ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr);
767     k   += (an-ak);
768   }
769   if (bk < bn) {
770     ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr);
771     k   += (bn-bk);
772   }
773   *n = k;
774   PetscFunctionReturn(0);
775 }
776 
777 /*@
778    PetscMergeIntArrayPair -     Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
779                                 The additional arrays are the same length as sorted arrays and are merged
780                                 in the order determined by the merging of the sorted pair.
781 
782    Not Collective
783 
784    Input Parameters:
785 +  an  - number of values in the first array
786 .  aI  - first sorted array of integers
787 .  aJ  - first additional array of integers
788 .  bn  - number of values in the second array
789 .  bI  - second array of integers
790 -  bJ  - second additional of integers
791 
792    Output Parameters:
793 +  n   - number of values in the merged array (== an + bn)
794 .  L   - merged sorted array
795 -  J   - merged additional array
796 
797    Notes:
798     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
799    Level: intermediate
800 
801 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
802 @*/
803 PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
804 {
805   PetscErrorCode ierr;
806   PetscInt       n_, *L_, *J_, ak, bk, k;
807 
808   PetscFunctionBegin;
809   PetscValidIntPointer(L,8);
810   PetscValidIntPointer(J,9);
811   n_ = an + bn;
812   *n = n_;
813   if (!*L) {
814     ierr = PetscMalloc1(n_, L);CHKERRQ(ierr);
815   }
816   L_ = *L;
817   if (!*J) {
818     ierr = PetscMalloc1(n_, J);CHKERRQ(ierr);
819   }
820   J_   = *J;
821   k = ak = bk = 0;
822   while (ak < an && bk < bn) {
823     if (aI[ak] <= bI[bk]) {
824       L_[k] = aI[ak];
825       J_[k] = aJ[ak];
826       ++ak;
827       ++k;
828     } else {
829       L_[k] = bI[bk];
830       J_[k] = bJ[bk];
831       ++bk;
832       ++k;
833     }
834   }
835   if (ak < an) {
836     ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr);
837     ierr = PetscArraycpy(J_+k,aJ+ak,an-ak);CHKERRQ(ierr);
838     k   += (an-ak);
839   }
840   if (bk < bn) {
841     ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr);
842     ierr = PetscArraycpy(J_+k,bJ+bk,bn-bk);CHKERRQ(ierr);
843   }
844   PetscFunctionReturn(0);
845 }
846 
847 /*@
848    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.
849 
850    Not Collective
851 
852    Input Parameters:
853 +  an  - number of values in the first array
854 .  aI  - first sorted array of integers
855 .  bn  - number of values in the second array
856 -  bI  - second array of integers
857 
858    Output Parameters:
859 +  n   - number of values in the merged array (<= an + bn)
860 -  L   - merged sorted array, allocated if address of NULL pointer is passed
861 
862    Level: intermediate
863 
864 .seealso: PetscIntSortSemiOrdered(), PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
865 @*/
866 PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
867 {
868   PetscErrorCode ierr;
869   PetscInt       ai,bi,k;
870 
871   PetscFunctionBegin;
872   if (!*L) {ierr = PetscMalloc1((an+bn),L);CHKERRQ(ierr);}
873   for (ai=0,bi=0,k=0; ai<an || bi<bn;) {
874     PetscInt t = -1;
875     for (; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
876     for (; bi<bn && bI[bi] == t; bi++);
877     for (; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
878     for (; ai<an && aI[ai] == t; ai++);
879   }
880   *n = k;
881   PetscFunctionReturn(0);
882 }
883 
884 /*@C
885    PetscProcessTree - Prepares tree data to be displayed graphically
886 
887    Not Collective
888 
889    Input Parameters:
890 +  n  - number of values
891 .  mask - indicates those entries in the tree, location 0 is always masked
892 -  parentid - indicates the parent of each entry
893 
894    Output Parameters:
895 +  Nlevels - the number of levels
896 .  Level - for each node tells its level
897 .  Levelcnts - the number of nodes on each level
898 .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
899 -  Column - for each id tells its column index
900 
901    Level: developer
902 
903    Notes:
904     This code is not currently used
905 
906 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
907 @*/
908 PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
909 {
910   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
911   PetscErrorCode ierr;
912   PetscBool      done = PETSC_FALSE;
913 
914   PetscFunctionBegin;
915   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
916   for (i=0; i<n; i++) {
917     if (mask[i]) continue;
918     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
919     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
920   }
921 
922   for (i=0; i<n; i++) {
923     if (!mask[i]) nmask++;
924   }
925 
926   /* determine the level in the tree of each node */
927   ierr = PetscCalloc1(n,&level);CHKERRQ(ierr);
928 
929   level[0] = 1;
930   while (!done) {
931     done = PETSC_TRUE;
932     for (i=0; i<n; i++) {
933       if (mask[i]) continue;
934       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
935       else if (!level[i]) done = PETSC_FALSE;
936     }
937   }
938   for (i=0; i<n; i++) {
939     level[i]--;
940     nlevels = PetscMax(nlevels,level[i]);
941   }
942 
943   /* count the number of nodes on each level and its max */
944   ierr = PetscCalloc1(nlevels,&levelcnt);CHKERRQ(ierr);
945   for (i=0; i<n; i++) {
946     if (mask[i]) continue;
947     levelcnt[level[i]-1]++;
948   }
949   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);
950 
951   /* for each level sort the ids by the parent id */
952   ierr = PetscMalloc2(levelmax,&workid,levelmax,&workparentid);CHKERRQ(ierr);
953   ierr = PetscMalloc1(nmask,&idbylevel);CHKERRQ(ierr);
954   for (j=1; j<=nlevels;j++) {
955     cnt = 0;
956     for (i=0; i<n; i++) {
957       if (mask[i]) continue;
958       if (level[i] != j) continue;
959       workid[cnt]         = i;
960       workparentid[cnt++] = parentid[i];
961     }
962     /*  PetscIntView(cnt,workparentid,0);
963     PetscIntView(cnt,workid,0);
964     ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr);
965     PetscIntView(cnt,workparentid,0);
966     PetscIntView(cnt,workid,0);*/
967     ierr  = PetscArraycpy(idbylevel+tcnt,workid,cnt);CHKERRQ(ierr);
968     tcnt += cnt;
969   }
970   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
971   ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr);
972 
973   /* for each node list its column */
974   ierr = PetscMalloc1(n,&column);CHKERRQ(ierr);
975   cnt = 0;
976   for (j=0; j<nlevels; j++) {
977     for (i=0; i<levelcnt[j]; i++) {
978       column[idbylevel[cnt++]] = i;
979     }
980   }
981 
982   *Nlevels   = nlevels;
983   *Level     = level;
984   *Levelcnt  = levelcnt;
985   *Idbylevel = idbylevel;
986   *Column    = column;
987   PetscFunctionReturn(0);
988 }
989 
990 /*@
991   PetscParallelSortedInt - Check whether an integer array, distributed over a communicator, is globally sorted.
992 
993   Collective
994 
995   Input Parameters:
996 + comm - the MPI communicator
997 . n - the local number of integers
998 - keys - the local array of integers
999 
1000   Output Parameters:
1001 . is_sorted - whether the array is globally sorted
1002 
1003   Level: developer
1004 
1005 .seealso: PetscParallelSortInt()
1006 @*/
1007 PetscErrorCode PetscParallelSortedInt(MPI_Comm comm, PetscInt n, const PetscInt keys[], PetscBool *is_sorted)
1008 {
1009   PetscBool      sorted;
1010   PetscInt       i, min, max, prevmax;
1011   PetscMPIInt    rank;
1012   PetscErrorCode ierr;
1013 
1014   PetscFunctionBegin;
1015   sorted = PETSC_TRUE;
1016   min    = PETSC_MAX_INT;
1017   max    = PETSC_MIN_INT;
1018   if (n) {
1019     min = keys[0];
1020     max = keys[0];
1021   }
1022   for (i = 1; i < n; i++) {
1023     if (keys[i] < keys[i - 1]) break;
1024     min = PetscMin(min,keys[i]);
1025     max = PetscMax(max,keys[i]);
1026   }
1027   if (i < n) sorted = PETSC_FALSE;
1028   prevmax = PETSC_MIN_INT;
1029   ierr = MPI_Exscan(&max, &prevmax, 1, MPIU_INT, MPI_MAX, comm);CHKERRMPI(ierr);
1030   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
1031   if (!rank) prevmax = PETSC_MIN_INT;
1032   if (prevmax > min) sorted = PETSC_FALSE;
1033   ierr = MPI_Allreduce(&sorted, is_sorted, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRMPI(ierr);
1034   PetscFunctionReturn(0);
1035 }
1036