Actual source code: sorti.c
petsc-3.4.5 2014-06-29
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: }
571: /*@
572: PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers.
573: The additional arrays are the same length as sorted arrays and are merged
574: in the order determined by the merging of the sorted pair.
576: Not Collective
578: Input Parameters:
579: + an - number of values in the first array
580: . aI - first sorted array of integers
581: . aJ - first additional array of integers
582: . bn - number of values in the second array
583: . bI - second array of integers
584: - bJ - second additional of integers
586: Output Parameters:
587: + n - number of values in the merged array (== an + bn)
588: . I - merged sorted array
589: - J - merged additional array
591: Level: intermediate
593: Concepts: merging^arrays
595: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
596: @*/
597: PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
598: {
600: PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k;
602: n_ = an + bn;
603: *n = n_;
604: if (!L_) {
605: PetscMalloc(n_*sizeof(PetscInt), L);
606: L_ = *L;
607: }
608: if (!J_) {
609: PetscMalloc(n_*sizeof(PetscInt), &J_);
610: J_ = *J;
611: }
612: k = ak = bk = 0;
613: while (ak < an && bk < bn) {
614: if (aI[ak] <= bI[bk]) {
615: L_[k] = aI[ak];
616: J_[k] = aJ[ak];
617: ++ak;
618: ++k;
619: } else {
620: L_[k] = bI[bk];
621: J_[k] = bJ[bk];
622: ++bk;
623: ++k;
624: }
625: }
626: if (ak < an) {
627: PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
628: PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));
629: k += (an-ak);
630: }
631: if (bk < bn) {
632: PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
633: PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));
634: }
635: return(0);
636: }
641: /*@
642: PetscProcessTree - Prepares tree data to be displayed graphically
644: Not Collective
646: Input Parameters:
647: + n - number of values
648: . mask - indicates those entries in the tree, location 0 is always masked
649: - parentid - indicates the parent of each entry
651: Output Parameters:
652: + Nlevels - the number of levels
653: . Level - for each node tells its level
654: . Levelcnts - the number of nodes on each level
655: . Idbylevel - a list of ids on each of the levels, first level followed by second etc
656: - Column - for each id tells its column index
658: Level: intermediate
661: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
662: @*/
663: PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
664: {
665: PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
667: PetscBool done = PETSC_FALSE;
670: if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
671: for (i=0; i<n; i++) {
672: if (mask[i]) continue;
673: if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
674: if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
675: }
677: for (i=0; i<n; i++) {
678: if (!mask[i]) nmask++;
679: }
681: /* determine the level in the tree of each node */
682: PetscMalloc(n*sizeof(PetscInt),&level);
683: PetscMemzero(level,n*sizeof(PetscInt));
685: level[0] = 1;
686: while (!done) {
687: done = PETSC_TRUE;
688: for (i=0; i<n; i++) {
689: if (mask[i]) continue;
690: if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
691: else if (!level[i]) done = PETSC_FALSE;
692: }
693: }
694: for (i=0; i<n; i++) {
695: level[i]--;
696: nlevels = PetscMax(nlevels,level[i]);
697: }
699: /* count the number of nodes on each level and its max */
700: PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);
701: PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));
702: for (i=0; i<n; i++) {
703: if (mask[i]) continue;
704: levelcnt[level[i]-1]++;
705: }
706: for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);
708: /* for each level sort the ids by the parent id */
709: PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);
710: PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);
711: for (j=1; j<=nlevels;j++) {
712: cnt = 0;
713: for (i=0; i<n; i++) {
714: if (mask[i]) continue;
715: if (level[i] != j) continue;
716: workid[cnt] = i;
717: workparentid[cnt++] = parentid[i];
718: }
719: /* PetscIntView(cnt,workparentid,0);
720: PetscIntView(cnt,workid,0);
721: PetscSortIntWithArray(cnt,workparentid,workid);
722: PetscIntView(cnt,workparentid,0);
723: PetscIntView(cnt,workid,0);*/
724: PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));
725: tcnt += cnt;
726: }
727: if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
728: PetscFree2(workid,workparentid);
730: /* for each node list its column */
731: PetscMalloc(n*sizeof(PetscInt),&column);
732: cnt = 0;
733: for (j=0; j<nlevels; j++) {
734: for (i=0; i<levelcnt[j]; i++) {
735: column[idbylevel[cnt++]] = i;
736: }
737: }
739: *Nlevels = nlevels;
740: *Level = level;
741: *Levelcnt = levelcnt;
742: *Idbylevel = idbylevel;
743: *Column = column;
744: return(0);
745: }