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) \
12: { \
13: t = a; \
14: a = b; \
15: b = t; \
16: }
18: /*@
19: PetscSortedReal - Determines whether the array of `PetscReal` is sorted.
21: Not Collective
23: Input Parameters:
24: + n - number of values
25: - X - array of integers
27: Output Parameters:
28: . sorted - flag whether the array is sorted
30: Level: intermediate
32: .seealso: `PetscSortReal()`, `PetscSortedInt()`, `PetscSortedMPIInt()`
33: @*/
34: PetscErrorCode PetscSortedReal(PetscInt n, const PetscReal X[], PetscBool *sorted)
35: {
36: PetscSorted(n, X, *sorted);
37: return 0;
38: }
40: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
41: static PetscErrorCode PetscSortReal_Private(PetscReal *v, PetscInt right)
42: {
43: PetscInt i, last;
44: PetscReal vl, tmp;
46: if (right <= 1) {
47: if (right == 1) {
48: if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
49: }
50: return 0;
51: }
52: SWAP(v[0], v[right / 2], tmp);
53: vl = v[0];
54: last = 0;
55: for (i = 1; i <= right; i++) {
56: if (v[i] < vl) {
57: last++;
58: SWAP(v[last], v[i], tmp);
59: }
60: }
61: SWAP(v[0], v[last], tmp);
62: PetscSortReal_Private(v, last - 1);
63: PetscSortReal_Private(v + last + 1, right - (last + 1));
64: return 0;
65: }
67: /*@
68: PetscSortReal - Sorts an array of `PetscReal` in place in increasing order.
70: Not Collective
72: Input Parameters:
73: + n - number of values
74: - v - array of doubles
76: Note:
77: This function serves as an alternative to `PetscRealSortSemiOrdered()`, and may perform faster especially if the array
78: is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
79: code to see which routine is fastest.
81: Level: intermediate
83: .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
84: @*/
85: PetscErrorCode PetscSortReal(PetscInt n, PetscReal v[])
86: {
87: PetscInt j, k;
88: PetscReal tmp, vk;
91: if (n < 8) {
92: for (k = 0; k < n; k++) {
93: vk = v[k];
94: for (j = k + 1; j < n; j++) {
95: if (vk > v[j]) {
96: SWAP(v[k], v[j], tmp);
97: vk = v[k];
98: }
99: }
100: }
101: } else PetscSortReal_Private(v, n - 1);
102: return 0;
103: }
105: #define SWAP2ri(a, b, c, d, rt, it) \
106: { \
107: rt = a; \
108: a = b; \
109: b = rt; \
110: it = c; \
111: c = d; \
112: d = it; \
113: }
115: /* modified from PetscSortIntWithArray_Private */
116: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscInt right)
117: {
118: PetscInt i, last, itmp;
119: PetscReal rvl, rtmp;
121: if (right <= 1) {
122: if (right == 1) {
123: if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
124: }
125: return 0;
126: }
127: SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
128: rvl = v[0];
129: last = 0;
130: for (i = 1; i <= right; i++) {
131: if (v[i] < rvl) {
132: last++;
133: SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
134: }
135: }
136: SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
137: PetscSortRealWithArrayInt_Private(v, V, last - 1);
138: PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1));
139: return 0;
140: }
141: /*@
142: PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
143: changes a second `PetscInt` array to match the sorted first array.
145: Not Collective
147: Input Parameters:
148: + n - number of values
149: . i - array of integers
150: - I - second array of integers
152: Level: intermediate
154: .seealso: `PetscSortReal()`
155: @*/
156: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n, PetscReal r[], PetscInt Ii[])
157: {
158: PetscInt j, k, itmp;
159: PetscReal rk, rtmp;
163: if (n < 8) {
164: for (k = 0; k < n; k++) {
165: rk = r[k];
166: for (j = k + 1; j < n; j++) {
167: if (rk > r[j]) {
168: SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
169: rk = r[k];
170: }
171: }
172: }
173: } else {
174: PetscSortRealWithArrayInt_Private(r, Ii, n - 1);
175: }
176: return 0;
177: }
179: /*@
180: PetscFindReal - Finds a PetscReal` in a sorted array of `PetscReal`s
182: Not Collective
184: Input Parameters:
185: + key - the value to locate
186: . n - number of values in the array
187: . ii - array of values
188: - eps - tolerance used to compare
190: Output Parameter:
191: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
193: Level: intermediate
195: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
196: @*/
197: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
198: {
199: PetscInt lo = 0, hi = n;
202: if (!n) {
203: *loc = -1;
204: return 0;
205: }
208: while (hi - lo > 1) {
209: PetscInt mid = lo + (hi - lo) / 2;
210: if (key < t[mid]) hi = mid;
211: else lo = mid;
212: }
213: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
214: return 0;
215: }
217: /*@
218: PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
220: Not Collective
222: Input Parameters:
223: + n - number of values
224: - v - array of doubles
226: Output Parameter:
227: . n - number of non-redundant values
229: Level: intermediate
231: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
232: @*/
233: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
234: {
235: PetscInt i, s = 0, N = *n, b = 0;
237: PetscSortReal(N, v);
238: for (i = 0; i < N - 1; i++) {
239: if (v[b + s + 1] != v[b]) {
240: v[b + 1] = v[b + s + 1];
241: b++;
242: } else s++;
243: }
244: *n = N - s;
245: return 0;
246: }
248: /*@
249: PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
251: Not Collective
253: Input Parameters:
254: + ncut - splitig index
255: - n - number of values to sort
257: Input/Output Parameters:
258: + a - array of values, on output the values are permuted such that its elements satisfy:
259: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
260: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
261: - idx - index for array a, on output permuted accordingly
263: Level: intermediate
265: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
266: @*/
267: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
268: {
269: PetscInt i, mid, last, itmp, j, first;
270: PetscScalar d, tmp;
271: PetscReal abskey;
273: first = 0;
274: last = n - 1;
275: if (ncut < first || ncut > last) return 0;
277: while (1) {
278: mid = first;
279: d = a[mid];
280: abskey = PetscAbsScalar(d);
281: i = last;
282: for (j = first + 1; j <= i; ++j) {
283: d = a[j];
284: if (PetscAbsScalar(d) >= abskey) {
285: ++mid;
286: /* interchange */
287: tmp = a[mid];
288: itmp = idx[mid];
289: a[mid] = a[j];
290: idx[mid] = idx[j];
291: a[j] = tmp;
292: idx[j] = itmp;
293: }
294: }
296: /* interchange */
297: tmp = a[mid];
298: itmp = idx[mid];
299: a[mid] = a[first];
300: idx[mid] = idx[first];
301: a[first] = tmp;
302: idx[first] = itmp;
304: /* test for while loop */
305: if (mid == ncut) break;
306: else if (mid > ncut) last = mid - 1;
307: else first = mid + 1;
308: }
309: return 0;
310: }
312: /*@
313: PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
315: Not Collective
317: Input Parameters:
318: + ncut - splitig index
319: - n - number of values to sort
321: Input/Output Parameters:
322: + a - array of values, on output the values are permuted such that its elements satisfy:
323: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
324: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
325: - idx - index for array a, on output permuted accordingly
327: Level: intermediate
329: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
330: @*/
331: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
332: {
333: PetscInt i, mid, last, itmp, j, first;
334: PetscReal d, tmp;
335: PetscReal abskey;
337: first = 0;
338: last = n - 1;
339: if (ncut < first || ncut > last) return 0;
341: while (1) {
342: mid = first;
343: d = a[mid];
344: abskey = PetscAbsReal(d);
345: i = last;
346: for (j = first + 1; j <= i; ++j) {
347: d = a[j];
348: if (PetscAbsReal(d) >= abskey) {
349: ++mid;
350: /* interchange */
351: tmp = a[mid];
352: itmp = idx[mid];
353: a[mid] = a[j];
354: idx[mid] = idx[j];
355: a[j] = tmp;
356: idx[j] = itmp;
357: }
358: }
360: /* interchange */
361: tmp = a[mid];
362: itmp = idx[mid];
363: a[mid] = a[first];
364: idx[mid] = idx[first];
365: a[first] = tmp;
366: idx[first] = itmp;
368: /* test for while loop */
369: if (mid == ncut) break;
370: else if (mid > ncut) last = mid - 1;
371: else first = mid + 1;
372: }
373: return 0;
374: }