Actual source code: sorti.c
petsc-3.12.5 2020-03-29
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>
7: #include <petsc/private/hashseti.h>
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))))
18: #define MEDIAN(v,right) MEDIAN3(v,right/4,right/2,right/4*3)
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)
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 { \
29: t1=a; a=b; b=t1; \
30: PetscMemcpy(t2,c,siz); \
31: PetscMemcpy(c,d,siz); \
32: PetscMemcpy(d,t2,siz); \
33: } while(0)
35: /*
36: Partition X[lo,hi] into two parts: X[lo,l) <= pivot; X[r,hi] > pivot
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
44: Output Parameters:
45: + l,r - of type PetscInt
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)
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)
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)
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: FuncName(l,X); \
112: FuncName(hi-r+1,X+r); \
113: } \
114: } while(0)
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: FuncName(l,X,Y); \
134: FuncName(hi-r+1,X+r,Y+r); \
135: } \
136: } while(0)
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: FuncName(l,X,Y,Z); \
156: FuncName(hi-r+1,X+r,Y+r,Z+r); \
157: } \
158: } while(0)
160: /*@
161: PetscSortInt - Sorts an array of integers in place in increasing order.
163: Not Collective
165: Input Parameters:
166: + n - number of values
167: - X - array of integers
169: Level: intermediate
171: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
172: @*/
173: PetscErrorCode PetscSortInt(PetscInt n,PetscInt X[])
174: {
176: PetscInt pivot,t1;
179: QuickSort1(PetscSortInt,X,n,pivot,t1,ierr);
180: return(0);
181: }
183: /*@
184: PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array
186: Not Collective
188: Input Parameters:
189: + n - number of values
190: - X - sorted array of integers
192: Output Parameter:
193: . n - number of non-redundant values
195: Level: intermediate
197: .seealso: PetscSortInt()
198: @*/
199: PetscErrorCode PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[])
200: {
201: PetscInt i,s = 0,N = *n, b = 0;
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: return(0);
212: }
214: /*@
215: PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
217: Not Collective
219: Input Parameters:
220: + n - number of values
221: - X - array of integers
223: Output Parameter:
224: . n - number of non-redundant values
226: Level: intermediate
228: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt()
229: @*/
230: PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[])
231: {
235: PetscSortInt(*n,X);
236: PetscSortedRemoveDupsInt(n,X);
237: return(0);
238: }
240: /*@
241: PetscFindInt - Finds integer in a sorted array of integers
243: Not Collective
245: Input Parameters:
246: + key - the integer to locate
247: . n - number of values in the array
248: - X - array of integers
250: Output Parameter:
251: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
253: Level: intermediate
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;
263: if (!n) {*loc = -1; return(0);}
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: return(0);
272: }
274: /*@
277: Not Collective
279: Input Parameters:
280: + n - number of values in the array
281: - X - array of integers
284: Output Parameter:
285: . dups - True if the array has dups, otherwise false
287: Level: intermediate
289: .seealso: PetscSortRemoveDupsInt()
290: @*/
292: {
294: PetscInt i;
295: PetscHSetI ht;
296: PetscBool missing;
300: *dups = PETSC_FALSE;
301: if (n > 1) {
302: PetscHSetICreate(&ht);
303: PetscHSetIResize(ht,n);
304: for (i=0; i<n; i++) {
305: PetscHSetIQueryAdd(ht,X[i],&missing);
306: if(!missing) {*dups = PETSC_TRUE; break;}
307: }
308: PetscHSetIDestroy(&ht);
309: }
310: return(0);
311: }
313: /*@
314: PetscFindMPIInt - Finds MPI integer in a sorted array of integers
316: Not Collective
318: Input Parameters:
319: + key - the integer to locate
320: . n - number of values in the array
321: - X - array of integers
323: Output Parameter:
324: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
326: Level: intermediate
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;
336: if (!n) {*loc = -1; return(0);}
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: return(0);
345: }
347: /*@
348: PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
349: changes a second array to match the sorted first array.
351: Not Collective
353: Input Parameters:
354: + n - number of values
355: . X - array of integers
356: - Y - second array of integers
358: Level: intermediate
360: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
361: @*/
362: PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[])
363: {
365: PetscInt pivot,t1,t2;
368: QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr);
369: return(0);
370: }
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.
376: Not Collective
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)
384: Level: intermediate
386: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
387: @*/
388: PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[])
389: {
391: PetscInt pivot,t1,t2,t3;
394: QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr);
395: return(0);
396: }
398: /*@
399: PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
401: Not Collective
403: Input Parameters:
404: + n - number of values
405: - X - array of integers
407: Level: intermediate
409: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
410: @*/
411: PetscErrorCode PetscSortMPIInt(PetscInt n,PetscMPIInt X[])
412: {
414: PetscMPIInt pivot,t1;
417: QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr);
418: return(0);
419: }
421: /*@
422: PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
424: Not Collective
426: Input Parameters:
427: + n - number of values
428: - X - array of integers
430: Output Parameter:
431: . n - number of non-redundant values
433: Level: intermediate
435: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
436: @*/
437: PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[])
438: {
440: PetscInt i,s = 0,N = *n, b = 0;
443: PetscSortMPIInt(N,X);
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: return(0);
451: }
453: /*@
454: PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
455: changes a second array to match the sorted first array.
457: Not Collective
459: Input Parameters:
460: + n - number of values
461: . X - array of integers
462: - Y - second array of integers
464: Level: intermediate
466: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
467: @*/
468: PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[])
469: {
471: PetscMPIInt pivot,t1,t2;
474: QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr);
475: return(0);
476: }
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.
482: Not Collective
484: Input Parameters:
485: + n - number of values
486: . X - array of MPI integers
487: - Y - second array of Petsc integers
489: Level: intermediate
491: Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays.
493: .seealso: PetscSortMPIIntWithArray()
494: @*/
495: PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[])
496: {
498: PetscMPIInt pivot,t1;
499: PetscInt t2;
502: QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr);
503: return(0);
504: }
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.
510: Not Collective
512: Input Parameters:
513: + n - number of values
514: . X - array of integers
515: - Y - second array of scalars
517: Level: intermediate
519: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
520: @*/
521: PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[])
522: {
524: PetscInt pivot,t1;
525: PetscScalar t2;
528: QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr);
529: return(0);
530: }
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.
537: Not Collective
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
546: Level: intermediate
548: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
549: @*/
550: PetscErrorCode PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2)
551: {
553: char *YY = (char*)Y;
554: PetscInt i,j,p,t1,pivot,hi=n-1,l,r;
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: PetscSortIntWithDataArray(l,X,Y,size,t2);
582: PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);
583: }
584: return(0);
585: }
587: /*@
588: PetscMergeIntArray - Merges two SORTED integer arrays, removes duplicate elements.
590: Not Collective
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
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
602: Level: intermediate
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: {
609: PetscInt *L_ = *L, ak, bk, k;
611: if (!L_) {
612: PetscMalloc1(an+bn, L);
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: PetscArraycpy(L_+k,aI+ak,an-ak);
634: k += (an-ak);
635: }
636: if (bk < bn) {
637: PetscArraycpy(L_+k,bI+bk,bn-bk);
638: k += (bn-bk);
639: }
640: *n = k;
641: return(0);
642: }
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.
649: Not Collective
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
659: Output Parameters:
660: + n - number of values in the merged array (== an + bn)
661: . L - merged sorted array
662: - J - merged additional array
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
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: {
673: PetscInt n_, *L_, *J_, ak, bk, k;
678: n_ = an + bn;
679: *n = n_;
680: if (!*L) {
681: PetscMalloc1(n_, L);
682: }
683: L_ = *L;
684: if (!*J) {
685: PetscMalloc1(n_, J);
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: PetscArraycpy(L_+k,aI+ak,an-ak);
704: PetscArraycpy(J_+k,aJ+ak,an-ak);
705: k += (an-ak);
706: }
707: if (bk < bn) {
708: PetscArraycpy(L_+k,bI+bk,bn-bk);
709: PetscArraycpy(J_+k,bJ+bk,bn-bk);
710: }
711: return(0);
712: }
714: /*@
715: PetscMergeMPIIntArray - Merges two SORTED integer arrays.
717: Not Collective
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
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
729: Level: intermediate
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: {
736: PetscInt ai,bi,k;
739: if (!*L) {PetscMalloc1((an+bn),L);}
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: return(0);
749: }
751: /*@C
752: PetscProcessTree - Prepares tree data to be displayed graphically
754: Not Collective
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
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
768: Level: developer
770: Notes:
771: This code is not currently used
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;
779: PetscBool done = PETSC_FALSE;
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: }
789: for (i=0; i<n; i++) {
790: if (!mask[i]) nmask++;
791: }
793: /* determine the level in the tree of each node */
794: PetscCalloc1(n,&level);
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: }
810: /* count the number of nodes on each level and its max */
811: PetscCalloc1(nlevels,&levelcnt);
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]);
818: /* for each level sort the ids by the parent id */
819: PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
820: PetscMalloc1(nmask,&idbylevel);
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: PetscSortIntWithArray(cnt,workparentid,workid);
832: PetscIntView(cnt,workparentid,0);
833: PetscIntView(cnt,workid,0);*/
834: PetscArraycpy(idbylevel+tcnt,workid,cnt);
835: tcnt += cnt;
836: }
837: if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
838: PetscFree2(workid,workparentid);
840: /* for each node list its column */
841: PetscMalloc1(n,&column);
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: }
849: *Nlevels = nlevels;
850: *Level = level;
851: *Levelcnt = levelcnt;
852: *Idbylevel = idbylevel;
853: *Column = column;
854: return(0);
855: }