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 comparison 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: {
31: PetscSorted(n,X,*sorted);
32: return 0;
33: }
35: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
36: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
37: {
38: PetscInt i,last;
39: PetscReal vl,tmp;
41: if (right <= 1) {
42: if (right == 1) {
43: if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
44: }
45: return 0;
46: }
47: SWAP(v[0],v[right/2],tmp);
48: vl = v[0];
49: last = 0;
50: for (i=1; i<=right; i++) {
51: if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
52: }
53: SWAP(v[0],v[last],tmp);
54: PetscSortReal_Private(v,last-1);
55: PetscSortReal_Private(v+last+1,right-(last+1));
56: return 0;
57: }
59: /*@
60: PetscSortReal - Sorts an array of doubles in place in increasing order.
62: Not Collective
64: Input Parameters:
65: + n - number of values
66: - v - array of doubles
68: Notes:
69: This function serves as an alternative to PetscRealSortSemiOrdered(), and may perform faster especially if the array
70: is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
71: code to see which routine is fastest.
73: Level: intermediate
75: .seealso: PetscRealSortSemiOrdered(), PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
76: @*/
77: PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
78: {
79: PetscInt j,k;
80: PetscReal tmp,vk;
83: if (n<8) {
84: for (k=0; k<n; k++) {
85: vk = v[k];
86: for (j=k+1; j<n; j++) {
87: if (vk > v[j]) {
88: SWAP(v[k],v[j],tmp);
89: vk = v[k];
90: }
91: }
92: }
93: } else PetscSortReal_Private(v,n-1);
94: return 0;
95: }
97: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
99: /* modified from PetscSortIntWithArray_Private */
100: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
101: {
102: PetscInt i,last,itmp;
103: PetscReal rvl,rtmp;
105: if (right <= 1) {
106: if (right == 1) {
107: if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
108: }
109: return 0;
110: }
111: SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
112: rvl = v[0];
113: last = 0;
114: for (i=1; i<=right; i++) {
115: if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
116: }
117: SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
118: PetscSortRealWithArrayInt_Private(v,V,last-1);
119: PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
120: return 0;
121: }
122: /*@
123: PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
124: changes a second PetscInt array to match the sorted first array.
126: Not Collective
128: Input Parameters:
129: + n - number of values
130: . i - array of integers
131: - I - second array of integers
133: Level: intermediate
135: .seealso: PetscSortReal()
136: @*/
137: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
138: {
139: PetscInt j,k,itmp;
140: PetscReal rk,rtmp;
144: if (n<8) {
145: for (k=0; k<n; k++) {
146: rk = r[k];
147: for (j=k+1; j<n; j++) {
148: if (rk > r[j]) {
149: SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
150: rk = r[k];
151: }
152: }
153: }
154: } else {
155: PetscSortRealWithArrayInt_Private(r,Ii,n-1);
156: }
157: return 0;
158: }
160: /*@
161: PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
163: Not Collective
165: Input Parameters:
166: + key - the value to locate
167: . n - number of values in the array
168: . ii - array of values
169: - eps - tolerance used to compare
171: Output Parameter:
172: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
174: Level: intermediate
176: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
177: @*/
178: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
179: {
180: PetscInt lo = 0,hi = n;
183: if (!n) {*loc = -1; return 0;}
186: while (hi - lo > 1) {
187: PetscInt mid = lo + (hi - lo)/2;
188: if (key < t[mid]) hi = mid;
189: else lo = mid;
190: }
191: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
192: return 0;
193: }
195: /*@
196: PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
198: Not Collective
200: Input Parameters:
201: + n - number of values
202: - v - array of doubles
204: Output Parameter:
205: . n - number of non-redundant values
207: Level: intermediate
209: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
210: @*/
211: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
212: {
213: PetscInt i,s = 0,N = *n, b = 0;
215: PetscSortReal(N,v);
216: for (i=0; i<N-1; i++) {
217: if (v[b+s+1] != v[b]) {
218: v[b+1] = v[b+s+1]; b++;
219: } else s++;
220: }
221: *n = N - s;
222: return 0;
223: }
225: /*@
226: PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
228: Not Collective
230: Input Parameters:
231: + ncut - splitig index
232: - n - number of values to sort
234: Input/Output Parameters:
235: + a - array of values, on output the values are permuted such that its elements satisfy:
236: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
237: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
238: - idx - index for array a, on output permuted accordingly
240: Level: intermediate
242: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
243: @*/
244: PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
245: {
246: PetscInt i,mid,last,itmp,j,first;
247: PetscScalar d,tmp;
248: PetscReal abskey;
250: first = 0;
251: last = n-1;
252: if (ncut < first || ncut > last) return 0;
254: while (1) {
255: mid = first;
256: d = a[mid];
257: abskey = PetscAbsScalar(d);
258: i = last;
259: for (j = first + 1; j <= i; ++j) {
260: d = a[j];
261: if (PetscAbsScalar(d) >= abskey) {
262: ++mid;
263: /* interchange */
264: tmp = a[mid]; itmp = idx[mid];
265: a[mid] = a[j]; idx[mid] = idx[j];
266: a[j] = tmp; idx[j] = itmp;
267: }
268: }
270: /* interchange */
271: tmp = a[mid]; itmp = idx[mid];
272: a[mid] = a[first]; idx[mid] = idx[first];
273: a[first] = tmp; idx[first] = itmp;
275: /* test for while loop */
276: if (mid == ncut) break;
277: else if (mid > ncut) last = mid - 1;
278: else first = mid + 1;
279: }
280: return 0;
281: }
283: /*@
284: PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
286: Not Collective
288: Input Parameters:
289: + ncut - splitig index
290: - n - number of values to sort
292: Input/Output Parameters:
293: + a - array of values, on output the values are permuted such that its elements satisfy:
294: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
295: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
296: - idx - index for array a, on output permuted accordingly
298: Level: intermediate
300: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
301: @*/
302: PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
303: {
304: PetscInt i,mid,last,itmp,j,first;
305: PetscReal d,tmp;
306: PetscReal abskey;
308: first = 0;
309: last = n-1;
310: if (ncut < first || ncut > last) return 0;
312: while (1) {
313: mid = first;
314: d = a[mid];
315: abskey = PetscAbsReal(d);
316: i = last;
317: for (j = first + 1; j <= i; ++j) {
318: d = a[j];
319: if (PetscAbsReal(d) >= abskey) {
320: ++mid;
321: /* interchange */
322: tmp = a[mid]; itmp = idx[mid];
323: a[mid] = a[j]; idx[mid] = idx[j];
324: a[j] = tmp; idx[j] = itmp;
325: }
326: }
328: /* interchange */
329: tmp = a[mid]; itmp = idx[mid];
330: a[mid] = a[first]; idx[mid] = idx[first];
331: a[first] = tmp; idx[first] = itmp;
333: /* test for while loop */
334: if (mid == ncut) break;
335: else if (mid > ncut) last = mid - 1;
336: else first = mid + 1;
337: }
338: return 0;
339: }