xref: /petsc/src/sys/utils/sorti.c (revision 63f3c55c12ae2f190c391f6df6c540efe018a44a)
1 
2 /*
3    This file contains routines for sorting integers. Values are sorted in place.
4    One can use src/sys/examples/tests/ex52.c for benchmarking.
5  */
6 #include <petsc/private/petscimpl.h>                /*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 #define TwoWayPartition2(X,Y,pivot,t1,t2,lo,hi,l,r)                              \
66   do {                                                                           \
67     l = lo;                                                                      \
68     r = hi;                                                                      \
69     while(1) {                                                                   \
70       while (X[l] < pivot) l++;                                                  \
71       while (X[r] > pivot) r--;                                                  \
72       if (l >= r) {r++; break;}                                                  \
73       SWAP2(X[l],X[r],Y[l],Y[r],t1,t2);                                          \
74       l++;                                                                       \
75       r--;                                                                       \
76     }                                                                            \
77   } while(0)
78 
79 #define TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,lo,hi,l,r)                         \
80   do {                                                                           \
81     l = lo;                                                                      \
82     r = hi;                                                                      \
83     while(1) {                                                                   \
84       while (X[l] < pivot) l++;                                                  \
85       while (X[r] > pivot) r--;                                                  \
86       if (l >= r) {r++; break;}                                                  \
87       SWAP3(X[l],X[r],Y[l],Y[r],Z[l],Z[r],t1,t2,t3);                             \
88       l++;                                                                       \
89       r--;                                                                       \
90     }                                                                            \
91   } while(0)
92 
93 /* Templates for similar functions used below */
94 #define QuickSort1(FuncName,X,n,pivot,t1,ierr)                                   \
95   do {                                                                           \
96     PetscInt i,j,p,l,r,hi=n-1;                                                   \
97     if (n < 8) {                                                                 \
98       for (i=0; i<n; i++) {                                                      \
99         pivot = X[i];                                                            \
100         for (j=i+1; j<n; j++) {                                                  \
101           if (pivot > X[j]) {                                                    \
102             SWAP1(X[i],X[j],t1);                                                 \
103             pivot = X[i];                                                        \
104           }                                                                      \
105         }                                                                        \
106       }                                                                          \
107     } else {                                                                     \
108       p     = MEDIAN(X,hi);                                                      \
109       pivot = X[p];                                                              \
110       TwoWayPartition1(X,pivot,t1,0,hi,l,r);                                     \
111       ierr  = FuncName(l,X);CHKERRQ(ierr);                                       \
112       ierr  = FuncName(hi-r+1,X+r);CHKERRQ(ierr);                                \
113     }                                                                            \
114   } while(0)
115 
116 #define QuickSort2(FuncName,X,Y,n,pivot,t1,t2,ierr)                              \
117   do {                                                                           \
118     PetscInt i,j,p,l,r,hi=n-1;                                                   \
119     if (n < 8) {                                                                 \
120       for (i=0; i<n; i++) {                                                      \
121         pivot = X[i];                                                            \
122         for (j=i+1; j<n; j++) {                                                  \
123           if (pivot > X[j]) {                                                    \
124             SWAP2(X[i],X[j],Y[i],Y[j],t1,t2);                                    \
125             pivot = X[i];                                                        \
126           }                                                                      \
127         }                                                                        \
128       }                                                                          \
129     } else {                                                                     \
130       p     = MEDIAN(X,hi);                                                      \
131       pivot = X[p];                                                              \
132       TwoWayPartition2(X,Y,pivot,t1,t2,0,hi,l,r);                                \
133       ierr  = FuncName(l,X,Y);CHKERRQ(ierr);                                     \
134       ierr  = FuncName(hi-r+1,X+r,Y+r);CHKERRQ(ierr);                            \
135     }                                                                            \
136   } while(0)
137 
138 #define QuickSort3(FuncName,X,Y,Z,n,pivot,t1,t2,t3,ierr)                         \
139   do {                                                                           \
140     PetscInt i,j,p,l,r,hi=n-1;                                                   \
141     if (n < 8) {                                                                 \
142       for (i=0; i<n; i++) {                                                      \
143         pivot = X[i];                                                            \
144         for (j=i+1; j<n; j++) {                                                  \
145           if (pivot > X[j]) {                                                    \
146             SWAP3(X[i],X[j],Y[i],Y[j],Z[i],Z[j],t1,t2,t3);                       \
147             pivot = X[i];                                                        \
148           }                                                                      \
149         }                                                                        \
150       }                                                                          \
151     } else {                                                                     \
152       p     = MEDIAN(X,hi);                                                      \
153       pivot = X[p];                                                              \
154       TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,0,hi,l,r);                           \
155       ierr  = FuncName(l,X,Y,Z);CHKERRQ(ierr);                                   \
156       ierr  = FuncName(hi-r+1,X+r,Y+r,Z+r);CHKERRQ(ierr);                        \
157     }                                                                            \
158   } while(0)
159 
160 /*@
161    PetscSortInt - Sorts an array of integers in place in increasing order.
162 
163    Not Collective
164 
165    Input Parameters:
166 +  n  - number of values
167 -  X  - array of integers
168 
169    Level: intermediate
170 
171 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
172 @*/
173 PetscErrorCode  PetscSortInt(PetscInt n,PetscInt X[])
174 {
175   PetscErrorCode ierr;
176   PetscInt       pivot,t1;
177 
178   PetscFunctionBegin;
179   QuickSort1(PetscSortInt,X,n,pivot,t1,ierr);
180   PetscFunctionReturn(0);
181 }
182 
183 /*@
184    PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array
185 
186    Not Collective
187 
188    Input Parameters:
189 +  n  - number of values
190 -  X  - sorted array of integers
191 
192    Output Parameter:
193 .  n - number of non-redundant values
194 
195    Level: intermediate
196 
197 .seealso: PetscSortInt()
198 @*/
199 PetscErrorCode  PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[])
200 {
201   PetscInt i,s = 0,N = *n, b = 0;
202 
203   PetscFunctionBegin;
204   for (i=0; i<N-1; i++) {
205     if (PetscUnlikely(X[b+s+1] < X[b])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Input array is not sorted");
206     if (X[b+s+1] != X[b]) {
207       X[b+1] = X[b+s+1]; b++;
208     } else s++;
209   }
210   *n = N - s;
211   PetscFunctionReturn(0);
212 }
213 
214 /*@
215    PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
216 
217    Not Collective
218 
219    Input Parameters:
220 +  n  - number of values
221 -  X  - array of integers
222 
223    Output Parameter:
224 .  n - number of non-redundant values
225 
226    Level: intermediate
227 
228 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt()
229 @*/
230 PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[])
231 {
232   PetscErrorCode ierr;
233 
234   PetscFunctionBegin;
235   ierr = PetscSortInt(*n,X);CHKERRQ(ierr);
236   ierr = PetscSortedRemoveDupsInt(n,X);CHKERRQ(ierr);
237   PetscFunctionReturn(0);
238 }
239 
240 /*@
241   PetscFindInt - Finds integer in a sorted array of integers
242 
243    Not Collective
244 
245    Input Parameters:
246 +  key - the integer to locate
247 .  n   - number of values in the array
248 -  X  - array of integers
249 
250    Output Parameter:
251 .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
252 
253    Level: intermediate
254 
255 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
256 @*/
257 PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt X[], PetscInt *loc)
258 {
259   PetscInt lo = 0,hi = n;
260 
261   PetscFunctionBegin;
262   PetscValidPointer(loc,4);
263   if (!n) {*loc = -1; PetscFunctionReturn(0);}
264   PetscValidPointer(X,3);
265   while (hi - lo > 1) {
266     PetscInt mid = lo + (hi - lo)/2;
267     if (key < X[mid]) hi = mid;
268     else               lo = mid;
269   }
270   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
271   PetscFunctionReturn(0);
272 }
273 
274 /*@
275   PetscCheckDupsInt - Checks if an integer array has duplicates
276 
277    Not Collective
278 
279    Input Parameters:
280 +  n  - number of values in the array
281 -  X  - array of integers
282 
283 
284    Output Parameter:
285 .  dups - True if the array has dups, otherwise false
286 
287    Level: intermediate
288 
289 .seealso: PetscSortRemoveDupsInt()
290 @*/
291 PetscErrorCode PetscCheckDupsInt(PetscInt n,const PetscInt X[],PetscBool *dups)
292 {
293   PetscErrorCode ierr;
294   PetscInt       i;
295   PetscHSetI     ht;
296   PetscBool      missing;
297 
298   PetscFunctionBegin;
299   PetscValidPointer(dups,3);
300   *dups = PETSC_FALSE;
301   if (n > 1) {
302     ierr = PetscHSetICreate(&ht);CHKERRQ(ierr);
303     ierr = PetscHSetIResize(ht,n);CHKERRQ(ierr);
304     for (i=0; i<n; i++) {
305       ierr = PetscHSetIQueryAdd(ht,X[i],&missing);CHKERRQ(ierr);
306       if(!missing) {*dups = PETSC_TRUE; break;}
307     }
308     ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr);
309   }
310   PetscFunctionReturn(0);
311 }
312 
313 /*@
314   PetscFindMPIInt - Finds MPI integer in a sorted array of integers
315 
316    Not Collective
317 
318    Input Parameters:
319 +  key - the integer to locate
320 .  n   - number of values in the array
321 -  X   - array of integers
322 
323    Output Parameter:
324 .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
325 
326    Level: intermediate
327 
328 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
329 @*/
330 PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc)
331 {
332   PetscInt lo = 0,hi = n;
333 
334   PetscFunctionBegin;
335   PetscValidPointer(loc,4);
336   if (!n) {*loc = -1; PetscFunctionReturn(0);}
337   PetscValidPointer(X,3);
338   while (hi - lo > 1) {
339     PetscInt mid = lo + (hi - lo)/2;
340     if (key < X[mid]) hi = mid;
341     else               lo = mid;
342   }
343   *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1);
344   PetscFunctionReturn(0);
345 }
346 
347 /*@
348    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
349        changes a second array to match the sorted first array.
350 
351    Not Collective
352 
353    Input Parameters:
354 +  n  - number of values
355 .  X  - array of integers
356 -  Y  - second array of integers
357 
358    Level: intermediate
359 
360 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
361 @*/
362 PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
363 {
364   PetscErrorCode ierr;
365   PetscInt       pivot,t1,t2;
366 
367   PetscFunctionBegin;
368   QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr);
369   PetscFunctionReturn(0);
370 }
371 
372 /*@
373    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
374        changes a pair of integer arrays to match the sorted first array.
375 
376    Not Collective
377 
378    Input Parameters:
379 +  n  - number of values
380 .  X  - array of integers
381 .  Y  - second array of integers (first array of the pair)
382 -  Z  - third array of integers  (second array of the pair)
383 
384    Level: intermediate
385 
386 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
387 @*/
388 PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
389 {
390   PetscErrorCode ierr;
391   PetscInt       pivot,t1,t2,t3;
392 
393   PetscFunctionBegin;
394   QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr);
395   PetscFunctionReturn(0);
396 }
397 
398 /*@
399    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
400 
401    Not Collective
402 
403    Input Parameters:
404 +  n  - number of values
405 -  X  - array of integers
406 
407    Level: intermediate
408 
409 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
410 @*/
411 PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
412 {
413   PetscErrorCode ierr;
414   PetscMPIInt    pivot,t1;
415 
416   PetscFunctionBegin;
417   QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
418   PetscFunctionReturn(0);
419 }
420 
421 /*@
422    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
423 
424    Not Collective
425 
426    Input Parameters:
427 +  n  - number of values
428 -  X  - array of integers
429 
430    Output Parameter:
431 .  n - number of non-redundant values
432 
433    Level: intermediate
434 
435 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
436 @*/
437 PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
438 {
439   PetscErrorCode ierr;
440   PetscInt       i,s = 0,N = *n, b = 0;
441 
442   PetscFunctionBegin;
443   ierr = PetscSortMPIInt(N,X);CHKERRQ(ierr);
444   for (i=0; i<N-1; i++) {
445     if (X[b+s+1] != X[b]) {
446       X[b+1] = X[b+s+1]; b++;
447     } else s++;
448   }
449   *n = N - s;
450   PetscFunctionReturn(0);
451 }
452 
453 /*@
454    PetscSortMPIIntWithArray - 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: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
467 @*/
468 PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
469 {
470   PetscErrorCode ierr;
471   PetscMPIInt    pivot,t1,t2;
472 
473   PetscFunctionBegin;
474   QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
475   PetscFunctionReturn(0);
476 }
477 
478 /*@
479    PetscSortMPIIntWithIntArray - Sorts an array of MPI integers in place in increasing order;
480        changes a second array of Petsc intergers to match the sorted first array.
481 
482    Not Collective
483 
484    Input Parameters:
485 +  n  - number of values
486 .  X  - array of MPI integers
487 -  Y  - second array of Petsc integers
488 
489    Level: intermediate
490 
491    Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays.
492 
493 .seealso: PetscSortMPIIntWithArray()
494 @*/
495 PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
496 {
497   PetscErrorCode ierr;
498   PetscMPIInt    pivot,t1;
499   PetscInt       t2;
500 
501   PetscFunctionBegin;
502   QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
503   PetscFunctionReturn(0);
504 }
505 
506 /*@
507    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
508        changes a second SCALAR array to match the sorted first INTEGER array.
509 
510    Not Collective
511 
512    Input Parameters:
513 +  n  - number of values
514 .  X  - array of integers
515 -  Y  - second array of scalars
516 
517    Level: intermediate
518 
519 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
520 @*/
521 PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
522 {
523   PetscErrorCode ierr;
524   PetscInt       pivot,t1;
525   PetscScalar    t2;
526 
527   PetscFunctionBegin;
528   QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
529   PetscFunctionReturn(0);
530 }
531 
532 /*@C
533    PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
534        changes a second array to match the sorted first INTEGER array.  Unlike other sort routines, the user must
535        provide workspace (the size of an element in the data array) to use when sorting.
536 
537    Not Collective
538 
539    Input Parameters:
540 +  n  - number of values
541 .  X  - array of integers
542 .  Y  - second array of data
543 .  size - sizeof elements in the data array in bytes
544 -  t2   - workspace of "size" bytes used when sorting
545 
546    Level: intermediate
547 
548 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
549 @*/
550 PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
551 {
552   PetscErrorCode ierr;
553   char           *YY = (char*)Y;
554   PetscInt       i,j,p,t1,pivot,hi=n-1,l,r;
555 
556   PetscFunctionBegin;
557   if (n<8) {
558     for (i=0; i<n; i++) {
559       pivot = X[i];
560       for (j=i+1; j<n; j++) {
561         if (pivot > X[j]) {
562           SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size);
563           pivot = X[i];
564         }
565       }
566     }
567   } else {
568     /* Two way partition */
569     p     = MEDIAN(X,hi);
570     pivot = X[p];
571     l     = 0;
572     r     = hi;
573     while(1) {
574       while (X[l] < pivot) l++;
575       while (X[r] > pivot) r--;
576       if (l >= r) {r++; break;}
577       SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size);
578       l++;
579       r--;
580     }
581     ierr = PetscSortIntWithDataArray(l,X,Y,size,t2);CHKERRQ(ierr);
582     ierr = PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);CHKERRQ(ierr);
583   }
584   PetscFunctionReturn(0);
585 }
586 
587 /*@
588    PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.
589 
590    Not Collective
591 
592    Input Parameters:
593 +  an  - number of values in the first array
594 .  aI  - first sorted array of integers
595 .  bn  - number of values in the second array
596 -  bI  - second array of integers
597 
598    Output Parameters:
599 +  n   - number of values in the merged array
600 -  L   - merged sorted array, this is allocated if an array is not provided
601 
602    Level: intermediate
603 
604 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
605 @*/
606 PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L)
607 {
608   PetscErrorCode ierr;
609   PetscInt       *L_ = *L, ak, bk, k;
610 
611   if (!L_) {
612     ierr = PetscMalloc1(an+bn, L);CHKERRQ(ierr);
613     L_   = *L;
614   }
615   k = ak = bk = 0;
616   while (ak < an && bk < bn) {
617     if (aI[ak] == bI[bk]) {
618       L_[k] = aI[ak];
619       ++ak;
620       ++bk;
621       ++k;
622     } else if (aI[ak] < bI[bk]) {
623       L_[k] = aI[ak];
624       ++ak;
625       ++k;
626     } else {
627       L_[k] = bI[bk];
628       ++bk;
629       ++k;
630     }
631   }
632   if (ak < an) {
633     ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr);
634     k   += (an-ak);
635   }
636   if (bk < bn) {
637     ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr);
638     k   += (bn-bk);
639   }
640   *n = k;
641   PetscFunctionReturn(0);
642 }
643 
644 /*@
645    PetscMergeIntArrayPair -     Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
646                                 The additional arrays are the same length as sorted arrays and are merged
647                                 in the order determined by the merging of the sorted pair.
648 
649    Not Collective
650 
651    Input Parameters:
652 +  an  - number of values in the first array
653 .  aI  - first sorted array of integers
654 .  aJ  - first additional array of integers
655 .  bn  - number of values in the second array
656 .  bI  - second array of integers
657 -  bJ  - second additional of integers
658 
659    Output Parameters:
660 +  n   - number of values in the merged array (== an + bn)
661 .  L   - merged sorted array
662 -  J   - merged additional array
663 
664    Notes:
665     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
666    Level: intermediate
667 
668 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
669 @*/
670 PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J)
671 {
672   PetscErrorCode ierr;
673   PetscInt       n_, *L_, *J_, ak, bk, k;
674 
675   PetscFunctionBegin;
676   PetscValidIntPointer(L,8);
677   PetscValidIntPointer(J,9);
678   n_ = an + bn;
679   *n = n_;
680   if (!*L) {
681     ierr = PetscMalloc1(n_, L);CHKERRQ(ierr);
682   }
683   L_ = *L;
684   if (!*J) {
685     ierr = PetscMalloc1(n_, J);CHKERRQ(ierr);
686   }
687   J_   = *J;
688   k = ak = bk = 0;
689   while (ak < an && bk < bn) {
690     if (aI[ak] <= bI[bk]) {
691       L_[k] = aI[ak];
692       J_[k] = aJ[ak];
693       ++ak;
694       ++k;
695     } else {
696       L_[k] = bI[bk];
697       J_[k] = bJ[bk];
698       ++bk;
699       ++k;
700     }
701   }
702   if (ak < an) {
703     ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr);
704     ierr = PetscArraycpy(J_+k,aJ+ak,an-ak);CHKERRQ(ierr);
705     k   += (an-ak);
706   }
707   if (bk < bn) {
708     ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr);
709     ierr = PetscArraycpy(J_+k,bJ+bk,bn-bk);CHKERRQ(ierr);
710   }
711   PetscFunctionReturn(0);
712 }
713 
714 /*@
715    PetscMergeMPIIntArray -     Merges two SORTED integer arrays.
716 
717    Not Collective
718 
719    Input Parameters:
720 +  an  - number of values in the first array
721 .  aI  - first sorted array of integers
722 .  bn  - number of values in the second array
723 -  bI  - second array of integers
724 
725    Output Parameters:
726 +  n   - number of values in the merged array (<= an + bn)
727 -  L   - merged sorted array, allocated if address of NULL pointer is passed
728 
729    Level: intermediate
730 
731 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
732 @*/
733 PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L)
734 {
735   PetscErrorCode ierr;
736   PetscInt       ai,bi,k;
737 
738   PetscFunctionBegin;
739   if (!*L) {ierr = PetscMalloc1((an+bn),L);CHKERRQ(ierr);}
740   for (ai=0,bi=0,k=0; ai<an || bi<bn; ) {
741     PetscInt t = -1;
742     for ( ; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai];
743     for ( ; bi<bn && bI[bi] == t; bi++);
744     for ( ; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi];
745     for ( ; ai<an && aI[ai] == t; ai++);
746   }
747   *n = k;
748   PetscFunctionReturn(0);
749 }
750 
751 /*@C
752    PetscProcessTree - Prepares tree data to be displayed graphically
753 
754    Not Collective
755 
756    Input Parameters:
757 +  n  - number of values
758 .  mask - indicates those entries in the tree, location 0 is always masked
759 -  parentid - indicates the parent of each entry
760 
761    Output Parameters:
762 +  Nlevels - the number of levels
763 .  Level - for each node tells its level
764 .  Levelcnts - the number of nodes on each level
765 .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
766 -  Column - for each id tells its column index
767 
768    Level: developer
769 
770    Notes:
771     This code is not currently used
772 
773 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
774 @*/
775 PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
776 {
777   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
778   PetscErrorCode ierr;
779   PetscBool      done = PETSC_FALSE;
780 
781   PetscFunctionBegin;
782   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
783   for (i=0; i<n; i++) {
784     if (mask[i]) continue;
785     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
786     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
787   }
788 
789   for (i=0; i<n; i++) {
790     if (!mask[i]) nmask++;
791   }
792 
793   /* determine the level in the tree of each node */
794   ierr = PetscCalloc1(n,&level);CHKERRQ(ierr);
795 
796   level[0] = 1;
797   while (!done) {
798     done = PETSC_TRUE;
799     for (i=0; i<n; i++) {
800       if (mask[i]) continue;
801       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
802       else if (!level[i]) done = PETSC_FALSE;
803     }
804   }
805   for (i=0; i<n; i++) {
806     level[i]--;
807     nlevels = PetscMax(nlevels,level[i]);
808   }
809 
810   /* count the number of nodes on each level and its max */
811   ierr = PetscCalloc1(nlevels,&levelcnt);CHKERRQ(ierr);
812   for (i=0; i<n; i++) {
813     if (mask[i]) continue;
814     levelcnt[level[i]-1]++;
815   }
816   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);
817 
818   /* for each level sort the ids by the parent id */
819   ierr = PetscMalloc2(levelmax,&workid,levelmax,&workparentid);CHKERRQ(ierr);
820   ierr = PetscMalloc1(nmask,&idbylevel);CHKERRQ(ierr);
821   for (j=1; j<=nlevels;j++) {
822     cnt = 0;
823     for (i=0; i<n; i++) {
824       if (mask[i]) continue;
825       if (level[i] != j) continue;
826       workid[cnt]         = i;
827       workparentid[cnt++] = parentid[i];
828     }
829     /*  PetscIntView(cnt,workparentid,0);
830     PetscIntView(cnt,workid,0);
831     ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr);
832     PetscIntView(cnt,workparentid,0);
833     PetscIntView(cnt,workid,0);*/
834     ierr  = PetscArraycpy(idbylevel+tcnt,workid,cnt);CHKERRQ(ierr);
835     tcnt += cnt;
836   }
837   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
838   ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr);
839 
840   /* for each node list its column */
841   ierr = PetscMalloc1(n,&column);CHKERRQ(ierr);
842   cnt = 0;
843   for (j=0; j<nlevels; j++) {
844     for (i=0; i<levelcnt[j]; i++) {
845       column[idbylevel[cnt++]] = i;
846     }
847   }
848 
849   *Nlevels   = nlevels;
850   *Level     = level;
851   *Levelcnt  = levelcnt;
852   *Idbylevel = idbylevel;
853   *Column    = column;
854   PetscFunctionReturn(0);
855 }
856