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(PetscCount 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, PetscCount right)
42: {
43: PetscCount 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(PetscCount n, PetscReal v[])
87: {
88: PetscFunctionBegin;
89: PetscAssertPointer(v, 2);
90: if (n < 8) {
91: PetscReal tmp, vk;
92: for (PetscCount k = 0; k < n; k++) {
93: vk = v[k];
94: for (PetscCount 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 {
102: PetscInt N;
104: PetscCall(PetscIntCast(n, &N));
105: PetscCall(PetscSortReal_Private(v, N - 1));
106: }
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: #define SWAP2ri(a, b, c, d, rt, it) \
111: do { \
112: rt = a; \
113: a = b; \
114: b = rt; \
115: it = c; \
116: c = d; \
117: d = it; \
118: } while (0)
120: /* modified from PetscSortIntWithArray_Private */
121: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscCount right)
122: {
123: PetscCount i, last;
124: PetscInt itmp;
125: PetscReal rvl, rtmp;
127: PetscFunctionBegin;
128: if (right <= 1) {
129: if (right == 1) {
130: if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
131: }
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
134: SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
135: rvl = v[0];
136: last = 0;
137: for (i = 1; i <= right; i++) {
138: if (v[i] < rvl) {
139: last++;
140: SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
141: }
142: }
143: SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
144: PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1));
145: PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1)));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
148: /*@
149: PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
150: changes a second `PetscInt` array to match the sorted first array.
152: Not Collective
154: Input Parameters:
155: + n - number of values
156: . Ii - array of integers
157: - r - second array of integers
159: Level: intermediate
161: .seealso: `PetscSortReal()`
162: @*/
163: PetscErrorCode PetscSortRealWithArrayInt(PetscCount n, PetscReal r[], PetscInt Ii[])
164: {
165: PetscCount j, k;
166: PetscInt itmp;
167: PetscReal rk, rtmp;
169: PetscFunctionBegin;
170: PetscAssertPointer(r, 2);
171: PetscAssertPointer(Ii, 3);
172: if (n < 8) {
173: for (k = 0; k < n; k++) {
174: rk = r[k];
175: for (j = k + 1; j < n; j++) {
176: if (rk > r[j]) {
177: SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
178: rk = r[k];
179: }
180: }
181: }
182: } else {
183: PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1));
184: }
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@
189: PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s
191: Not Collective
193: Input Parameters:
194: + key - the value to locate
195: . n - number of values in the array
196: . t - array of values
197: - eps - tolerance used to compare
199: Output Parameter:
200: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
202: Level: intermediate
204: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
205: @*/
206: PetscErrorCode PetscFindReal(PetscReal key, PetscCount n, const PetscReal t[], PetscReal eps, PetscInt *loc)
207: {
208: PetscInt lo = 0, hi;
210: PetscFunctionBegin;
211: PetscAssertPointer(loc, 5);
212: if (!n) {
213: *loc = -1;
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
216: PetscAssertPointer(t, 3);
217: PetscCall(PetscIntCast(n, &hi));
218: while (hi - lo > 1) {
219: PetscInt mid = lo + (hi - lo) / 2;
220: 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]);
221: if (key < t[mid]) hi = mid;
222: else lo = mid;
223: }
224: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
231: Not Collective
233: Input Parameters:
234: + n - number of values
235: - v - array of doubles
237: Output Parameter:
238: . n - number of non-redundant values
240: Level: intermediate
242: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
243: @*/
244: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
245: {
246: PetscInt i, s = 0, N = *n, b = 0;
248: PetscFunctionBegin;
249: PetscCall(PetscSortReal(N, v));
250: for (i = 0; i < N - 1; i++) {
251: if (v[b + s + 1] != v[b]) {
252: v[b + 1] = v[b + s + 1];
253: b++;
254: } else s++;
255: }
256: *n = N - s;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
263: Not Collective
265: Input Parameters:
266: + ncut - splitting index
267: - n - number of values to sort
269: Input/Output Parameters:
270: + a - array of values, on output the values are permuted such that its elements satisfy:
271: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
272: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
273: - idx - index for array a, on output permuted accordingly
275: Level: intermediate
277: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
278: @*/
279: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
280: {
281: PetscInt i, mid, last, itmp, j, first;
282: PetscScalar d, tmp;
283: PetscReal abskey;
285: PetscFunctionBegin;
286: first = 0;
287: last = n - 1;
288: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
290: while (1) {
291: mid = first;
292: d = a[mid];
293: abskey = PetscAbsScalar(d);
294: i = last;
295: for (j = first + 1; j <= i; ++j) {
296: d = a[j];
297: if (PetscAbsScalar(d) >= abskey) {
298: ++mid;
299: /* interchange */
300: tmp = a[mid];
301: itmp = idx[mid];
302: a[mid] = a[j];
303: idx[mid] = idx[j];
304: a[j] = tmp;
305: idx[j] = itmp;
306: }
307: }
309: /* interchange */
310: tmp = a[mid];
311: itmp = idx[mid];
312: a[mid] = a[first];
313: idx[mid] = idx[first];
314: a[first] = tmp;
315: idx[first] = itmp;
317: /* test for while loop */
318: if (mid == ncut) break;
319: else if (mid > ncut) last = mid - 1;
320: else first = mid + 1;
321: }
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
328: Not Collective
330: Input Parameters:
331: + ncut - splitting index
332: - n - number of values to sort
334: Input/Output Parameters:
335: + a - array of values, on output the values are permuted such that its elements satisfy:
336: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
337: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
338: - idx - index for array a, on output permuted accordingly
340: Level: intermediate
342: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
343: @*/
344: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
345: {
346: PetscInt i, mid, last, itmp, j, first;
347: PetscReal d, tmp;
348: PetscReal abskey;
350: PetscFunctionBegin;
351: first = 0;
352: last = n - 1;
353: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
355: while (1) {
356: mid = first;
357: d = a[mid];
358: abskey = PetscAbsReal(d);
359: i = last;
360: for (j = first + 1; j <= i; ++j) {
361: d = a[j];
362: if (PetscAbsReal(d) >= abskey) {
363: ++mid;
364: /* interchange */
365: tmp = a[mid];
366: itmp = idx[mid];
367: a[mid] = a[j];
368: idx[mid] = idx[j];
369: a[j] = tmp;
370: idx[j] = itmp;
371: }
372: }
374: /* interchange */
375: tmp = a[mid];
376: itmp = idx[mid];
377: a[mid] = a[first];
378: idx[mid] = idx[first];
379: a[first] = tmp;
380: idx[first] = itmp;
382: /* test for while loop */
383: if (mid == ncut) break;
384: else if (mid > ncut) last = mid - 1;
385: else first = mid + 1;
386: }
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }