Actual source code: sorti.c
petsc-3.6.1 2015-08-06
2: /*
3: This file contains routines for sorting integers. Values are sorted in place.
4: */
5: #include <petsc/private/petscimpl.h> /*I "petscsys.h" I*/
7: #define SWAP(a,b,t) {t=a;a=b;b=t;}
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: /* -----------------------------------------------------------------------*/
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;
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: }
55: /*@
56: PetscSortInt - Sorts an array of integers in place in increasing order.
58: Not Collective
60: Input Parameters:
61: + n - number of values
62: - i - array of integers
64: Level: intermediate
66: Concepts: sorting^ints
68: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
69: @*/
70: PetscErrorCode PetscSortInt(PetscInt n,PetscInt i[])
71: {
72: PetscInt j,k,tmp,ik;
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 PetscSortInt_Private(i,n-1);
86: return(0);
87: }
91: /*@
92: PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
94: Not Collective
96: Input Parameters:
97: + n - number of values
98: - ii - array of integers
100: Output Parameter:
101: . n - number of non-redundant values
103: Level: intermediate
105: Concepts: sorting^ints
107: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
108: @*/
109: PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[])
110: {
112: PetscInt i,s = 0,N = *n, b = 0;
115: PetscSortInt(N,ii);
116: for (i=0; i<N-1; i++) {
117: if (ii[b+s+1] != ii[b]) {
118: ii[b+1] = ii[b+s+1]; b++;
119: } else s++;
120: }
121: *n = N - s;
122: return(0);
123: }
127: /*@
128: PetscFindInt - Finds integer in a sorted array of integers
130: Not Collective
132: Input Parameters:
133: + key - the integer to locate
134: . n - number of values in the array
135: - ii - array of integers
137: Output Parameter:
138: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
140: Level: intermediate
142: Concepts: sorting^ints
144: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
145: @*/
146: PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt ii[], PetscInt *loc)
147: {
148: PetscInt lo = 0,hi = n;
152: if (!n) {*loc = -1; return(0);}
154: while (hi - lo > 1) {
155: PetscInt mid = lo + (hi - lo)/2;
156: if (key < ii[mid]) hi = mid;
157: else lo = mid;
158: }
159: *loc = key == ii[lo] ? lo : -(lo + (key > ii[lo]) + 1);
160: return(0);
161: }
164: /* -----------------------------------------------------------------------*/
165: #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;}
169: /*
170: A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
171: Assumes 0 origin for v, number of elements = right+1 (right is index of
172: right-most member).
173: */
174: static PetscErrorCode PetscSortIntWithArray_Private(PetscInt *v,PetscInt *V,PetscInt right)
175: {
177: PetscInt i,vl,last,tmp;
180: if (right <= 1) {
181: if (right == 1) {
182: if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
183: }
184: return(0);
185: }
186: SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
187: vl = v[0];
188: last = 0;
189: for (i=1; i<=right; i++) {
190: if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
191: }
192: SWAP2(v[0],v[last],V[0],V[last],tmp);
193: PetscSortIntWithArray_Private(v,V,last-1);
194: PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
195: return(0);
196: }
200: /*@
201: PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
202: changes a second array to match the sorted first array.
204: Not Collective
206: Input Parameters:
207: + n - number of values
208: . i - array of integers
209: - I - second array of integers
211: Level: intermediate
213: Concepts: sorting^ints with array
215: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
216: @*/
217: PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[])
218: {
220: PetscInt j,k,tmp,ik;
223: if (n<8) {
224: for (k=0; k<n; k++) {
225: ik = i[k];
226: for (j=k+1; j<n; j++) {
227: if (ik > i[j]) {
228: SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
229: ik = i[k];
230: }
231: }
232: }
233: } else {
234: PetscSortIntWithArray_Private(i,Ii,n-1);
235: }
236: return(0);
237: }
240: #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;}
244: /*
245: A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
246: Assumes 0 origin for v, number of elements = right+1 (right is index of
247: right-most member).
248: */
249: static PetscErrorCode PetscSortIntWithArrayPair_Private(PetscInt *L,PetscInt *J, PetscInt *K, PetscInt right)
250: {
252: PetscInt i,vl,last,tmp;
255: if (right <= 1) {
256: if (right == 1) {
257: if (L[0] > L[1]) SWAP3(L[0],L[1],J[0],J[1],K[0],K[1],tmp);
258: }
259: return(0);
260: }
261: SWAP3(L[0],L[right/2],J[0],J[right/2],K[0],K[right/2],tmp);
262: vl = L[0];
263: last = 0;
264: for (i=1; i<=right; i++) {
265: if (L[i] < vl) {last++; SWAP3(L[last],L[i],J[last],J[i],K[last],K[i],tmp);}
266: }
267: SWAP3(L[0],L[last],J[0],J[last],K[0],K[last],tmp);
268: PetscSortIntWithArrayPair_Private(L,J,K,last-1);
269: PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));
270: return(0);
271: }
275: /*@
276: PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
277: changes a pair of integer arrays to match the sorted first array.
279: Not Collective
281: Input Parameters:
282: + n - number of values
283: . I - array of integers
284: . J - second array of integers (first array of the pair)
285: - K - third array of integers (second array of the pair)
287: Level: intermediate
289: Concepts: sorting^ints with array pair
291: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
292: @*/
293: PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K)
294: {
296: PetscInt j,k,tmp,ik;
299: if (n<8) {
300: for (k=0; k<n; k++) {
301: ik = L[k];
302: for (j=k+1; j<n; j++) {
303: if (ik > L[j]) {
304: SWAP3(L[k],L[j],J[k],J[j],K[k],K[j],tmp);
305: ik = L[k];
306: }
307: }
308: }
309: } else {
310: PetscSortIntWithArrayPair_Private(L,J,K,n-1);
311: }
312: return(0);
313: }
317: /*
318: A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
319: Assumes 0 origin for v, number of elements = right+1 (right is index of
320: right-most member).
321: */
322: static void PetscSortMPIInt_Private(PetscMPIInt *v,PetscInt right)
323: {
324: PetscInt i,j;
325: PetscMPIInt pivot,tmp;
327: if (right <= 1) {
328: if (right == 1) {
329: if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
330: }
331: return;
332: }
333: i = MEDIAN(v,right); /* Choose a pivot */
334: SWAP(v[0],v[i],tmp); /* Move it out of the way */
335: pivot = v[0];
336: for (i=0,j=right+1;; ) {
337: while (++i < j && v[i] <= pivot) ; /* Scan from the left */
338: while (v[--j] > pivot) ; /* Scan from the right */
339: if (i >= j) break;
340: SWAP(v[i],v[j],tmp);
341: }
342: SWAP(v[0],v[j],tmp); /* Put pivot back in place. */
343: PetscSortMPIInt_Private(v,j-1);
344: PetscSortMPIInt_Private(v+j+1,right-(j+1));
345: }
349: /*@
350: PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.
352: Not Collective
354: Input Parameters:
355: + n - number of values
356: - i - array of integers
358: Level: intermediate
360: Concepts: sorting^ints
362: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
363: @*/
364: PetscErrorCode PetscSortMPIInt(PetscInt n,PetscMPIInt i[])
365: {
366: PetscInt j,k;
367: PetscMPIInt tmp,ik;
370: if (n<8) {
371: for (k=0; k<n; k++) {
372: ik = i[k];
373: for (j=k+1; j<n; j++) {
374: if (ik > i[j]) {
375: SWAP(i[k],i[j],tmp);
376: ik = i[k];
377: }
378: }
379: }
380: } else PetscSortMPIInt_Private(i,n-1);
381: return(0);
382: }
386: /*@
387: PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries
389: Not Collective
391: Input Parameters:
392: + n - number of values
393: - ii - array of integers
395: Output Parameter:
396: . n - number of non-redundant values
398: Level: intermediate
400: Concepts: sorting^ints
402: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
403: @*/
404: PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[])
405: {
407: PetscInt i,s = 0,N = *n, b = 0;
410: PetscSortMPIInt(N,ii);
411: for (i=0; i<N-1; i++) {
412: if (ii[b+s+1] != ii[b]) {
413: ii[b+1] = ii[b+s+1]; b++;
414: } else s++;
415: }
416: *n = N - s;
417: return(0);
418: }
422: /*
423: A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
424: Assumes 0 origin for v, number of elements = right+1 (right is index of
425: right-most member).
426: */
427: static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right)
428: {
430: PetscMPIInt i,vl,last,tmp;
433: if (right <= 1) {
434: if (right == 1) {
435: if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
436: }
437: return(0);
438: }
439: SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
440: vl = v[0];
441: last = 0;
442: for (i=1; i<=right; i++) {
443: if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
444: }
445: SWAP2(v[0],v[last],V[0],V[last],tmp);
446: PetscSortMPIIntWithArray_Private(v,V,last-1);
447: PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
448: return(0);
449: }
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: . i - array of integers
462: - I - second array of integers
464: Level: intermediate
466: Concepts: sorting^ints with array
468: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
469: @*/
470: PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[])
471: {
473: PetscMPIInt j,k,tmp,ik;
476: if (n<8) {
477: for (k=0; k<n; k++) {
478: ik = i[k];
479: for (j=k+1; j<n; j++) {
480: if (ik > i[j]) {
481: SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
482: ik = i[k];
483: }
484: }
485: }
486: } else {
487: PetscSortMPIIntWithArray_Private(i,Ii,n-1);
488: }
489: return(0);
490: }
492: /* -----------------------------------------------------------------------*/
493: #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;}
497: /*
498: Modified from PetscSortIntWithArray_Private().
499: */
500: static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right)
501: {
503: PetscInt i,vl,last,tmp;
504: PetscScalar stmp;
507: if (right <= 1) {
508: if (right == 1) {
509: if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp);
510: }
511: return(0);
512: }
513: SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp);
514: vl = v[0];
515: last = 0;
516: for (i=1; i<=right; i++) {
517: if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);}
518: }
519: SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp);
520: PetscSortIntWithScalarArray_Private(v,V,last-1);
521: PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));
522: return(0);
523: }
527: /*@
528: PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
529: changes a second SCALAR array to match the sorted first INTEGER array.
531: Not Collective
533: Input Parameters:
534: + n - number of values
535: . i - array of integers
536: - I - second array of scalars
538: Level: intermediate
540: Concepts: sorting^ints with array
542: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
543: @*/
544: PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[])
545: {
547: PetscInt j,k,tmp,ik;
548: PetscScalar stmp;
551: if (n<8) {
552: for (k=0; k<n; k++) {
553: ik = i[k];
554: for (j=k+1; j<n; j++) {
555: if (ik > i[j]) {
556: SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp);
557: ik = i[k];
558: }
559: }
560: }
561: } else {
562: PetscSortIntWithScalarArray_Private(i,Ii,n-1);
563: }
564: return(0);
565: }
567: #define SWAP2IntData(a,b,c,d,t,td,siz) \
568: do { \
569: PetscErrorCode _ierr; \
570: t=a;a=b;b=t; \
571: _PetscMemcpy(td,c,siz);CHKERRQ(_ierr); \
572: _PetscMemcpy(c,d,siz);CHKERRQ(_ierr); \
573: _PetscMemcpy(d,td,siz);CHKERRQ(_ierr); \
574: } while(0)
578: /*
579: Modified from PetscSortIntWithArray_Private().
580: */
581: static PetscErrorCode PetscSortIntWithDataArray_Private(PetscInt *v,char *V,PetscInt right,size_t size,void *work)
582: {
584: PetscInt i,vl,last,tmp;
587: if (right <= 1) {
588: if (right == 1) {
589: if (v[0] > v[1]) SWAP2IntData(v[0],v[1],V,V+size,tmp,work,size);
590: }
591: return(0);
592: }
593: SWAP2IntData(v[0],v[right/2],V,V+size*(right/2),tmp,work,size);
594: vl = v[0];
595: last = 0;
596: for (i=1; i<=right; i++) {
597: if (v[i] < vl) {last++; SWAP2IntData(v[last],v[i],V+size*last,V+size*i,tmp,work,size);}
598: }
599: SWAP2IntData(v[0],v[last],V,V + size*last,tmp,work,size);
600: PetscSortIntWithDataArray_Private(v,V,last-1,size,work);
601: PetscSortIntWithDataArray_Private(v+last+1,V+size*(last+1),right-(last+1),size,work);
602: return(0);
603: }
607: /*@
608: PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
609: changes a second array to match the sorted first INTEGER array. Unlike other sort routines, the user must
610: provide workspace (the size of an element in the data array) to use when sorting.
612: Not Collective
614: Input Parameters:
615: + n - number of values
616: . i - array of integers
617: . Ii - second array of data
618: . size - sizeof elements in the data array in bytes
619: - work - workspace of "size" bytes used when sorting
621: Level: intermediate
623: Concepts: sorting^ints with array
625: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
626: @*/
627: PetscErrorCode PetscSortIntWithDataArray(PetscInt n,PetscInt i[],void *Ii,size_t size,void *work)
628: {
629: char *V = (char *) Ii;
631: PetscInt j,k,tmp,ik;
634: if (n<8) {
635: for (k=0; k<n; k++) {
636: ik = i[k];
637: for (j=k+1; j<n; j++) {
638: if (ik > i[j]) {
639: SWAP2IntData(i[k],i[j],V+size*k,V+size*j,tmp,work,size);
640: ik = i[k];
641: }
642: }
643: }
644: } else {
645: PetscSortIntWithDataArray_Private(i,V,n-1,size,work);
646: }
647: return(0);
648: }
653: /*@
654: PetscMergeIntArray - Merges two SORTED integer arrays, removes duplicate elements.
656: Not Collective
658: Input Parameters:
659: + an - number of values in the first array
660: . aI - first sorted array of integers
661: . bn - number of values in the second array
662: - bI - second array of integers
664: Output Parameters:
665: + n - number of values in the merged array
666: - I - merged sorted array, this is allocated if an array is not provided
668: Level: intermediate
670: Concepts: merging^arrays
672: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
673: @*/
674: PetscErrorCode PetscMergeIntArray(PetscInt an,const PetscInt *aI, PetscInt bn, const PetscInt *bI, PetscInt *n, PetscInt **L)
675: {
677: PetscInt *L_ = *L, ak, bk, k;
679: if (!L_) {
680: PetscMalloc1(an+bn, L);
681: L_ = *L;
682: }
683: k = ak = bk = 0;
684: while (ak < an && bk < bn) {
685: if (aI[ak] == bI[bk]) {
686: L_[k] = aI[ak];
687: ++ak;
688: ++bk;
689: ++k;
690: } else if (aI[ak] < bI[bk]) {
691: L_[k] = aI[ak];
692: ++ak;
693: ++k;
694: } else {
695: L_[k] = bI[bk];
696: ++bk;
697: ++k;
698: }
699: }
700: if (ak < an) {
701: PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
702: k += (an-ak);
703: }
704: if (bk < bn) {
705: PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
706: k += (bn-bk);
707: }
708: *n = k;
709: return(0);
710: }
714: /*@
715: PetscMergeIntArrayPair - Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
716: The additional arrays are the same length as sorted arrays and are merged
717: in the order determined by the merging of the sorted pair.
719: Not Collective
721: Input Parameters:
722: + an - number of values in the first array
723: . aI - first sorted array of integers
724: . aJ - first additional array of integers
725: . bn - number of values in the second array
726: . bI - second array of integers
727: - bJ - second additional of integers
729: Output Parameters:
730: + n - number of values in the merged array (== an + bn)
731: . I - merged sorted array
732: - J - merged additional array
734: Level: intermediate
736: Concepts: merging^arrays
738: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
739: @*/
740: PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
741: {
743: PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k;
745: n_ = an + bn;
746: *n = n_;
747: if (!L_) {
748: PetscMalloc1(n_, L);
749: L_ = *L;
750: }
751: if (!J_) {
752: PetscMalloc1(n_, &J_);
753: J_ = *J;
754: }
755: k = ak = bk = 0;
756: while (ak < an && bk < bn) {
757: if (aI[ak] <= bI[bk]) {
758: L_[k] = aI[ak];
759: J_[k] = aJ[ak];
760: ++ak;
761: ++k;
762: } else {
763: L_[k] = bI[bk];
764: J_[k] = bJ[bk];
765: ++bk;
766: ++k;
767: }
768: }
769: if (ak < an) {
770: PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
771: PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));
772: k += (an-ak);
773: }
774: if (bk < bn) {
775: PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
776: PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));
777: }
778: return(0);
779: }
784: /*@
785: PetscProcessTree - Prepares tree data to be displayed graphically
787: Not Collective
789: Input Parameters:
790: + n - number of values
791: . mask - indicates those entries in the tree, location 0 is always masked
792: - parentid - indicates the parent of each entry
794: Output Parameters:
795: + Nlevels - the number of levels
796: . Level - for each node tells its level
797: . Levelcnts - the number of nodes on each level
798: . Idbylevel - a list of ids on each of the levels, first level followed by second etc
799: - Column - for each id tells its column index
801: Level: intermediate
804: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
805: @*/
806: PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
807: {
808: PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
810: PetscBool done = PETSC_FALSE;
813: if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
814: for (i=0; i<n; i++) {
815: if (mask[i]) continue;
816: if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
817: if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
818: }
820: for (i=0; i<n; i++) {
821: if (!mask[i]) nmask++;
822: }
824: /* determine the level in the tree of each node */
825: PetscCalloc1(n,&level);
827: level[0] = 1;
828: while (!done) {
829: done = PETSC_TRUE;
830: for (i=0; i<n; i++) {
831: if (mask[i]) continue;
832: if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
833: else if (!level[i]) done = PETSC_FALSE;
834: }
835: }
836: for (i=0; i<n; i++) {
837: level[i]--;
838: nlevels = PetscMax(nlevels,level[i]);
839: }
841: /* count the number of nodes on each level and its max */
842: PetscCalloc1(nlevels,&levelcnt);
843: for (i=0; i<n; i++) {
844: if (mask[i]) continue;
845: levelcnt[level[i]-1]++;
846: }
847: for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);
849: /* for each level sort the ids by the parent id */
850: PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
851: PetscMalloc1(nmask,&idbylevel);
852: for (j=1; j<=nlevels;j++) {
853: cnt = 0;
854: for (i=0; i<n; i++) {
855: if (mask[i]) continue;
856: if (level[i] != j) continue;
857: workid[cnt] = i;
858: workparentid[cnt++] = parentid[i];
859: }
860: /* PetscIntView(cnt,workparentid,0);
861: PetscIntView(cnt,workid,0);
862: PetscSortIntWithArray(cnt,workparentid,workid);
863: PetscIntView(cnt,workparentid,0);
864: PetscIntView(cnt,workid,0);*/
865: PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));
866: tcnt += cnt;
867: }
868: if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
869: PetscFree2(workid,workparentid);
871: /* for each node list its column */
872: PetscMalloc1(n,&column);
873: cnt = 0;
874: for (j=0; j<nlevels; j++) {
875: for (i=0; i<levelcnt[j]; i++) {
876: column[idbylevel[cnt++]] = i;
877: }
878: }
880: *Nlevels = nlevels;
881: *Level = level;
882: *Levelcnt = levelcnt;
883: *Idbylevel = idbylevel;
884: *Column = column;
885: return(0);
886: }