Actual source code: sortd.c
2: /*
3: This file contains routines for sorting doubles. Values are sorted in place.
4: These are provided because the general sort routines incur a great deal
5: of overhead in calling the comparision routines.
7: */
8: #include <petscsys.h>
9: #include <petsc/private/petscimpl.h>
11: #define SWAP(a,b,t) {t=a;a=b;b=t;}
13: /*@
14: PetscSortedReal - Determines whether the array is sorted.
16: Not Collective
18: Input Parameters:
19: + n - number of values
20: - X - array of integers
22: Output Parameters:
23: . sorted - flag whether the array is sorted
25: Level: intermediate
27: .seealso: PetscSortReal(), PetscSortedInt(), PetscSortedMPIInt()
28: @*/
29: PetscErrorCode PetscSortedReal(PetscInt n,const PetscReal X[],PetscBool *sorted)
30: {
32: PetscSorted(n,X,*sorted);
33: return(0);
34: }
36: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
37: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
38: {
39: PetscInt i,last;
40: PetscReal vl,tmp;
43: if (right <= 1) {
44: if (right == 1) {
45: if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
46: }
47: return(0);
48: }
49: SWAP(v[0],v[right/2],tmp);
50: vl = v[0];
51: last = 0;
52: for (i=1; i<=right; i++) {
53: if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
54: }
55: SWAP(v[0],v[last],tmp);
56: PetscSortReal_Private(v,last-1);
57: PetscSortReal_Private(v+last+1,right-(last+1));
58: return(0);
59: }
61: /*@
62: PetscSortReal - Sorts an array of doubles in place in increasing order.
64: Not Collective
66: Input Parameters:
67: + n - number of values
68: - v - array of doubles
70: Notes:
71: This function serves as an alternative to PetscRealSortSemiOrdered(), and may perform faster especially if the array
72: is completely random. There are exceptions to this and so it is __highly__ recomended that the user benchmark their
73: code to see which routine is fastest.
75: Level: intermediate
77: .seealso: PetscRealSortSemiOrdered(), PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
78: @*/
79: PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
80: {
81: PetscInt j,k;
82: PetscReal tmp,vk;
86: if (n<8) {
87: for (k=0; k<n; k++) {
88: vk = v[k];
89: for (j=k+1; j<n; j++) {
90: if (vk > v[j]) {
91: SWAP(v[k],v[j],tmp);
92: vk = v[k];
93: }
94: }
95: }
96: } else PetscSortReal_Private(v,n-1);
97: return(0);
98: }
100: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
102: /* modified from PetscSortIntWithArray_Private */
103: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
104: {
106: PetscInt i,last,itmp;
107: PetscReal rvl,rtmp;
110: if (right <= 1) {
111: if (right == 1) {
112: if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
113: }
114: return(0);
115: }
116: SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
117: rvl = v[0];
118: last = 0;
119: for (i=1; i<=right; i++) {
120: if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
121: }
122: SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
123: PetscSortRealWithArrayInt_Private(v,V,last-1);
124: PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
125: return(0);
126: }
127: /*@
128: PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
129: changes a second PetscInt array to match the sorted first array.
131: Not Collective
133: Input Parameters:
134: + n - number of values
135: . i - array of integers
136: - I - second array of integers
138: Level: intermediate
140: .seealso: PetscSortReal()
141: @*/
142: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
143: {
145: PetscInt j,k,itmp;
146: PetscReal rk,rtmp;
151: if (n<8) {
152: for (k=0; k<n; k++) {
153: rk = r[k];
154: for (j=k+1; j<n; j++) {
155: if (rk > r[j]) {
156: SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
157: rk = r[k];
158: }
159: }
160: }
161: } else {
162: PetscSortRealWithArrayInt_Private(r,Ii,n-1);
163: }
164: return(0);
165: }
167: /*@
168: PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
170: Not Collective
172: Input Parameters:
173: + key - the value to locate
174: . n - number of values in the array
175: . ii - array of values
176: - eps - tolerance used to compare
178: Output Parameter:
179: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
181: Level: intermediate
183: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
184: @*/
185: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
186: {
187: PetscInt lo = 0,hi = n;
191: if (!n) {*loc = -1; return(0);}
194: while (hi - lo > 1) {
195: PetscInt mid = lo + (hi - lo)/2;
196: if (key < t[mid]) hi = mid;
197: else lo = mid;
198: }
199: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
200: return(0);
201: }
203: /*@
204: PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
206: Not Collective
208: Input Parameters:
209: + n - number of values
210: - v - array of doubles
212: Output Parameter:
213: . n - number of non-redundant values
215: Level: intermediate
217: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
218: @*/
219: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
220: {
222: PetscInt i,s = 0,N = *n, b = 0;
225: PetscSortReal(N,v);
226: for (i=0; i<N-1; i++) {
227: if (v[b+s+1] != v[b]) {
228: v[b+1] = v[b+s+1]; b++;
229: } else s++;
230: }
231: *n = N - s;
232: return(0);
233: }
235: /*@
236: PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
238: Not Collective
240: Input Parameters:
241: + ncut - splitig index
242: . n - number of values to sort
243: . a - array of values
244: - idx - index for array a
246: Output Parameters:
247: + a - permuted array of values such that its elements satisfy:
248: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
249: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
250: - idx - permuted index of array a
252: Level: intermediate
254: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
255: @*/
256: PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
257: {
258: PetscInt i,mid,last,itmp,j,first;
259: PetscScalar d,tmp;
260: PetscReal abskey;
263: first = 0;
264: last = n-1;
265: if (ncut < first || ncut > last) return(0);
267: while (1) {
268: mid = first;
269: d = a[mid];
270: abskey = PetscAbsScalar(d);
271: i = last;
272: for (j = first + 1; j <= i; ++j) {
273: d = a[j];
274: if (PetscAbsScalar(d) >= abskey) {
275: ++mid;
276: /* interchange */
277: tmp = a[mid]; itmp = idx[mid];
278: a[mid] = a[j]; idx[mid] = idx[j];
279: a[j] = tmp; idx[j] = itmp;
280: }
281: }
283: /* interchange */
284: tmp = a[mid]; itmp = idx[mid];
285: a[mid] = a[first]; idx[mid] = idx[first];
286: a[first] = tmp; idx[first] = itmp;
288: /* test for while loop */
289: if (mid == ncut) break;
290: else if (mid > ncut) last = mid - 1;
291: else first = mid + 1;
292: }
293: return(0);
294: }
296: /*@
297: PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
299: Not Collective
301: Input Parameters:
302: + ncut - splitig index
303: . n - number of values to sort
304: . a - array of values in PetscReal
305: - idx - index for array a
307: Output Parameters:
308: + a - permuted array of real values such that its elements satisfy:
309: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
310: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
311: - idx - permuted index of array a
313: Level: intermediate
315: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
316: @*/
317: PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
318: {
319: PetscInt i,mid,last,itmp,j,first;
320: PetscReal d,tmp;
321: PetscReal abskey;
324: first = 0;
325: last = n-1;
326: if (ncut < first || ncut > last) return(0);
328: while (1) {
329: mid = first;
330: d = a[mid];
331: abskey = PetscAbsReal(d);
332: i = last;
333: for (j = first + 1; j <= i; ++j) {
334: d = a[j];
335: if (PetscAbsReal(d) >= abskey) {
336: ++mid;
337: /* interchange */
338: tmp = a[mid]; itmp = idx[mid];
339: a[mid] = a[j]; idx[mid] = idx[j];
340: a[j] = tmp; idx[j] = itmp;
341: }
342: }
344: /* interchange */
345: tmp = a[mid]; itmp = idx[mid];
346: a[mid] = a[first]; idx[mid] = idx[first];
347: a[first] = tmp; idx[first] = itmp;
349: /* test for while loop */
350: if (mid == ncut) break;
351: else if (mid > ncut) last = mid - 1;
352: else first = mid + 1;
353: }
354: return(0);
355: }