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