xref: /petsc/src/sys/utils/sorti.c (revision b43011057fd0aadbcf9699df876a41721d17264b)
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 
127 /* -----------------------------------------------------------------------*/
128 #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;}
129 
130 #undef __FUNCT__
131 #define __FUNCT__ "PetscSortIntWithArray_Private"
132 /*
133    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
134    Assumes 0 origin for v, number of elements = right+1 (right is index of
135    right-most member).
136 */
137 static PetscErrorCode PetscSortIntWithArray_Private(PetscInt *v,PetscInt *V,PetscInt right)
138 {
139   PetscErrorCode ierr;
140   PetscInt       i,vl,last,tmp;
141 
142   PetscFunctionBegin;
143   if (right <= 1) {
144     if (right == 1) {
145       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
146     }
147     PetscFunctionReturn(0);
148   }
149   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
150   vl   = v[0];
151   last = 0;
152   for (i=1; i<=right; i++) {
153     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
154   }
155   SWAP2(v[0],v[last],V[0],V[last],tmp);
156   ierr = PetscSortIntWithArray_Private(v,V,last-1);CHKERRQ(ierr);
157   ierr = PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
158   PetscFunctionReturn(0);
159 }
160 
161 #undef __FUNCT__
162 #define __FUNCT__ "PetscSortIntWithArray"
163 /*@
164    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
165        changes a second array to match the sorted first array.
166 
167    Not Collective
168 
169    Input Parameters:
170 +  n  - number of values
171 .  i  - array of integers
172 -  I - second array of integers
173 
174    Level: intermediate
175 
176    Concepts: sorting^ints with array
177 
178 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
179 @*/
180 PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[])
181 {
182   PetscErrorCode ierr;
183   PetscInt       j,k,tmp,ik;
184 
185   PetscFunctionBegin;
186   if (n<8) {
187     for (k=0; k<n; k++) {
188       ik = i[k];
189       for (j=k+1; j<n; j++) {
190 	if (ik > i[j]) {
191 	  SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
192 	  ik = i[k];
193 	}
194       }
195     }
196   } else {
197     ierr = PetscSortIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr);
198   }
199   PetscFunctionReturn(0);
200 }
201 
202 
203 #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;}
204 
205 #undef __FUNCT__
206 #define __FUNCT__ "PetscSortIntWithArrayPair_Private"
207 /*
208    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
209    Assumes 0 origin for v, number of elements = right+1 (right is index of
210    right-most member).
211 */
212 static PetscErrorCode PetscSortIntWithArrayPair_Private(PetscInt *L,PetscInt *J, PetscInt *K, PetscInt right)
213 {
214   PetscErrorCode ierr;
215   PetscInt       i,vl,last,tmp;
216 
217   PetscFunctionBegin;
218   if (right <= 1) {
219     if (right == 1) {
220       if (L[0] > L[1]) SWAP3(L[0],L[1],J[0],J[1],K[0],K[1],tmp);
221     }
222     PetscFunctionReturn(0);
223   }
224   SWAP3(L[0],L[right/2],J[0],J[right/2],K[0],K[right/2],tmp);
225   vl   = L[0];
226   last = 0;
227   for (i=1; i<=right; i++) {
228     if (L[i] < vl) {last++; SWAP3(L[last],L[i],J[last],J[i],K[last],K[i],tmp);}
229   }
230   SWAP3(L[0],L[last],J[0],J[last],K[0],K[last],tmp);
231   ierr = PetscSortIntWithArrayPair_Private(L,J,K,last-1);CHKERRQ(ierr);
232   ierr = PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));CHKERRQ(ierr);
233   PetscFunctionReturn(0);
234 }
235 
236 #undef __FUNCT__
237 #define __FUNCT__ "PetscSortIntWithArrayPair"
238 /*@
239    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
240        changes a pair of integer arrays to match the sorted first array.
241 
242    Not Collective
243 
244    Input Parameters:
245 +  n  - number of values
246 .  I  - array of integers
247 .  J  - second array of integers (first array of the pair)
248 -  K  - third array of integers  (second array of the pair)
249 
250    Level: intermediate
251 
252    Concepts: sorting^ints with array pair
253 
254 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
255 @*/
256 PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K)
257 {
258   PetscErrorCode ierr;
259   PetscInt       j,k,tmp,ik;
260 
261   PetscFunctionBegin;
262   if (n<8) {
263     for (k=0; k<n; k++) {
264       ik = L[k];
265       for (j=k+1; j<n; j++) {
266 	if (ik > L[j]) {
267 	  SWAP3(L[k],L[j],J[k],J[j],K[k],K[j],tmp);
268 	  ik = L[k];
269 	}
270       }
271     }
272   } else {
273     ierr = PetscSortIntWithArrayPair_Private(L,J,K,n-1);CHKERRQ(ierr);
274   }
275   PetscFunctionReturn(0);
276 }
277 
278 
279 #undef __FUNCT__
280 #define __FUNCT__ "PetscSortMPIIntWithArray_Private"
281 /*
282    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
283    Assumes 0 origin for v, number of elements = right+1 (right is index of
284    right-most member).
285 */
286 static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right)
287 {
288   PetscErrorCode ierr;
289   PetscMPIInt    i,vl,last,tmp;
290 
291   PetscFunctionBegin;
292   if (right <= 1) {
293     if (right == 1) {
294       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
295     }
296     PetscFunctionReturn(0);
297   }
298   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
299   vl   = v[0];
300   last = 0;
301   for (i=1; i<=right; i++) {
302     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
303   }
304   SWAP2(v[0],v[last],V[0],V[last],tmp);
305   ierr = PetscSortMPIIntWithArray_Private(v,V,last-1);CHKERRQ(ierr);
306   ierr = PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
307   PetscFunctionReturn(0);
308 }
309 
310 #undef __FUNCT__
311 #define __FUNCT__ "PetscSortMPIIntWithArray"
312 /*@
313    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
314        changes a second array to match the sorted first array.
315 
316    Not Collective
317 
318    Input Parameters:
319 +  n  - number of values
320 .  i  - array of integers
321 -  I - second array of integers
322 
323    Level: intermediate
324 
325    Concepts: sorting^ints with array
326 
327 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
328 @*/
329 PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[])
330 {
331   PetscErrorCode ierr;
332   PetscMPIInt    j,k,tmp,ik;
333 
334   PetscFunctionBegin;
335   if (n<8) {
336     for (k=0; k<n; k++) {
337       ik = i[k];
338       for (j=k+1; j<n; j++) {
339 	if (ik > i[j]) {
340 	  SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
341 	  ik = i[k];
342 	}
343       }
344     }
345   } else {
346     ierr = PetscSortMPIIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr);
347   }
348   PetscFunctionReturn(0);
349 }
350 
351 /* -----------------------------------------------------------------------*/
352 #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;}
353 
354 #undef __FUNCT__
355 #define __FUNCT__ "PetscSortIntWithScalarArray_Private"
356 /*
357    Modified from PetscSortIntWithArray_Private().
358 */
359 static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right)
360 {
361   PetscErrorCode ierr;
362   PetscInt       i,vl,last,tmp;
363   PetscScalar    stmp;
364 
365   PetscFunctionBegin;
366   if (right <= 1) {
367     if (right == 1) {
368       if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp);
369     }
370     PetscFunctionReturn(0);
371   }
372   SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp);
373   vl   = v[0];
374   last = 0;
375   for (i=1; i<=right; i++) {
376     if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);}
377   }
378   SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp);
379   ierr = PetscSortIntWithScalarArray_Private(v,V,last-1);CHKERRQ(ierr);
380   ierr = PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
381   PetscFunctionReturn(0);
382 }
383 
384 #undef __FUNCT__
385 #define __FUNCT__ "PetscSortIntWithScalarArray"
386 /*@
387    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
388        changes a second SCALAR array to match the sorted first INTEGER array.
389 
390    Not Collective
391 
392    Input Parameters:
393 +  n  - number of values
394 .  i  - array of integers
395 -  I - second array of scalars
396 
397    Level: intermediate
398 
399    Concepts: sorting^ints with array
400 
401 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
402 @*/
403 PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[])
404 {
405   PetscErrorCode ierr;
406   PetscInt       j,k,tmp,ik;
407   PetscScalar    stmp;
408 
409   PetscFunctionBegin;
410   if (n<8) {
411     for (k=0; k<n; k++) {
412       ik = i[k];
413       for (j=k+1; j<n; j++) {
414 	if (ik > i[j]) {
415 	  SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp);
416 	  ik = i[k];
417 	}
418       }
419     }
420   } else {
421     ierr = PetscSortIntWithScalarArray_Private(i,Ii,n-1);CHKERRQ(ierr);
422   }
423   PetscFunctionReturn(0);
424 }
425 
426 #undef __FUNCT__
427 #define __FUNCT__ "PetscMergeIntArrays"
428 /*@
429    PetscMergeIntArray - Merges n SORTED integer arrays
430 
431    Not Collective
432 
433    Input Parameters:
434 +   n - number of arrays
435 .   lens - length of each array
436 -   arr - pointer to each array
437 
438    Output Parameters:
439 +  N - number of sorted values
440 -  values - sorted values
441 
442    Level: intermediate
443 
444    Concepts: merging^arrays
445 
446    Developer Notes: How do we estimate how many entries there are in the sorted list? When an input arr[i] runs out how do we avoid checking repeatedly again?
447 
448 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
449 @*/
450 PetscErrorCode  PetscMergeIntArrays(PetscInt n,const PetscInt lens[], const PetscInt **arr, PetscInt *N, PetscInt **values)
451 {
452   PetscFunctionBegin;
453   PetscFunctionReturn(0);
454 }
455 
456 #undef __FUNCT__
457 #define __FUNCT__ "PetscMergeIntArrayPair"
458 /*@
459    PetscMergeIntArrayPair -     Merges two SORTED integer arrays along with an additional array of integers.
460                                 The additional arrays are the same length as sorted arrays and are merged
461                                 in the order determined by the merging of the sorted pair.
462 
463    Not Collective
464 
465    Input Parameters:
466 +  an  - number of values in the first array
467 .  aI  - first sorted array of integers
468 .  aJ  - first additional array of integers
469 .  bn  - number of values in the second array
470 .  bI  - second array of integers
471 -  bJ  - second additional of integers
472 
473    Output Parameters:
474 +  n   - number of values in the merged array (== an + bn)
475 .  I   - merged sorted array
476 -  J   - merged additional array
477 
478    Level: intermediate
479 
480    Concepts: merging^arrays
481 
482 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
483 @*/
484 PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
485 {
486   PetscErrorCode ierr;
487   PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k;
488 
489   n_ = an + bn;
490   *n = n_;
491   if(!L_) {
492     ierr = PetscMalloc(n_*sizeof(PetscInt), L); CHKERRQ(ierr);
493     L_ = *L;
494   }
495   if(!J_){
496     ierr = PetscMalloc(n_*sizeof(PetscInt), &J_); CHKERRQ(ierr);
497     J_ = *J;
498   }
499   k = ak = bk = 0;
500   while(ak < an && bk < bn) {
501     if(aI[ak] <= bI[bk]) {
502       L_[k] = aI[ak];
503       J_[k] = aJ[ak];
504       ++ak;
505       ++k;
506     }
507     else {
508       L_[k] = bI[bk];
509       J_[k] = bJ[bk];
510       ++bk;
511       ++k;
512     }
513   }
514   if(ak < an) {
515     ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr);
516     ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr);
517     k += (an-ak);
518   }
519   if(bk < bn) {
520     ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr);
521     ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr);
522     k += (bn-bk);
523   }
524   PetscFunctionReturn(0);
525 }
526 
527 
528 #undef __FUNCT__
529 #define __FUNCT__ "PetscProcessTree"
530 /*@
531    PetscProcessTree - Prepares tree data to be displayed graphically
532 
533    Not Collective
534 
535    Input Parameters:
536 +  n  - number of values
537 .  mask - indicates those entries in the tree, location 0 is always masked
538 -  parentid - indicates the parent of each entry
539 
540    Output Parameters:
541 +  Nlevels - the number of levels
542 .  Level - for each node tells its level
543 .  Levelcnts - the number of nodes on each level
544 .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
545 -  Column - for each id tells its column index
546 
547    Level: intermediate
548 
549 
550 .seealso: PetscSortReal(), PetscSortIntWithPermutation()
551 @*/
552 PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool  mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
553 {
554   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
555   PetscErrorCode ierr;
556   PetscBool      done = PETSC_FALSE;
557 
558   PetscFunctionBegin;
559   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
560   for (i=0; i<n; i++) {
561     if (mask[i]) continue;
562     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
563     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
564   }
565 
566 
567   for (i=0; i<n; i++) {
568     if (!mask[i]) nmask++;
569   }
570 
571   /* determine the level in the tree of each node */
572   ierr = PetscMalloc(n*sizeof(PetscInt),&level);CHKERRQ(ierr);
573   ierr = PetscMemzero(level,n*sizeof(PetscInt));CHKERRQ(ierr);
574   level[0] = 1;
575   while (!done) {
576     done = PETSC_TRUE;
577     for (i=0; i<n; i++) {
578       if (mask[i]) continue;
579       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
580       else if (!level[i]) done = PETSC_FALSE;
581     }
582   }
583   for (i=0; i<n; i++) {
584     level[i]--;
585     nlevels = PetscMax(nlevels,level[i]);
586   }
587 
588   /* count the number of nodes on each level and its max */
589   ierr = PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);CHKERRQ(ierr);
590   ierr = PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));CHKERRQ(ierr);
591   for (i=0; i<n; i++) {
592     if (mask[i]) continue;
593     levelcnt[level[i]-1]++;
594   }
595   for (i=0; i<nlevels;i++) {
596     levelmax = PetscMax(levelmax,levelcnt[i]);
597   }
598 
599   /* for each level sort the ids by the parent id */
600   ierr = PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);CHKERRQ(ierr);
601   ierr = PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);CHKERRQ(ierr);
602   for (j=1; j<=nlevels;j++) {
603     cnt = 0;
604     for (i=0; i<n; i++) {
605       if (mask[i]) continue;
606       if (level[i] != j) continue;
607       workid[cnt]         = i;
608       workparentid[cnt++] = parentid[i];
609     }
610     /*  PetscIntView(cnt,workparentid,0);
611     PetscIntView(cnt,workid,0);
612     ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr);
613     PetscIntView(cnt,workparentid,0);
614     PetscIntView(cnt,workid,0);*/
615     ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr);
616     tcnt += cnt;
617   }
618   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
619   ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr);
620 
621   /* for each node list its column */
622   ierr = PetscMalloc(n*sizeof(PetscInt),&column);CHKERRQ(ierr);
623   cnt = 0;
624   for (j=0; j<nlevels; j++) {
625     for (i=0; i<levelcnt[j]; i++) {
626       column[idbylevel[cnt++]] = i;
627     }
628   }
629 
630   *Nlevels   = nlevels;
631   *Level     = level;
632   *Levelcnt  = levelcnt;
633   *Idbylevel = idbylevel;
634   *Column    = column;
635   PetscFunctionReturn(0);
636 }
637