Actual source code: sorti.c
petsc-3.3-p7 2013-05-11
2: /*
3: This file contains routines for sorting integers. Values are sorted in place.
4: */
5: #include <petscsys.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 {
86: PetscSortInt_Private(i,n-1);
87: }
88: return(0);
89: }
93: /*@
94: PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries
96: Not Collective
98: Input Parameters:
99: + n - number of values
100: - ii - array of integers
102: Output Parameter:
103: . n - number of non-redundant values
105: Level: intermediate
107: Concepts: sorting^ints
109: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
110: @*/
111: PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[])
112: {
114: PetscInt i,s = 0,N = *n, b = 0;
117: PetscSortInt(N,ii);
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: return(0);
124: }
127: /* -----------------------------------------------------------------------*/
128: #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;}
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: {
140: PetscInt i,vl,last,tmp;
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: return(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: PetscSortIntWithArray_Private(v,V,last-1);
157: PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
158: return(0);
159: }
163: /*@
164: PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
165: changes a second array to match the sorted first array.
167: Not Collective
169: Input Parameters:
170: + n - number of values
171: . i - array of integers
172: - I - second array of integers
174: Level: intermediate
176: Concepts: sorting^ints with array
178: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
179: @*/
180: PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[])
181: {
183: PetscInt j,k,tmp,ik;
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: PetscSortIntWithArray_Private(i,Ii,n-1);
198: }
199: return(0);
200: }
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;}
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: {
215: PetscInt i,vl,last,tmp;
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: return(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: PetscSortIntWithArrayPair_Private(L,J,K,last-1);
232: PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));
233: return(0);
234: }
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.
242: Not Collective
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)
250: Level: intermediate
252: Concepts: sorting^ints with array pair
254: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
255: @*/
256: PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K)
257: {
259: PetscInt j,k,tmp,ik;
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: PetscSortIntWithArrayPair_Private(L,J,K,n-1);
274: }
275: return(0);
276: }
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: {
289: PetscMPIInt i,vl,last,tmp;
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: return(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: PetscSortMPIIntWithArray_Private(v,V,last-1);
306: PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
307: return(0);
308: }
312: /*@
313: PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
314: changes a second array to match the sorted first array.
316: Not Collective
318: Input Parameters:
319: + n - number of values
320: . i - array of integers
321: - I - second array of integers
323: Level: intermediate
325: Concepts: sorting^ints with array
327: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
328: @*/
329: PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[])
330: {
332: PetscMPIInt j,k,tmp,ik;
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: PetscSortMPIIntWithArray_Private(i,Ii,n-1);
347: }
348: return(0);
349: }
351: /* -----------------------------------------------------------------------*/
352: #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;}
356: /*
357: Modified from PetscSortIntWithArray_Private().
358: */
359: static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right)
360: {
362: PetscInt i,vl,last,tmp;
363: PetscScalar stmp;
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: return(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: PetscSortIntWithScalarArray_Private(v,V,last-1);
380: PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));
381: return(0);
382: }
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.
390: Not Collective
392: Input Parameters:
393: + n - number of values
394: . i - array of integers
395: - I - second array of scalars
397: Level: intermediate
399: Concepts: sorting^ints with array
401: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
402: @*/
403: PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[])
404: {
406: PetscInt j,k,tmp,ik;
407: PetscScalar stmp;
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: PetscSortIntWithScalarArray_Private(i,Ii,n-1);
422: }
423: return(0);
424: }
430: /*@
431: PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers.
432: The additional arrays are the same length as sorted arrays and are merged
433: in the order determined by the merging of the sorted pair.
435: Not Collective
437: Input Parameters:
438: + an - number of values in the first array
439: . aI - first sorted array of integers
440: . aJ - first additional array of integers
441: . bn - number of values in the second array
442: . bI - second array of integers
443: - bJ - second additional of integers
445: Output Parameters:
446: + n - number of values in the merged array (== an + bn)
447: . I - merged sorted array
448: - J - merged additional array
450: Level: intermediate
452: Concepts: merging^arrays
454: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
455: @*/
456: PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
457: {
459: PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k;
461: n_ = an + bn;
462: *n = n_;
463: if(!L_) {
464: PetscMalloc(n_*sizeof(PetscInt), L);
465: L_ = *L;
466: }
467: if(!J_){
468: PetscMalloc(n_*sizeof(PetscInt), &J_);
469: J_ = *J;
470: }
471: k = ak = bk = 0;
472: while(ak < an && bk < bn) {
473: if(aI[ak] <= bI[bk]) {
474: L_[k] = aI[ak];
475: J_[k] = aJ[ak];
476: ++ak;
477: ++k;
478: }
479: else {
480: L_[k] = bI[bk];
481: J_[k] = bJ[bk];
482: ++bk;
483: ++k;
484: }
485: }
486: if(ak < an) {
487: PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
488: PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));
489: k += (an-ak);
490: }
491: if(bk < bn) {
492: PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
493: PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));
494: k += (bn-bk);
495: }
496: return(0);
497: }
502: /*@
503: PetscProcessTree - Prepares tree data to be displayed graphically
505: Not Collective
507: Input Parameters:
508: + n - number of values
509: . mask - indicates those entries in the tree, location 0 is always masked
510: - parentid - indicates the parent of each entry
512: Output Parameters:
513: + Nlevels - the number of levels
514: . Level - for each node tells its level
515: . Levelcnts - the number of nodes on each level
516: . Idbylevel - a list of ids on each of the levels, first level followed by second etc
517: - Column - for each id tells its column index
519: Level: intermediate
522: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
523: @*/
524: PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
525: {
526: PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
528: PetscBool done = PETSC_FALSE;
531: if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
532: for (i=0; i<n; i++) {
533: if (mask[i]) continue;
534: if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
535: if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
536: }
539: for (i=0; i<n; i++) {
540: if (!mask[i]) nmask++;
541: }
543: /* determine the level in the tree of each node */
544: PetscMalloc(n*sizeof(PetscInt),&level);
545: PetscMemzero(level,n*sizeof(PetscInt));
546: level[0] = 1;
547: while (!done) {
548: done = PETSC_TRUE;
549: for (i=0; i<n; i++) {
550: if (mask[i]) continue;
551: if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
552: else if (!level[i]) done = PETSC_FALSE;
553: }
554: }
555: for (i=0; i<n; i++) {
556: level[i]--;
557: nlevels = PetscMax(nlevels,level[i]);
558: }
560: /* count the number of nodes on each level and its max */
561: PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);
562: PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));
563: for (i=0; i<n; i++) {
564: if (mask[i]) continue;
565: levelcnt[level[i]-1]++;
566: }
567: for (i=0; i<nlevels;i++) {
568: levelmax = PetscMax(levelmax,levelcnt[i]);
569: }
571: /* for each level sort the ids by the parent id */
572: PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);
573: PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);
574: for (j=1; j<=nlevels;j++) {
575: cnt = 0;
576: for (i=0; i<n; i++) {
577: if (mask[i]) continue;
578: if (level[i] != j) continue;
579: workid[cnt] = i;
580: workparentid[cnt++] = parentid[i];
581: }
582: /* PetscIntView(cnt,workparentid,0);
583: PetscIntView(cnt,workid,0);
584: PetscSortIntWithArray(cnt,workparentid,workid);
585: PetscIntView(cnt,workparentid,0);
586: PetscIntView(cnt,workid,0);*/
587: PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));
588: tcnt += cnt;
589: }
590: if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
591: PetscFree2(workid,workparentid);
593: /* for each node list its column */
594: PetscMalloc(n*sizeof(PetscInt),&column);
595: cnt = 0;
596: for (j=0; j<nlevels; j++) {
597: for (i=0; i<levelcnt[j]; i++) {
598: column[idbylevel[cnt++]] = i;
599: }
600: }
602: *Nlevels = nlevels;
603: *Level = level;
604: *Levelcnt = levelcnt;
605: *Idbylevel = idbylevel;
606: *Column = column;
607: return(0);
608: }