xref: /petsc/src/sys/utils/sorti.c (revision 6636e97a9dbd9dfa567780c5ec4fb4018c0e427e)
1 
2 /*
3    This file contains routines for sorting integers. Values are sorted in place.
4  */
5 #include <petscsys.h>                /*I  "petscsys.h"  I*/
6 
7 #define SWAP(a,b,t) {t=a;a=b;b=t;}
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 /* -----------------------------------------------------------------------*/
21 
22 #undef __FUNCT__
23 #define __FUNCT__ "PetscSortInt_Private"
24 /*
25    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
26    Assumes 0 origin for v, number of elements = right+1 (right is index of
27    right-most member).
28 */
29 static void PetscSortInt_Private(PetscInt *v,PetscInt right)
30 {
31   PetscInt       i,j,pivot,tmp;
32 
33   if (right <= 1) {
34     if (right == 1) {
35       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
36     }
37     return;
38   }
39   i = MEDIAN(v,right);          /* Choose a pivot */
40   SWAP(v[0],v[i],tmp);          /* Move it out of the way */
41   pivot = v[0];
42   for (i=0,j=right+1;;) {
43     while (++i < j && v[i] <= pivot); /* Scan from the left */
44     while (v[--j] > pivot);           /* Scan from the right */
45     if (i >= j) break;
46     SWAP(v[i],v[j],tmp);
47   }
48   SWAP(v[0],v[j],tmp);          /* Put pivot back in place. */
49   PetscSortInt_Private(v,j-1);
50   PetscSortInt_Private(v+j+1,right-(j+1));
51 }
52 
53 #undef __FUNCT__
54 #define __FUNCT__ "PetscSortInt"
55 /*@
56    PetscSortInt - Sorts an array of integers in place in increasing order.
57 
58    Not Collective
59 
60    Input Parameters:
61 +  n  - number of values
62 -  i  - array of integers
63 
64    Level: intermediate
65 
66    Concepts: sorting^ints
67 
68 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
69 @*/
70 PetscErrorCode  PetscSortInt(PetscInt n,PetscInt i[])
71 {
72   PetscInt       j,k,tmp,ik;
73 
74   PetscFunctionBegin;
75   if (n<8) {
76     for (k=0; k<n; k++) {
77       ik = i[k];
78       for (j=k+1; j<n; j++) {
79 	if (ik > i[j]) {
80 	  SWAP(i[k],i[j],tmp);
81 	  ik = i[k];
82 	}
83       }
84     }
85   } else {
86     PetscSortInt_Private(i,n-1);
87   }
88   PetscFunctionReturn(0);
89 }
90 
91 #undef __FUNCT__
92 #define __FUNCT__ "PetscSortRemoveDupsInt"
93 /*@
94    PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
95 
96    Not Collective
97 
98    Input Parameters:
99 +  n  - number of values
100 -  ii  - array of integers
101 
102    Output Parameter:
103 .  n - number of non-redundant values
104 
105    Level: intermediate
106 
107    Concepts: sorting^ints
108 
109 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
110 @*/
111 PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[])
112 {
113   PetscErrorCode ierr;
114   PetscInt       i,s = 0,N = *n, b = 0;
115 
116   PetscFunctionBegin;
117   ierr = PetscSortInt(N,ii);CHKERRQ(ierr);
118   for (i=0; i<N-1; i++) {
119     if (ii[b+s+1] != ii[b]) {ii[b+1] = ii[b+s+1]; b++;}
120     else s++;
121   }
122   *n = N - s;
123   PetscFunctionReturn(0);
124 }
125 
126 #undef __FUNCT__
127 #define __FUNCT__ "PetscFindInt"
128 /*@
129   PetscFindInt - Finds and integers in a sorted array of integers
130 
131    Not Collective
132 
133    Input Parameters:
134 +  key - the integer to locate
135 .  n   - number of value in the array
136 -  ii  - array of integers
137 
138    Output Parameter:
139 .  loc - the location, or -1 if the value is not found
140 
141    Level: intermediate
142 
143    Concepts: sorting^ints
144 
145 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
146 @*/
147 PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt ii[], PetscInt *loc)
148 {
149   /* inclusive indices
150      0 <= imin when using truncate toward zero divide
151        imid = (imin+imax)/2;
152      imin unrestricted when using truncate toward minus infinity divide
153        imid = (imin+imax)>>1; or
154        imid = (int)floor((imin+imax)/2.0); */
155   PetscInt imin = 0, imax = n-1;
156 
157   PetscFunctionBegin;
158   PetscValidPointer(ii, 3);
159   PetscValidPointer(loc, 4);
160   /* continually narrow search until just one element remains */
161   while(imin < imax) {
162     PetscInt imid = ((unsigned PetscInt) imin + (unsigned PetscInt) imax) >> 1;
163 
164     /* code must guarantee the interval is reduced at each iteration
165     assert(imid < imax); */
166     /* reduce the search */
167     if (ii[imid] < key) {
168       imin = imid + 1;
169     } else {
170       imax = imid;
171     }
172   }
173   /* At exit of while:
174      if ii[] is empty, then imax < imin
175      otherwise imax == imin */
176   /* deferred test for equality */
177   if ((imax == imin) && (ii[imin] == key)) {
178     *loc = imin;
179   } else {
180     *loc = -1;
181   }
182   PetscFunctionReturn(0);
183 }
184 
185 
186 /* -----------------------------------------------------------------------*/
187 #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;}
188 
189 #undef __FUNCT__
190 #define __FUNCT__ "PetscSortIntWithArray_Private"
191 /*
192    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
193    Assumes 0 origin for v, number of elements = right+1 (right is index of
194    right-most member).
195 */
196 static PetscErrorCode PetscSortIntWithArray_Private(PetscInt *v,PetscInt *V,PetscInt right)
197 {
198   PetscErrorCode ierr;
199   PetscInt       i,vl,last,tmp;
200 
201   PetscFunctionBegin;
202   if (right <= 1) {
203     if (right == 1) {
204       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
205     }
206     PetscFunctionReturn(0);
207   }
208   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
209   vl   = v[0];
210   last = 0;
211   for (i=1; i<=right; i++) {
212     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
213   }
214   SWAP2(v[0],v[last],V[0],V[last],tmp);
215   ierr = PetscSortIntWithArray_Private(v,V,last-1);CHKERRQ(ierr);
216   ierr = PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
217   PetscFunctionReturn(0);
218 }
219 
220 #undef __FUNCT__
221 #define __FUNCT__ "PetscSortIntWithArray"
222 /*@
223    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
224        changes a second array to match the sorted first array.
225 
226    Not Collective
227 
228    Input Parameters:
229 +  n  - number of values
230 .  i  - array of integers
231 -  I - second array of integers
232 
233    Level: intermediate
234 
235    Concepts: sorting^ints with array
236 
237 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
238 @*/
239 PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[])
240 {
241   PetscErrorCode ierr;
242   PetscInt       j,k,tmp,ik;
243 
244   PetscFunctionBegin;
245   if (n<8) {
246     for (k=0; k<n; k++) {
247       ik = i[k];
248       for (j=k+1; j<n; j++) {
249 	if (ik > i[j]) {
250 	  SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
251 	  ik = i[k];
252 	}
253       }
254     }
255   } else {
256     ierr = PetscSortIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr);
257   }
258   PetscFunctionReturn(0);
259 }
260 
261 
262 #define SWAP3(a,b,c,d,e,f,t) {t=a;a=b;b=t;t=c;c=d;d=t;t=e;e=f;f=t;}
263 
264 #undef __FUNCT__
265 #define __FUNCT__ "PetscSortIntWithArrayPair_Private"
266 /*
267    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
268    Assumes 0 origin for v, number of elements = right+1 (right is index of
269    right-most member).
270 */
271 static PetscErrorCode PetscSortIntWithArrayPair_Private(PetscInt *L,PetscInt *J, PetscInt *K, PetscInt right)
272 {
273   PetscErrorCode ierr;
274   PetscInt       i,vl,last,tmp;
275 
276   PetscFunctionBegin;
277   if (right <= 1) {
278     if (right == 1) {
279       if (L[0] > L[1]) SWAP3(L[0],L[1],J[0],J[1],K[0],K[1],tmp);
280     }
281     PetscFunctionReturn(0);
282   }
283   SWAP3(L[0],L[right/2],J[0],J[right/2],K[0],K[right/2],tmp);
284   vl   = L[0];
285   last = 0;
286   for (i=1; i<=right; i++) {
287     if (L[i] < vl) {last++; SWAP3(L[last],L[i],J[last],J[i],K[last],K[i],tmp);}
288   }
289   SWAP3(L[0],L[last],J[0],J[last],K[0],K[last],tmp);
290   ierr = PetscSortIntWithArrayPair_Private(L,J,K,last-1);CHKERRQ(ierr);
291   ierr = PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));CHKERRQ(ierr);
292   PetscFunctionReturn(0);
293 }
294 
295 #undef __FUNCT__
296 #define __FUNCT__ "PetscSortIntWithArrayPair"
297 /*@
298    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
299        changes a pair of integer arrays to match the sorted first array.
300 
301    Not Collective
302 
303    Input Parameters:
304 +  n  - number of values
305 .  I  - array of integers
306 .  J  - second array of integers (first array of the pair)
307 -  K  - third array of integers  (second array of the pair)
308 
309    Level: intermediate
310 
311    Concepts: sorting^ints with array pair
312 
313 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
314 @*/
315 PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K)
316 {
317   PetscErrorCode ierr;
318   PetscInt       j,k,tmp,ik;
319 
320   PetscFunctionBegin;
321   if (n<8) {
322     for (k=0; k<n; k++) {
323       ik = L[k];
324       for (j=k+1; j<n; j++) {
325 	if (ik > L[j]) {
326 	  SWAP3(L[k],L[j],J[k],J[j],K[k],K[j],tmp);
327 	  ik = L[k];
328 	}
329       }
330     }
331   } else {
332     ierr = PetscSortIntWithArrayPair_Private(L,J,K,n-1);CHKERRQ(ierr);
333   }
334   PetscFunctionReturn(0);
335 }
336 
337 #undef __FUNCT__
338 #define __FUNCT__ "PetscSortMPIInt_Private"
339 /*
340    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
341    Assumes 0 origin for v, number of elements = right+1 (right is index of
342    right-most member).
343 */
344 static void PetscSortMPIInt_Private(PetscMPIInt *v,PetscInt right)
345 {
346   PetscInt          i,j;
347   PetscMPIInt       pivot,tmp;
348 
349   if (right <= 1) {
350     if (right == 1) {
351       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
352     }
353     return;
354   }
355   i = MEDIAN(v,right);          /* Choose a pivot */
356   SWAP(v[0],v[i],tmp);          /* Move it out of the way */
357   pivot = v[0];
358   for (i=0,j=right+1;;) {
359     while (++i < j && v[i] <= pivot); /* Scan from the left */
360     while (v[--j] > pivot);           /* Scan from the right */
361     if (i >= j) break;
362     SWAP(v[i],v[j],tmp);
363   }
364   SWAP(v[0],v[j],tmp);          /* Put pivot back in place. */
365   PetscSortMPIInt_Private(v,j-1);
366   PetscSortMPIInt_Private(v+j+1,right-(j+1));
367 }
368 
369 #undef __FUNCT__
370 #define __FUNCT__ "PetscSortMPIInt"
371 /*@
372    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
373 
374    Not Collective
375 
376    Input Parameters:
377 +  n  - number of values
378 -  i  - array of integers
379 
380    Level: intermediate
381 
382    Concepts: sorting^ints
383 
384 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
385 @*/
386 PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt i[])
387 {
388   PetscInt       j,k;
389   PetscMPIInt    tmp,ik;
390 
391   PetscFunctionBegin;
392   if (n<8) {
393     for (k=0; k<n; k++) {
394       ik = i[k];
395       for (j=k+1; j<n; j++) {
396 	if (ik > i[j]) {
397 	  SWAP(i[k],i[j],tmp);
398 	  ik = i[k];
399 	}
400       }
401     }
402   } else {
403     PetscSortMPIInt_Private(i,n-1);
404   }
405   PetscFunctionReturn(0);
406 }
407 
408 #undef __FUNCT__
409 #define __FUNCT__ "PetscSortRemoveDupsMPIInt"
410 /*@
411    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
412 
413    Not Collective
414 
415    Input Parameters:
416 +  n  - number of values
417 -  ii  - array of integers
418 
419    Output Parameter:
420 .  n - number of non-redundant values
421 
422    Level: intermediate
423 
424    Concepts: sorting^ints
425 
426 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
427 @*/
428 PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[])
429 {
430   PetscErrorCode ierr;
431   PetscInt       i,s = 0,N = *n, b = 0;
432 
433   PetscFunctionBegin;
434   ierr = PetscSortMPIInt(N,ii);CHKERRQ(ierr);
435   for (i=0; i<N-1; i++) {
436     if (ii[b+s+1] != ii[b]) {ii[b+1] = ii[b+s+1]; b++;}
437     else s++;
438   }
439   *n = N - s;
440   PetscFunctionReturn(0);
441 }
442 
443 #undef __FUNCT__
444 #define __FUNCT__ "PetscSortMPIIntWithArray_Private"
445 /*
446    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
447    Assumes 0 origin for v, number of elements = right+1 (right is index of
448    right-most member).
449 */
450 static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right)
451 {
452   PetscErrorCode ierr;
453   PetscMPIInt    i,vl,last,tmp;
454 
455   PetscFunctionBegin;
456   if (right <= 1) {
457     if (right == 1) {
458       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
459     }
460     PetscFunctionReturn(0);
461   }
462   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
463   vl   = v[0];
464   last = 0;
465   for (i=1; i<=right; i++) {
466     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
467   }
468   SWAP2(v[0],v[last],V[0],V[last],tmp);
469   ierr = PetscSortMPIIntWithArray_Private(v,V,last-1);CHKERRQ(ierr);
470   ierr = PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
471   PetscFunctionReturn(0);
472 }
473 
474 #undef __FUNCT__
475 #define __FUNCT__ "PetscSortMPIIntWithArray"
476 /*@
477    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
478        changes a second array to match the sorted first array.
479 
480    Not Collective
481 
482    Input Parameters:
483 +  n  - number of values
484 .  i  - array of integers
485 -  I - second array of integers
486 
487    Level: intermediate
488 
489    Concepts: sorting^ints with array
490 
491 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
492 @*/
493 PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[])
494 {
495   PetscErrorCode ierr;
496   PetscMPIInt    j,k,tmp,ik;
497 
498   PetscFunctionBegin;
499   if (n<8) {
500     for (k=0; k<n; k++) {
501       ik = i[k];
502       for (j=k+1; j<n; j++) {
503 	if (ik > i[j]) {
504 	  SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
505 	  ik = i[k];
506 	}
507       }
508     }
509   } else {
510     ierr = PetscSortMPIIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr);
511   }
512   PetscFunctionReturn(0);
513 }
514 
515 /* -----------------------------------------------------------------------*/
516 #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;}
517 
518 #undef __FUNCT__
519 #define __FUNCT__ "PetscSortIntWithScalarArray_Private"
520 /*
521    Modified from PetscSortIntWithArray_Private().
522 */
523 static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right)
524 {
525   PetscErrorCode ierr;
526   PetscInt       i,vl,last,tmp;
527   PetscScalar    stmp;
528 
529   PetscFunctionBegin;
530   if (right <= 1) {
531     if (right == 1) {
532       if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp);
533     }
534     PetscFunctionReturn(0);
535   }
536   SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp);
537   vl   = v[0];
538   last = 0;
539   for (i=1; i<=right; i++) {
540     if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);}
541   }
542   SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp);
543   ierr = PetscSortIntWithScalarArray_Private(v,V,last-1);CHKERRQ(ierr);
544   ierr = PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
545   PetscFunctionReturn(0);
546 }
547 
548 #undef __FUNCT__
549 #define __FUNCT__ "PetscSortIntWithScalarArray"
550 /*@
551    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
552        changes a second SCALAR array to match the sorted first INTEGER array.
553 
554    Not Collective
555 
556    Input Parameters:
557 +  n  - number of values
558 .  i  - array of integers
559 -  I - second array of scalars
560 
561    Level: intermediate
562 
563    Concepts: sorting^ints with array
564 
565 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
566 @*/
567 PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[])
568 {
569   PetscErrorCode ierr;
570   PetscInt       j,k,tmp,ik;
571   PetscScalar    stmp;
572 
573   PetscFunctionBegin;
574   if (n<8) {
575     for (k=0; k<n; k++) {
576       ik = i[k];
577       for (j=k+1; j<n; j++) {
578 	if (ik > i[j]) {
579 	  SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp);
580 	  ik = i[k];
581 	}
582       }
583     }
584   } else {
585     ierr = PetscSortIntWithScalarArray_Private(i,Ii,n-1);CHKERRQ(ierr);
586   }
587   PetscFunctionReturn(0);
588 }
589 
590 
591 
592 #undef __FUNCT__
593 #define __FUNCT__ "PetscMergeIntArrayPair"
594 /*@
595    PetscMergeIntArrayPair -     Merges two SORTED integer arrays along with an additional array of integers.
596                                 The additional arrays are the same length as sorted arrays and are merged
597                                 in the order determined by the merging of the sorted pair.
598 
599    Not Collective
600 
601    Input Parameters:
602 +  an  - number of values in the first array
603 .  aI  - first sorted array of integers
604 .  aJ  - first additional array of integers
605 .  bn  - number of values in the second array
606 .  bI  - second array of integers
607 -  bJ  - second additional of integers
608 
609    Output Parameters:
610 +  n   - number of values in the merged array (== an + bn)
611 .  I   - merged sorted array
612 -  J   - merged additional array
613 
614    Level: intermediate
615 
616    Concepts: merging^arrays
617 
618 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
619 @*/
620 PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
621 {
622   PetscErrorCode ierr;
623   PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k;
624 
625   n_ = an + bn;
626   *n = n_;
627   if (!L_) {
628     ierr = PetscMalloc(n_*sizeof(PetscInt), L); CHKERRQ(ierr);
629     L_ = *L;
630   }
631   if (!J_){
632     ierr = PetscMalloc(n_*sizeof(PetscInt), &J_); CHKERRQ(ierr);
633     J_ = *J;
634   }
635   k = ak = bk = 0;
636   while(ak < an && bk < bn) {
637     if (aI[ak] <= bI[bk]) {
638       L_[k] = aI[ak];
639       J_[k] = aJ[ak];
640       ++ak;
641       ++k;
642     }
643     else {
644       L_[k] = bI[bk];
645       J_[k] = bJ[bk];
646       ++bk;
647       ++k;
648     }
649   }
650   if (ak < an) {
651     ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr);
652     ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr);
653     k += (an-ak);
654   }
655   if (bk < bn) {
656     ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr);
657     ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr);
658     k += (bn-bk);
659   }
660   PetscFunctionReturn(0);
661 }
662 
663 
664 #undef __FUNCT__
665 #define __FUNCT__ "PetscProcessTree"
666 /*@
667    PetscProcessTree - Prepares tree data to be displayed graphically
668 
669    Not Collective
670 
671    Input Parameters:
672 +  n  - number of values
673 .  mask - indicates those entries in the tree, location 0 is always masked
674 -  parentid - indicates the parent of each entry
675 
676    Output Parameters:
677 +  Nlevels - the number of levels
678 .  Level - for each node tells its level
679 .  Levelcnts - the number of nodes on each level
680 .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
681 -  Column - for each id tells its column index
682 
683    Level: intermediate
684 
685 
686 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
687 @*/
688 PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool  mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
689 {
690   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
691   PetscErrorCode ierr;
692   PetscBool      done = PETSC_FALSE;
693 
694   PetscFunctionBegin;
695   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
696   for (i=0; i<n; i++) {
697     if (mask[i]) continue;
698     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
699     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
700   }
701 
702 
703   for (i=0; i<n; i++) {
704     if (!mask[i]) nmask++;
705   }
706 
707   /* determine the level in the tree of each node */
708   ierr = PetscMalloc(n*sizeof(PetscInt),&level);CHKERRQ(ierr);
709   ierr = PetscMemzero(level,n*sizeof(PetscInt));CHKERRQ(ierr);
710   level[0] = 1;
711   while (!done) {
712     done = PETSC_TRUE;
713     for (i=0; i<n; i++) {
714       if (mask[i]) continue;
715       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
716       else if (!level[i]) done = PETSC_FALSE;
717     }
718   }
719   for (i=0; i<n; i++) {
720     level[i]--;
721     nlevels = PetscMax(nlevels,level[i]);
722   }
723 
724   /* count the number of nodes on each level and its max */
725   ierr = PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);CHKERRQ(ierr);
726   ierr = PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));CHKERRQ(ierr);
727   for (i=0; i<n; i++) {
728     if (mask[i]) continue;
729     levelcnt[level[i]-1]++;
730   }
731   for (i=0; i<nlevels;i++) {
732     levelmax = PetscMax(levelmax,levelcnt[i]);
733   }
734 
735   /* for each level sort the ids by the parent id */
736   ierr = PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);CHKERRQ(ierr);
737   ierr = PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);CHKERRQ(ierr);
738   for (j=1; j<=nlevels;j++) {
739     cnt = 0;
740     for (i=0; i<n; i++) {
741       if (mask[i]) continue;
742       if (level[i] != j) continue;
743       workid[cnt]         = i;
744       workparentid[cnt++] = parentid[i];
745     }
746     /*  PetscIntView(cnt,workparentid,0);
747     PetscIntView(cnt,workid,0);
748     ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr);
749     PetscIntView(cnt,workparentid,0);
750     PetscIntView(cnt,workid,0);*/
751     ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr);
752     tcnt += cnt;
753   }
754   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
755   ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr);
756 
757   /* for each node list its column */
758   ierr = PetscMalloc(n*sizeof(PetscInt),&column);CHKERRQ(ierr);
759   cnt = 0;
760   for (j=0; j<nlevels; j++) {
761     for (i=0; i<levelcnt[j]; i++) {
762       column[idbylevel[cnt++]] = i;
763     }
764   }
765 
766   *Nlevels   = nlevels;
767   *Level     = level;
768   *Levelcnt  = levelcnt;
769   *Idbylevel = idbylevel;
770   *Column    = column;
771   PetscFunctionReturn(0);
772 }
773