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 <petscsys.h>
8: #include <petsc/private/petscimpl.h>
10: #define SWAP(a, b, t) \
11: do { \
12: t = a; \
13: a = b; \
14: b = t; \
15: } while (0)
17: /*@
18: PetscSortedReal - Determines whether the array of `PetscReal` is sorted.
20: Not Collective
22: Input Parameters:
23: + n - number of values
24: - X - array of integers
26: Output Parameter:
27: . sorted - flag whether the array is sorted
29: Level: intermediate
31: .seealso: `PetscSortReal()`, `PetscSortedInt()`, `PetscSortedMPIInt()`
32: @*/
33: PetscErrorCode PetscSortedReal(PetscInt n, const PetscReal X[], PetscBool *sorted)
34: {
35: PetscFunctionBegin;
36: PetscSorted(n, X, *sorted);
37: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
47: if (right <= 1) {
48: if (right == 1) {
49: if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
50: }
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
53: SWAP(v[0], v[right / 2], tmp);
54: vl = v[0];
55: last = 0;
56: for (i = 1; i <= right; i++) {
57: if (v[i] < vl) {
58: last++;
59: SWAP(v[last], v[i], tmp);
60: }
61: }
62: SWAP(v[0], v[last], tmp);
63: PetscCall(PetscSortReal_Private(v, last - 1));
64: PetscCall(PetscSortReal_Private(v + last + 1, right - (last + 1)));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*@
69: PetscSortReal - Sorts an array of `PetscReal` in place in increasing order.
71: Not Collective
73: Input Parameters:
74: + n - number of values
75: - v - array of doubles
77: Level: intermediate
79: Note:
80: This function serves as an alternative to `PetscRealSortSemiOrdered()`, and may perform faster especially if the array
81: is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
82: code to see which routine is fastest.
84: .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
85: @*/
86: PetscErrorCode PetscSortReal(PetscInt n, PetscReal v[])
87: {
88: PetscInt j, k;
89: PetscReal tmp, vk;
91: PetscFunctionBegin;
92: PetscAssertPointer(v, 2);
93: if (n < 8) {
94: for (k = 0; k < n; k++) {
95: vk = v[k];
96: for (j = k + 1; j < n; j++) {
97: if (vk > v[j]) {
98: SWAP(v[k], v[j], tmp);
99: vk = v[k];
100: }
101: }
102: }
103: } else PetscCall(PetscSortReal_Private(v, n - 1));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: #define SWAP2ri(a, b, c, d, rt, it) \
108: do { \
109: rt = a; \
110: a = b; \
111: b = rt; \
112: it = c; \
113: c = d; \
114: d = it; \
115: } while (0)
117: /* modified from PetscSortIntWithArray_Private */
118: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscInt right)
119: {
120: PetscInt i, last, itmp;
121: PetscReal rvl, rtmp;
123: PetscFunctionBegin;
124: if (right <= 1) {
125: if (right == 1) {
126: if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
130: SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
131: rvl = v[0];
132: last = 0;
133: for (i = 1; i <= right; i++) {
134: if (v[i] < rvl) {
135: last++;
136: SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
137: }
138: }
139: SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
140: PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1));
141: PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1)));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
144: /*@
145: PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
146: changes a second `PetscInt` array to match the sorted first array.
148: Not Collective
150: Input Parameters:
151: + n - number of values
152: . Ii - array of integers
153: - r - second array of integers
155: Level: intermediate
157: .seealso: `PetscSortReal()`
158: @*/
159: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n, PetscReal r[], PetscInt Ii[])
160: {
161: PetscInt j, k, itmp;
162: PetscReal rk, rtmp;
164: PetscFunctionBegin;
165: PetscAssertPointer(r, 2);
166: PetscAssertPointer(Ii, 3);
167: if (n < 8) {
168: for (k = 0; k < n; k++) {
169: rk = r[k];
170: for (j = k + 1; j < n; j++) {
171: if (rk > r[j]) {
172: SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
173: rk = r[k];
174: }
175: }
176: }
177: } else {
178: PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1));
179: }
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s
186: Not Collective
188: Input Parameters:
189: + key - the value to locate
190: . n - number of values in the array
191: . t - array of values
192: - eps - tolerance used to compare
194: Output Parameter:
195: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
197: Level: intermediate
199: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
200: @*/
201: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
202: {
203: PetscInt lo = 0, hi = n;
205: PetscFunctionBegin;
206: PetscAssertPointer(loc, 5);
207: if (!n) {
208: *loc = -1;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
211: PetscAssertPointer(t, 3);
212: while (hi - lo > 1) {
213: PetscInt mid = lo + (hi - lo) / 2;
214: 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]);
215: if (key < t[mid]) hi = mid;
216: else lo = mid;
217: }
218: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
225: Not Collective
227: Input Parameters:
228: + n - number of values
229: - v - array of doubles
231: Output Parameter:
232: . n - number of non-redundant values
234: Level: intermediate
236: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
237: @*/
238: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
239: {
240: PetscInt i, s = 0, N = *n, b = 0;
242: PetscFunctionBegin;
243: PetscCall(PetscSortReal(N, v));
244: for (i = 0; i < N - 1; i++) {
245: if (v[b + s + 1] != v[b]) {
246: v[b + 1] = v[b + s + 1];
247: b++;
248: } else s++;
249: }
250: *n = N - s;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
257: Not Collective
259: Input Parameters:
260: + ncut - splitting index
261: - n - number of values to sort
263: Input/Output Parameters:
264: + a - array of values, on output the values are permuted such that its elements satisfy:
265: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
266: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
267: - idx - index for array a, on output permuted accordingly
269: Level: intermediate
271: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
272: @*/
273: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
274: {
275: PetscInt i, mid, last, itmp, j, first;
276: PetscScalar d, tmp;
277: PetscReal abskey;
279: PetscFunctionBegin;
280: first = 0;
281: last = n - 1;
282: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
284: while (1) {
285: mid = first;
286: d = a[mid];
287: abskey = PetscAbsScalar(d);
288: i = last;
289: for (j = first + 1; j <= i; ++j) {
290: d = a[j];
291: if (PetscAbsScalar(d) >= abskey) {
292: ++mid;
293: /* interchange */
294: tmp = a[mid];
295: itmp = idx[mid];
296: a[mid] = a[j];
297: idx[mid] = idx[j];
298: a[j] = tmp;
299: idx[j] = itmp;
300: }
301: }
303: /* interchange */
304: tmp = a[mid];
305: itmp = idx[mid];
306: a[mid] = a[first];
307: idx[mid] = idx[first];
308: a[first] = tmp;
309: idx[first] = itmp;
311: /* test for while loop */
312: if (mid == ncut) break;
313: else if (mid > ncut) last = mid - 1;
314: else first = mid + 1;
315: }
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
322: Not Collective
324: Input Parameters:
325: + ncut - splitting index
326: - n - number of values to sort
328: Input/Output Parameters:
329: + a - array of values, on output the values are permuted such that its elements satisfy:
330: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
331: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
332: - idx - index for array a, on output permuted accordingly
334: Level: intermediate
336: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
337: @*/
338: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
339: {
340: PetscInt i, mid, last, itmp, j, first;
341: PetscReal d, tmp;
342: PetscReal abskey;
344: PetscFunctionBegin;
345: first = 0;
346: last = n - 1;
347: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
349: while (1) {
350: mid = first;
351: d = a[mid];
352: abskey = PetscAbsReal(d);
353: i = last;
354: for (j = first + 1; j <= i; ++j) {
355: d = a[j];
356: if (PetscAbsReal(d) >= abskey) {
357: ++mid;
358: /* interchange */
359: tmp = a[mid];
360: itmp = idx[mid];
361: a[mid] = a[j];
362: idx[mid] = idx[j];
363: a[j] = tmp;
364: idx[j] = itmp;
365: }
366: }
368: /* interchange */
369: tmp = a[mid];
370: itmp = idx[mid];
371: a[mid] = a[first];
372: idx[mid] = idx[first];
373: a[first] = tmp;
374: idx[first] = itmp;
376: /* test for while loop */
377: if (mid == ncut) break;
378: else if (mid > ncut) last = mid - 1;
379: else first = mid + 1;
380: }
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }