Actual source code: sortd.c
1: /*
2: This file contains routines for sorting doubles. Values are sorted in place.
3: These are provided because the general sort routines incur a great deal
4: of overhead in calling the comparison routines.
6: */
7: #include <petsc/private/petscimpl.h>
9: #define SWAP(a, b, t) \
10: do { \
11: t = a; \
12: a = b; \
13: b = t; \
14: } while (0)
16: /*@
17: PetscSortedReal - Determines whether the array of `PetscReal` is sorted.
19: Not Collective
21: Input Parameters:
22: + n - number of values
23: - X - array of integers
25: Output Parameter:
26: . sorted - flag whether the array is sorted
28: Level: intermediate
30: .seealso: `PetscSortReal()`, `PetscSortedInt()`, `PetscSortedMPIInt()`
31: @*/
32: PetscErrorCode PetscSortedReal(PetscCount n, const PetscReal X[], PetscBool *sorted)
33: {
34: PetscFunctionBegin;
35: PetscSorted(n, X, *sorted);
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
40: static PetscErrorCode PetscSortReal_Private(PetscReal *v, PetscCount right)
41: {
42: PetscCount i, last;
43: PetscReal vl, tmp;
45: PetscFunctionBegin;
46: if (right <= 1) {
47: if (right == 1) {
48: if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscCall(PetscSortReal_Private(v, last - 1));
63: PetscCall(PetscSortReal_Private(v + last + 1, right - (last + 1)));
64: PetscFunctionReturn(PETSC_SUCCESS);
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: Level: intermediate
78: Note:
79: This function serves as an alternative to `PetscRealSortSemiOrdered()`, and may perform faster especially if the array
80: is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
81: code to see which routine is fastest.
83: .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
84: @*/
85: PetscErrorCode PetscSortReal(PetscCount n, PetscReal v[])
86: {
87: PetscFunctionBegin;
88: PetscAssertPointer(v, 2);
89: if (n < 8) {
90: PetscReal tmp, vk;
91: for (PetscCount k = 0; k < n; k++) {
92: vk = v[k];
93: for (PetscCount j = k + 1; j < n; j++) {
94: if (vk > v[j]) {
95: SWAP(v[k], v[j], tmp);
96: vk = v[k];
97: }
98: }
99: }
100: } else {
101: PetscInt N;
103: PetscCall(PetscIntCast(n, &N));
104: PetscCall(PetscSortReal_Private(v, N - 1));
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: #define SWAP2ri(a, b, c, d, rt, it) \
110: do { \
111: rt = a; \
112: a = b; \
113: b = rt; \
114: it = c; \
115: c = d; \
116: d = it; \
117: } while (0)
119: /* modified from PetscSortIntWithArray_Private */
120: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscCount right)
121: {
122: PetscCount i, last;
123: PetscInt itmp;
124: PetscReal rvl, rtmp;
126: PetscFunctionBegin;
127: if (right <= 1) {
128: if (right == 1) {
129: if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
130: }
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
133: SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
134: rvl = v[0];
135: last = 0;
136: for (i = 1; i <= right; i++) {
137: if (v[i] < rvl) {
138: last++;
139: SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
140: }
141: }
142: SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
143: PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1));
144: PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1)));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
147: /*@
148: PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
149: changes a second `PetscInt` array to match the sorted first array.
151: Not Collective
153: Input Parameters:
154: + n - number of values
155: . Ii - array of integers
156: - r - second array of integers
158: Level: intermediate
160: .seealso: `PetscSortReal()`
161: @*/
162: PetscErrorCode PetscSortRealWithArrayInt(PetscCount n, PetscReal r[], PetscInt Ii[])
163: {
164: PetscCount j, k;
165: PetscInt itmp;
166: PetscReal rk, rtmp;
168: PetscFunctionBegin;
169: PetscAssertPointer(r, 2);
170: PetscAssertPointer(Ii, 3);
171: if (n < 8) {
172: for (k = 0; k < n; k++) {
173: rk = r[k];
174: for (j = k + 1; j < n; j++) {
175: if (rk > r[j]) {
176: SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
177: rk = r[k];
178: }
179: }
180: }
181: } else {
182: PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1));
183: }
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*@
188: PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s
190: Not Collective
192: Input Parameters:
193: + key - the value to locate
194: . n - number of values in the array
195: . t - array of values
196: - eps - tolerance used to compare
198: Output Parameter:
199: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
201: Level: intermediate
203: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
204: @*/
205: PetscErrorCode PetscFindReal(PetscReal key, PetscCount n, const PetscReal t[], PetscReal eps, PetscInt *loc)
206: {
207: PetscInt lo = 0, hi;
209: PetscFunctionBegin;
210: PetscAssertPointer(loc, 5);
211: if (!n) {
212: *loc = -1;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
215: PetscAssertPointer(t, 3);
216: PetscCall(PetscIntCast(n, &hi));
217: while (hi - lo > 1) {
218: PetscInt mid = lo + (hi - lo) / 2;
219: PetscAssert(t[lo] <= t[mid] && t[mid] <= t[hi - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input array was not sorted: (%g, %g, %g)", (double)t[lo], (double)t[mid], (double)t[hi - 1]);
220: if (key < t[mid]) hi = mid;
221: else lo = mid;
222: }
223: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: /*@
228: PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
230: Not Collective
232: Input Parameters:
233: + n - initial number of values
234: - v - array of values
236: Note:
237: On output both `n` and `v` are modified with non-redundant values.
239: Level: intermediate
241: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
242: @*/
243: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
244: {
245: PetscInt i, s = 0, N = *n, b = 0;
247: PetscFunctionBegin;
248: PetscCall(PetscSortReal(N, v));
249: for (i = 0; i < N - 1; i++) {
250: if (v[b + s + 1] != v[b]) {
251: v[b + 1] = v[b + s + 1];
252: b++;
253: } else s++;
254: }
255: *n = N - s;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@
260: PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
262: Not Collective
264: Input Parameters:
265: + ncut - splitting index
266: - n - number of values to sort
268: Input/Output Parameters:
269: + a - array of values, on output the values are permuted such that its elements satisfy:
270: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
271: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
272: - idx - index for array a, on output permuted accordingly
274: Level: intermediate
276: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
277: @*/
278: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
279: {
280: PetscInt i, mid, last, itmp, j, first;
281: PetscScalar d, tmp;
282: PetscReal abskey;
284: PetscFunctionBegin;
285: first = 0;
286: last = n - 1;
287: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
289: while (1) {
290: mid = first;
291: d = a[mid];
292: abskey = PetscAbsScalar(d);
293: i = last;
294: for (j = first + 1; j <= i; ++j) {
295: d = a[j];
296: if (PetscAbsScalar(d) >= abskey) {
297: ++mid;
298: /* interchange */
299: tmp = a[mid];
300: itmp = idx[mid];
301: a[mid] = a[j];
302: idx[mid] = idx[j];
303: a[j] = tmp;
304: idx[j] = itmp;
305: }
306: }
308: /* interchange */
309: tmp = a[mid];
310: itmp = idx[mid];
311: a[mid] = a[first];
312: idx[mid] = idx[first];
313: a[first] = tmp;
314: idx[first] = itmp;
316: /* test for while loop */
317: if (mid == ncut) break;
318: else if (mid > ncut) last = mid - 1;
319: else first = mid + 1;
320: }
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
327: Not Collective
329: Input Parameters:
330: + ncut - splitting index
331: - n - number of values to sort
333: Input/Output Parameters:
334: + a - array of values, on output the values are permuted such that its elements satisfy:
335: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
336: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
337: - idx - index for array a, on output permuted accordingly
339: Level: intermediate
341: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
342: @*/
343: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
344: {
345: PetscInt i, mid, last, itmp, j, first;
346: PetscReal d, tmp;
347: PetscReal abskey;
349: PetscFunctionBegin;
350: first = 0;
351: last = n - 1;
352: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
354: while (1) {
355: mid = first;
356: d = a[mid];
357: abskey = PetscAbsReal(d);
358: i = last;
359: for (j = first + 1; j <= i; ++j) {
360: d = a[j];
361: if (PetscAbsReal(d) >= abskey) {
362: ++mid;
363: /* interchange */
364: tmp = a[mid];
365: itmp = idx[mid];
366: a[mid] = a[j];
367: idx[mid] = idx[j];
368: a[j] = tmp;
369: idx[j] = itmp;
370: }
371: }
373: /* interchange */
374: tmp = a[mid];
375: itmp = idx[mid];
376: a[mid] = a[first];
377: idx[mid] = idx[first];
378: a[first] = tmp;
379: idx[first] = itmp;
381: /* test for while loop */
382: if (mid == ncut) break;
383: else if (mid > ncut) last = mid - 1;
384: else first = mid + 1;
385: }
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }