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