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.
5
6 */
7 #include <petscsys.h> /*I "petscsys.h" I*/
8 #include <petsc/private/petscimpl.h>
9
10 #define SWAP(a, b, t) \
11 do { \
12 t = a; \
13 a = b; \
14 b = t; \
15 } while (0)
16
17 /*@
18 PetscSortedReal - Determines whether the array of `PetscReal` is sorted.
19
20 Not Collective
21
22 Input Parameters:
23 + n - number of values
24 - X - array of integers
25
26 Output Parameter:
27 . sorted - flag whether the array is sorted
28
29 Level: intermediate
30
31 .seealso: `PetscSortReal()`, `PetscSortedInt()`, `PetscSortedMPIInt()`
32 @*/
PetscSortedReal(PetscCount n,const PetscReal X[],PetscBool * sorted)33 PetscErrorCode PetscSortedReal(PetscCount n, const PetscReal X[], PetscBool *sorted)
34 {
35 PetscFunctionBegin;
36 PetscSorted(n, X, *sorted);
37 PetscFunctionReturn(PETSC_SUCCESS);
38 }
39
40 /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
PetscSortReal_Private(PetscReal * v,PetscCount right)41 static PetscErrorCode PetscSortReal_Private(PetscReal *v, PetscCount right)
42 {
43 PetscCount i, last;
44 PetscReal vl, tmp;
45
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 }
67
68 /*@
69 PetscSortReal - Sorts an array of `PetscReal` in place in increasing order.
70
71 Not Collective
72
73 Input Parameters:
74 + n - number of values
75 - v - array of doubles
76
77 Level: intermediate
78
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.
83
84 .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
85 @*/
PetscSortReal(PetscCount n,PetscReal v[])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;
103
104 PetscCall(PetscIntCast(n, &N));
105 PetscCall(PetscSortReal_Private(v, N - 1));
106 }
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
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)
119
120 /* modified from PetscSortIntWithArray_Private */
PetscSortRealWithArrayInt_Private(PetscReal * v,PetscInt * V,PetscCount right)121 static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscCount right)
122 {
123 PetscCount i, last;
124 PetscInt itmp;
125 PetscReal rvl, rtmp;
126
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.
151
152 Not Collective
153
154 Input Parameters:
155 + n - number of values
156 . Ii - array of integers
157 - r - second array of integers
158
159 Level: intermediate
160
161 .seealso: `PetscSortReal()`
162 @*/
PetscSortRealWithArrayInt(PetscCount n,PetscReal r[],PetscInt Ii[])163 PetscErrorCode PetscSortRealWithArrayInt(PetscCount n, PetscReal r[], PetscInt Ii[])
164 {
165 PetscCount j, k;
166 PetscInt itmp;
167 PetscReal rk, rtmp;
168
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 }
187
188 /*@
189 PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s
190
191 Not Collective
192
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
198
199 Output Parameter:
200 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
201
202 Level: intermediate
203
204 .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
205 @*/
PetscFindReal(PetscReal key,PetscCount n,const PetscReal t[],PetscReal eps,PetscInt * loc)206 PetscErrorCode PetscFindReal(PetscReal key, PetscCount n, const PetscReal t[], PetscReal eps, PetscInt *loc)
207 {
208 PetscInt lo = 0, hi;
209
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 }
227
228 /*@
229 PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries
230
231 Not Collective
232
233 Input Parameters:
234 + n - initial number of values
235 - v - array of values
236
237 Note:
238 On output both `n` and `v` are modified with non-redundant values.
239
240 Level: intermediate
241
242 .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
243 @*/
PetscSortRemoveDupsReal(PetscInt * n,PetscReal v[])244 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
245 {
246 PetscInt i, s = 0, N = *n, b = 0;
247
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 }
259
260 /*@
261 PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.
262
263 Not Collective
264
265 Input Parameters:
266 + ncut - splitting index
267 - n - number of values to sort
268
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
274
275 Level: intermediate
276
277 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
278 @*/
PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])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;
284
285 PetscFunctionBegin;
286 first = 0;
287 last = n - 1;
288 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
289
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 }
308
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;
316
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 }
324
325 /*@
326 PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.
327
328 Not Collective
329
330 Input Parameters:
331 + ncut - splitting index
332 - n - number of values to sort
333
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
339
340 Level: intermediate
341
342 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
343 @*/
PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])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;
349
350 PetscFunctionBegin;
351 first = 0;
352 last = n - 1;
353 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);
354
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 }
373
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;
381
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 }
389