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 Parameter:
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: PetscFunctionBegin;
37: PetscSorted(n, X, *sorted);
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
42: static PetscErrorCode PetscSortReal_Private(PetscReal *v, PetscInt right)
43: {
44: PetscInt i, last;
45: PetscReal vl, tmp;
47: PetscFunctionBegin;
48: if (right <= 1) {
49: if (right == 1) {
50: if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
51: }
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
54: SWAP(v[0], v[right / 2], tmp);
55: vl = v[0];
56: last = 0;
57: for (i = 1; i <= right; i++) {
58: if (v[i] < vl) {
59: last++;
60: SWAP(v[last], v[i], tmp);
61: }
62: }
63: SWAP(v[0], v[last], tmp);
64: PetscCall(PetscSortReal_Private(v, last - 1));
65: PetscCall(PetscSortReal_Private(v + last + 1, right - (last + 1)));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: /*@
70: PetscSortReal - Sorts an array of `PetscReal` in place in increasing order.
72: Not Collective
74: Input Parameters:
75: + n - number of values
76: - v - array of doubles
78: Level: intermediate
80: Note:
81: This function serves as an alternative to `PetscRealSortSemiOrdered()`, and may perform faster especially if the array
82: is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
83: code to see which routine is fastest.
85: .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
86: @*/
87: PetscErrorCode PetscSortReal(PetscInt n, PetscReal v[])
88: {
89: PetscInt j, k;
90: PetscReal tmp, vk;
92: PetscFunctionBegin;
94: if (n < 8) {
95: for (k = 0; k < n; k++) {
96: vk = v[k];
97: for (j = k + 1; j < n; j++) {
98: if (vk > v[j]) {
99: SWAP(v[k], v[j], tmp);
100: vk = v[k];
101: }
102: }
103: }
104: } else PetscCall(PetscSortReal_Private(v, n - 1));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: #define SWAP2ri(a, b, c, d, rt, it) \
109: { \
110: rt = a; \
111: a = b; \
112: b = rt; \
113: it = c; \
114: c = d; \
115: d = it; \
116: }
118: /* modified from PetscSortIntWithArray_Private */
119: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscInt right)
120: {
121: PetscInt i, last, itmp;
122: PetscReal rvl, rtmp;
124: PetscFunctionBegin;
125: if (right <= 1) {
126: if (right == 1) {
127: if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
128: }
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
131: SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
132: rvl = v[0];
133: last = 0;
134: for (i = 1; i <= right; i++) {
135: if (v[i] < rvl) {
136: last++;
137: SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
138: }
139: }
140: SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
141: PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1));
142: PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1)));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
145: /*@
146: PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
147: changes a second `PetscInt` array to match the sorted first array.
149: Not Collective
151: Input Parameters:
152: + n - number of values
153: . i - array of integers
154: - I - second array of integers
156: Level: intermediate
158: .seealso: `PetscSortReal()`
159: @*/
160: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n, PetscReal r[], PetscInt Ii[])
161: {
162: PetscInt j, k, itmp;
163: PetscReal rk, rtmp;
165: PetscFunctionBegin;
168: if (n < 8) {
169: for (k = 0; k < n; k++) {
170: rk = r[k];
171: for (j = k + 1; j < n; j++) {
172: if (rk > r[j]) {
173: SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
174: rk = r[k];
175: }
176: }
177: }
178: } else {
179: PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1));
180: }
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@
185: PetscFindReal - Finds a PetscReal` in a sorted array of `PetscReal`s
187: Not Collective
189: Input Parameters:
190: + key - the value to locate
191: . n - number of values in the array
192: . ii - array of values
193: - eps - tolerance used to compare
195: Output Parameter:
196: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
198: Level: intermediate
200: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
201: @*/
202: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
203: {
204: PetscInt lo = 0, hi = n;
206: PetscFunctionBegin;
208: if (!n) {
209: *loc = -1;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscCheckSorted(n, t);
214: while (hi - lo > 1) {
215: PetscInt mid = lo + (hi - lo) / 2;
216: if (key < t[mid]) hi = mid;
217: else lo = mid;
218: }
219: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
226: Not Collective
228: Input Parameters:
229: + n - number of values
230: - v - array of doubles
232: Output Parameter:
233: . n - number of non-redundant values
235: Level: intermediate
237: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
238: @*/
239: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
240: {
241: PetscInt i, s = 0, N = *n, b = 0;
243: PetscFunctionBegin;
244: PetscCall(PetscSortReal(N, v));
245: for (i = 0; i < N - 1; i++) {
246: if (v[b + s + 1] != v[b]) {
247: v[b + 1] = v[b + s + 1];
248: b++;
249: } else s++;
250: }
251: *n = N - s;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
258: Not Collective
260: Input Parameters:
261: + ncut - splitig index
262: - n - number of values to sort
264: Input/Output Parameters:
265: + a - array of values, on output the values are permuted such that its elements satisfy:
266: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
267: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
268: - idx - index for array a, on output permuted accordingly
270: Level: intermediate
272: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
273: @*/
274: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
275: {
276: PetscInt i, mid, last, itmp, j, first;
277: PetscScalar d, tmp;
278: PetscReal abskey;
280: PetscFunctionBegin;
281: first = 0;
282: last = n - 1;
283: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
285: while (1) {
286: mid = first;
287: d = a[mid];
288: abskey = PetscAbsScalar(d);
289: i = last;
290: for (j = first + 1; j <= i; ++j) {
291: d = a[j];
292: if (PetscAbsScalar(d) >= abskey) {
293: ++mid;
294: /* interchange */
295: tmp = a[mid];
296: itmp = idx[mid];
297: a[mid] = a[j];
298: idx[mid] = idx[j];
299: a[j] = tmp;
300: idx[j] = itmp;
301: }
302: }
304: /* interchange */
305: tmp = a[mid];
306: itmp = idx[mid];
307: a[mid] = a[first];
308: idx[mid] = idx[first];
309: a[first] = tmp;
310: idx[first] = itmp;
312: /* test for while loop */
313: if (mid == ncut) break;
314: else if (mid > ncut) last = mid - 1;
315: else first = mid + 1;
316: }
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@
321: PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
323: Not Collective
325: Input Parameters:
326: + ncut - splitig index
327: - n - number of values to sort
329: Input/Output Parameters:
330: + a - array of values, on output the values are permuted such that its elements satisfy:
331: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
332: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
333: - idx - index for array a, on output permuted accordingly
335: Level: intermediate
337: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
338: @*/
339: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
340: {
341: PetscInt i, mid, last, itmp, j, first;
342: PetscReal d, tmp;
343: PetscReal abskey;
345: PetscFunctionBegin;
346: first = 0;
347: last = n - 1;
348: if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
350: while (1) {
351: mid = first;
352: d = a[mid];
353: abskey = PetscAbsReal(d);
354: i = last;
355: for (j = first + 1; j <= i; ++j) {
356: d = a[j];
357: if (PetscAbsReal(d) >= abskey) {
358: ++mid;
359: /* interchange */
360: tmp = a[mid];
361: itmp = idx[mid];
362: a[mid] = a[j];
363: idx[mid] = idx[j];
364: a[j] = tmp;
365: idx[j] = itmp;
366: }
367: }
369: /* interchange */
370: tmp = a[mid];
371: itmp = idx[mid];
372: a[mid] = a[first];
373: idx[mid] = idx[first];
374: a[first] = tmp;
375: idx[first] = itmp;
377: /* test for while loop */
378: if (mid == ncut) break;
379: else if (mid > ncut) last = mid - 1;
380: else first = mid + 1;
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }