Actual source code: sortso.c
1: #include <petsc/private/petscimpl.h>
3: static inline int Compare_PetscMPIInt_Private(const void *left, const void *right, PETSC_UNUSED PetscCtx ctx)
4: {
5: PetscMPIInt l = *(PetscMPIInt *)left, r = *(PetscMPIInt *)right;
6: return (l < r) ? -1 : (l > r);
7: }
9: static inline int Compare_PetscInt_Private(const void *left, const void *right, PETSC_UNUSED PetscCtx ctx)
10: {
11: PetscInt l = *(PetscInt *)left, r = *(PetscInt *)right;
12: return (l < r) ? -1 : (l > r);
13: }
15: static inline int Compare_PetscReal_Private(const void *left, const void *right, PETSC_UNUSED PetscCtx ctx)
16: {
17: PetscReal l = *(PetscReal *)left, r = *(PetscReal *)right;
18: return (l < r) ? -1 : (l > r);
19: }
21: #define MIN_GALLOP_CONST_GLOBAL 8
22: static PetscInt MIN_GALLOP_GLOBAL = MIN_GALLOP_CONST_GLOBAL;
23: typedef int (*CompFunc)(const void *, const void *, void *);
25: /* Mostly to force clang uses the builtins, but can't hurt to have gcc compilers also do the same */
26: #if !PetscDefined(USE_DEBUG) && defined(__GNUC__)
27: static inline void COPYSWAPPY(char *a, char *b, char *t, size_t size)
28: {
29: __builtin_memcpy(t, b, size);
30: __builtin_memmove(b, a, size);
31: __builtin_memcpy(a, t, size);
32: return;
33: }
35: static inline void COPYSWAPPY2(char *al, char *ar, size_t asize, char *bl, char *br, size_t bsize, char *t)
36: {
37: __builtin_memcpy(t, ar, asize);
38: __builtin_memmove(ar, al, asize);
39: __builtin_memcpy(al, t, asize);
40: __builtin_memcpy(t, br, bsize);
41: __builtin_memmove(br, bl, bsize);
42: __builtin_memcpy(bl, t, bsize);
43: return;
44: }
46: static inline void Petsc_memcpy(char *dest, const char *src, size_t size)
47: {
48: __builtin_memcpy(dest, src, size);
49: return;
50: }
52: static inline void Petsc_memcpy2(char *adest, const char *asrc, size_t asize, char *bdest, const char *bsrc, size_t bsize)
53: {
54: __builtin_memcpy(adest, asrc, asize);
55: __builtin_memcpy(bdest, bsrc, bsize);
56: return;
57: }
59: static inline void Petsc_memmove(char *dest, const char *src, size_t size)
60: {
61: __builtin_memmove(dest, src, size);
62: return;
63: }
65: static inline void Petsc_memmove2(char *adest, const char *asrc, size_t asize, char *bdest, const char *bsrc, size_t bsize)
66: {
67: __builtin_memmove(adest, asrc, asize);
68: __builtin_memmove(bdest, bsrc, bsize);
69: return;
70: }
71: #else
72: static inline void COPYSWAPPY(char *a, char *b, char *t, size_t size)
73: {
74: PetscFunctionBegin;
75: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(t, b, size));
76: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(b, a, size));
77: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(a, t, size));
78: PetscFunctionReturnVoid();
79: }
81: static inline void COPYSWAPPY2(char *al, char *ar, size_t asize, char *bl, char *br, size_t bsize, char *t)
82: {
83: PetscFunctionBegin;
84: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(t, ar, asize));
85: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(ar, al, asize));
86: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(al, t, asize));
87: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(t, br, bsize));
88: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(br, bl, bsize));
89: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(bl, t, bsize));
90: PetscFunctionReturnVoid();
91: }
93: static inline void Petsc_memcpy(char *dest, const char *src, size_t size)
94: {
95: PetscFunctionBegin;
96: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(dest, src, size));
97: PetscFunctionReturnVoid();
98: }
100: static inline void Petsc_memcpy2(char *adest, const char *asrc, size_t asize, char *bdest, const char *bsrc, size_t bsize)
101: {
102: PetscFunctionBegin;
103: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(adest, asrc, asize));
104: PetscCallAbort(PETSC_COMM_SELF, PetscMemcpy(bdest, bsrc, bsize));
105: PetscFunctionReturnVoid();
106: }
108: static inline void Petsc_memmove(char *dest, const char *src, size_t size)
109: {
110: PetscFunctionBegin;
111: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(dest, src, size));
112: PetscFunctionReturnVoid();
113: }
115: static inline void Petsc_memmove2(char *adest, const char *asrc, size_t asize, char *bdest, const char *bsrc, size_t bsize)
116: {
117: PetscFunctionBegin;
118: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(adest, asrc, asize));
119: PetscCallAbort(PETSC_COMM_SELF, PetscMemmove(bdest, bsrc, bsize));
120: PetscFunctionReturnVoid();
121: }
122: #endif
124: /* Start left look right. Looking for e.g. B[0] in A or mergelo. l inclusive, r inclusive. Returns first m such that arr[m] >
125: x. Output also inclusive.
127: NOTE: Both gallopsearch functions CANNOT distinguish between inserting AFTER the entire array vs inserting at entry n!! For example for an array:
129: {0,1,2,3,4,5}
131: when looking to insert "5" this routine will return *m = 6, but when looking to insert "6" it will ALSO return *m = 6.
132: */
133: static inline PetscErrorCode PetscGallopSearchLeft_Private(const char *arr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt l, PetscInt r, const char *x, PetscInt *m)
134: {
135: PetscInt last = l, k = 1, mid, cur = l + 1;
137: PetscFunctionBegin;
138: *m = l;
139: PetscAssert(r >= l, PETSC_COMM_SELF, PETSC_ERR_PLIB, "r %" PetscInt_FMT " < l %" PetscInt_FMT " in PetscGallopSearchLeft", r, l);
140: if ((*cmp)(x, arr + r * size, ctx) >= 0) {
141: *m = r;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
144: if ((*cmp)(x, arr + l * size, ctx) < 0 || PetscUnlikely(!(r - l))) PetscFunctionReturn(PETSC_SUCCESS);
145: while (PETSC_TRUE) {
146: if (PetscUnlikely(cur > r)) {
147: cur = r;
148: break;
149: }
150: if ((*cmp)(x, arr + cur * size, ctx) < 0) break;
151: last = cur;
152: cur += (k <<= 1) + 1;
153: ++k;
154: }
155: /* standard binary search but take last 0 mid 0 cur 1 into account*/
156: while (cur > last + 1) {
157: mid = last + ((cur - last) >> 1);
158: if ((*cmp)(x, arr + mid * size, ctx) < 0) {
159: cur = mid;
160: } else {
161: last = mid;
162: }
163: }
164: *m = cur;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /* Start right look left. Looking for e.g. A[-1] in B or mergehi. l inclusive, r inclusive. Returns last m such that arr[m]
169: < x. Output also inclusive */
170: static inline PetscErrorCode PetscGallopSearchRight_Private(const char *arr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt l, PetscInt r, const char *x, PetscInt *m)
171: {
172: PetscInt last = r, k = 1, mid, cur = r - 1;
174: PetscFunctionBegin;
175: *m = r;
176: PetscAssert(r >= l, PETSC_COMM_SELF, PETSC_ERR_PLIB, "r %" PetscInt_FMT " < l %" PetscInt_FMT " in PetscGallopSearchRight", r, l);
177: if ((*cmp)(x, arr + l * size, ctx) <= 0) {
178: *m = l;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
181: if ((*cmp)(x, arr + r * size, ctx) > 0 || PetscUnlikely(!(r - l))) PetscFunctionReturn(PETSC_SUCCESS);
182: while (PETSC_TRUE) {
183: if (PetscUnlikely(cur < l)) {
184: cur = l;
185: break;
186: }
187: if ((*cmp)(x, arr + cur * size, ctx) > 0) break;
188: last = cur;
189: cur -= (k <<= 1) + 1;
190: ++k;
191: }
192: /* standard binary search but take last r-1 mid r-1 cur r-2 into account*/
193: while (last > cur + 1) {
194: mid = last - ((last - cur) >> 1);
195: if ((*cmp)(x, arr + mid * size, ctx) > 0) {
196: cur = mid;
197: } else {
198: last = mid;
199: }
200: }
201: *m = cur;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /* Mergesort where size of left half <= size of right half, so mergesort is done left to right. Arr should be pointer to
206: complete array, left is first index of left array, mid is first index of right array, right is last index of right
207: array */
208: static inline PetscErrorCode PetscTimSortMergeLo_Private(char *arr, char *tarr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt mid, PetscInt right)
209: {
210: PetscInt i = 0, j = mid, k = left, gallopleft = 0, gallopright = 0;
211: const PetscInt llen = mid - left;
213: PetscFunctionBegin;
214: Petsc_memcpy(tarr, arr + (left * size), llen * size);
215: while ((i < llen) && (j <= right)) {
216: if ((*cmp)(tarr + (i * size), arr + (j * size), ctx) < 0) {
217: Petsc_memcpy(arr + (k * size), tarr + (i * size), size);
218: ++k;
219: gallopright = 0;
220: if (++i < llen && ++gallopleft >= MIN_GALLOP_GLOBAL) {
221: PetscInt l1, l2, diff1, diff2;
222: do {
223: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
224: /* search temp for right[j], can move up to that of temp into arr immediately */
225: PetscCall(PetscGallopSearchLeft_Private(tarr, size, cmp, ctx, i, llen - 1, arr + (j * size), &l1));
226: diff1 = l1 - i;
227: Petsc_memcpy(arr + (k * size), tarr + (i * size), diff1 * size);
228: k += diff1;
229: i = l1;
230: /* search right for temp[i], can move up to that many of right into arr */
231: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, j, right, tarr + (i * size), &l2));
232: diff2 = l2 - j;
233: Petsc_memmove(arr + k * size, arr + j * size, diff2 * size);
234: k += diff2;
235: j = l2;
236: if (i >= llen || j > right) break;
237: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
238: ++MIN_GALLOP_GLOBAL;
239: }
240: } else {
241: Petsc_memmove(arr + (k * size), arr + (j * size), size);
242: ++k;
243: gallopleft = 0;
244: if (++j <= right && ++gallopright >= MIN_GALLOP_GLOBAL) {
245: PetscInt l1, l2, diff1, diff2;
246: do {
247: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
248: /* search right for temp[i], can move up to that many of right into arr */
249: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, j, right, tarr + (i * size), &l2));
250: diff2 = l2 - j;
251: Petsc_memmove(arr + (k * size), arr + (j * size), diff2 * size);
252: k += diff2;
253: j = l2;
254: /* search temp for right[j], can copy up to that of temp into arr immediately */
255: PetscCall(PetscGallopSearchLeft_Private(tarr, size, cmp, ctx, i, llen - 1, arr + (j * size), &l1));
256: diff1 = l1 - i;
257: Petsc_memcpy(arr + (k * size), tarr + (i * size), diff1 * size);
258: k += diff1;
259: i = l1;
260: if (i >= llen || j > right) break;
261: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
262: ++MIN_GALLOP_GLOBAL;
263: }
264: }
265: }
266: if (i < llen) Petsc_memcpy(arr + (k * size), tarr + (i * size), (llen - i) * size);
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static inline PetscErrorCode PetscTimSortMergeLoWithArray_Private(char *arr, char *atarr, size_t asize, char *barr, char *btarr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt mid, PetscInt right)
271: {
272: PetscInt i = 0, j = mid, k = left, gallopleft = 0, gallopright = 0;
273: const PetscInt llen = mid - left;
275: PetscFunctionBegin;
276: Petsc_memcpy2(atarr, arr + (left * asize), llen * asize, btarr, barr + (left * bsize), llen * bsize);
277: while ((i < llen) && (j <= right)) {
278: if ((*cmp)(atarr + (i * asize), arr + (j * asize), ctx) < 0) {
279: Petsc_memcpy2(arr + (k * asize), atarr + (i * asize), asize, barr + (k * bsize), btarr + (i * bsize), bsize);
280: ++k;
281: gallopright = 0;
282: if (++i < llen && ++gallopleft >= MIN_GALLOP_GLOBAL) {
283: PetscInt l1, l2, diff1, diff2;
284: do {
285: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
286: /* search temp for right[j], can move up to that of temp into arr immediately */
287: PetscCall(PetscGallopSearchLeft_Private(atarr, asize, cmp, ctx, i, llen - 1, arr + (j * asize), &l1));
288: diff1 = l1 - i;
289: Petsc_memcpy2(arr + (k * asize), atarr + (i * asize), diff1 * asize, barr + (k * bsize), btarr + (i * bsize), diff1 * bsize);
290: k += diff1;
291: i = l1;
292: /* search right for temp[i], can move up to that many of right into arr */
293: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, j, right, atarr + (i * asize), &l2));
294: diff2 = l2 - j;
295: Petsc_memmove2(arr + (k * asize), arr + (j * asize), diff2 * asize, barr + (k * bsize), barr + (j * bsize), diff2 * bsize);
296: k += diff2;
297: j = l2;
298: if (i >= llen || j > right) break;
299: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
300: ++MIN_GALLOP_GLOBAL;
301: }
302: } else {
303: Petsc_memmove2(arr + (k * asize), arr + (j * asize), asize, barr + (k * bsize), barr + (j * bsize), bsize);
304: ++k;
305: gallopleft = 0;
306: if (++j <= right && ++gallopright >= MIN_GALLOP_GLOBAL) {
307: PetscInt l1, l2, diff1, diff2;
308: do {
309: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
310: /* search right for temp[i], can move up to that many of right into arr */
311: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, j, right, atarr + (i * asize), &l2));
312: diff2 = l2 - j;
313: Petsc_memmove2(arr + (k * asize), arr + (j * asize), diff2 * asize, barr + (k * bsize), barr + (j * bsize), diff2 * bsize);
314: k += diff2;
315: j = l2;
316: /* search temp for right[j], can copy up to that of temp into arr immediately */
317: PetscCall(PetscGallopSearchLeft_Private(atarr, asize, cmp, ctx, i, llen - 1, arr + (j * asize), &l1));
318: diff1 = l1 - i;
319: Petsc_memcpy2(arr + (k * asize), atarr + (i * asize), diff1 * asize, barr + (k * bsize), btarr + (i * bsize), diff1 * bsize);
320: k += diff1;
321: i = l1;
322: if (i >= llen || j > right) break;
323: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
324: ++MIN_GALLOP_GLOBAL;
325: }
326: }
327: }
328: if (i < llen) Petsc_memcpy2(arr + (k * asize), atarr + (i * asize), (llen - i) * asize, barr + (k * bsize), btarr + (i * bsize), (llen - i) * bsize);
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /* Mergesort where size of right half < size of left half, so mergesort is done right to left. Arr should be pointer to
333: complete array, left is first index of left array, mid is first index of right array, right is last index of right
334: array */
335: static inline PetscErrorCode PetscTimSortMergeHi_Private(char *arr, char *tarr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt mid, PetscInt right)
336: {
337: PetscInt i = right - mid, j = mid - 1, k = right, gallopleft = 0, gallopright = 0;
338: const PetscInt rlen = right - mid + 1;
340: PetscFunctionBegin;
341: Petsc_memcpy(tarr, arr + mid * size, rlen * size);
342: while ((i >= 0) && (j >= left)) {
343: if ((*cmp)(tarr + i * size, arr + j * size, ctx) > 0) {
344: Petsc_memcpy(arr + k * size, tarr + i * size, size);
345: --k;
346: gallopleft = 0;
347: if (--i >= 0 && ++gallopright >= MIN_GALLOP_GLOBAL) {
348: PetscInt l1, l2, diff1, diff2;
349: do {
350: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
351: /* search temp for left[j], can copy up to that many of temp into arr */
352: PetscCall(PetscGallopSearchRight_Private(tarr, size, cmp, ctx, 0, i, (arr) + j * size, &l1));
353: diff1 = i - l1;
354: Petsc_memcpy(arr + (k - diff1 + 1) * size, tarr + (l1 + 1) * size, diff1 * size);
355: k -= diff1;
356: i = l1;
357: /* search left for temp[i], can move up to that many of left up arr */
358: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, left, j, (tarr) + i * size, &l2));
359: diff2 = j - l2;
360: Petsc_memmove(arr + (k - diff2 + 1) * size, arr + (l2 + 1) * size, diff2 * size);
361: k -= diff2;
362: j = l2;
363: if (i < 0 || j < left) break;
364: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
365: ++MIN_GALLOP_GLOBAL;
366: }
367: } else {
368: Petsc_memmove(arr + k * size, arr + j * size, size);
369: --k;
370: gallopright = 0;
371: if (--j >= left && ++gallopleft >= MIN_GALLOP_GLOBAL) {
372: PetscInt l1, l2, diff1, diff2;
373: do {
374: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
375: /* search left for temp[i], can move up to that many of left up arr */
376: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, left, j, (tarr) + i * size, &l2));
377: diff2 = j - l2;
378: Petsc_memmove(arr + (k - diff2 + 1) * size, arr + (l2 + 1) * size, diff2 * size);
379: k -= diff2;
380: j = l2;
381: /* search temp for left[j], can copy up to that many of temp into arr */
382: PetscCall(PetscGallopSearchRight_Private(tarr, size, cmp, ctx, 0, i, (arr) + j * size, &l1));
383: diff1 = i - l1;
384: Petsc_memcpy(arr + (k - diff1 + 1) * size, tarr + (l1 + 1) * size, diff1 * size);
385: k -= diff1;
386: i = l1;
387: if (i < 0 || j < left) break;
388: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
389: ++MIN_GALLOP_GLOBAL;
390: }
391: }
392: }
393: if (i >= 0) Petsc_memcpy(arr + left * size, tarr, (i + 1) * size);
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static inline PetscErrorCode PetscTimSortMergeHiWithArray_Private(char *arr, char *atarr, size_t asize, char *barr, char *btarr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt mid, PetscInt right)
398: {
399: PetscInt i = right - mid, j = mid - 1, k = right, gallopleft = 0, gallopright = 0;
400: const PetscInt rlen = right - mid + 1;
402: PetscFunctionBegin;
403: Petsc_memcpy2(atarr, arr + mid * asize, rlen * asize, btarr, barr + mid * bsize, rlen * bsize);
404: while ((i >= 0) && (j >= left)) {
405: if ((*cmp)(atarr + i * asize, arr + j * asize, ctx) > 0) {
406: Petsc_memcpy2(arr + k * asize, atarr + i * asize, asize, barr + k * bsize, btarr + i * bsize, bsize);
407: --k;
408: gallopleft = 0;
409: if (--i >= 0 && ++gallopright >= MIN_GALLOP_GLOBAL) {
410: PetscInt l1, l2, diff1, diff2;
411: do {
412: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
413: /* search temp for left[j], can copy up to that many of temp into arr */
414: PetscCall(PetscGallopSearchRight_Private(atarr, asize, cmp, ctx, 0, i, (arr) + j * asize, &l1));
415: diff1 = i - l1;
416: Petsc_memcpy2(arr + (k - diff1 + 1) * asize, atarr + (l1 + 1) * asize, diff1 * asize, barr + (k - diff1 + 1) * bsize, btarr + (l1 + 1) * bsize, diff1 * bsize);
417: k -= diff1;
418: i = l1;
419: /* search left for temp[i], can move up to that many of left up arr */
420: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, left, j, (atarr) + i * asize, &l2));
421: diff2 = j - l2;
422: Petsc_memmove2(arr + (k - diff2 + 1) * asize, arr + (l2 + 1) * asize, diff2 * asize, barr + (k - diff2 + 1) * bsize, barr + (l2 + 1) * bsize, diff2 * bsize);
423: k -= diff2;
424: j = l2;
425: if (i < 0 || j < left) break;
426: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
427: ++MIN_GALLOP_GLOBAL;
428: }
429: } else {
430: Petsc_memmove2(arr + k * asize, arr + j * asize, asize, barr + k * bsize, barr + j * bsize, bsize);
431: --k;
432: gallopright = 0;
433: if (--j >= left && ++gallopleft >= MIN_GALLOP_GLOBAL) {
434: PetscInt l1, l2, diff1, diff2;
435: do {
436: if (MIN_GALLOP_GLOBAL > 1) --MIN_GALLOP_GLOBAL;
437: /* search left for temp[i], can move up to that many of left up arr */
438: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, left, j, (atarr) + i * asize, &l2));
439: diff2 = j - l2;
440: Petsc_memmove2(arr + (k - diff2 + 1) * asize, arr + (l2 + 1) * asize, diff2 * asize, barr + (k - diff2 + 1) * bsize, barr + (l2 + 1) * bsize, diff2 * bsize);
441: k -= diff2;
442: j = l2;
443: /* search temp for left[j], can copy up to that many of temp into arr */
444: PetscCall(PetscGallopSearchRight_Private(atarr, asize, cmp, ctx, 0, i, (arr) + j * asize, &l1));
445: diff1 = i - l1;
446: Petsc_memcpy2(arr + (k - diff1 + 1) * asize, atarr + (l1 + 1) * asize, diff1 * asize, barr + (k - diff1 + 1) * bsize, btarr + (l1 + 1) * bsize, diff1 * bsize);
447: k -= diff1;
448: i = l1;
449: if (i < 0 || j < left) break;
450: } while (diff1 > MIN_GALLOP_GLOBAL || diff2 > MIN_GALLOP_GLOBAL);
451: ++MIN_GALLOP_GLOBAL;
452: }
453: }
454: }
455: if (i >= 0) Petsc_memcpy2(arr + left * asize, atarr, (i + 1) * asize, barr + left * bsize, btarr, (i + 1) * bsize);
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /* Left is inclusive lower bound of array slice, start is start location of unsorted section, right is inclusive upper
460: bound of array slice. If unsure of where unsorted section starts or if entire length is unsorted pass start = left */
461: static inline PetscErrorCode PetscInsertionSort_Private(char *arr, char *tarr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt start, PetscInt right)
462: {
463: PetscInt i = start == left ? start + 1 : start;
465: PetscFunctionBegin;
466: for (; i <= right; ++i) {
467: PetscInt j = i - 1;
468: Petsc_memcpy(tarr, arr + (i * size), size);
469: while ((j >= left) && ((*cmp)(tarr, arr + j * size, ctx) < 0)) {
470: COPYSWAPPY(arr + (j + 1) * size, arr + j * size, tarr + size, size);
471: --j;
472: }
473: Petsc_memcpy(arr + (j + 1) * size, tarr, size);
474: }
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static inline PetscErrorCode PetscInsertionSortWithArray_Private(char *arr, char *atarr, size_t asize, char *barr, char *btarr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt start, PetscInt right)
479: {
480: PetscInt i = start == left ? start + 1 : start;
482: PetscFunctionBegin;
483: for (; i <= right; ++i) {
484: PetscInt j = i - 1;
485: Petsc_memcpy2(atarr, arr + (i * asize), asize, btarr, barr + (i * bsize), bsize);
486: while ((j >= left) && ((*cmp)(atarr, arr + (j * asize), ctx) < 0)) {
487: COPYSWAPPY2(arr + (j + 1) * asize, arr + j * asize, asize, barr + (j + 1) * bsize, barr + j * bsize, bsize, atarr + asize);
488: --j;
489: }
490: Petsc_memcpy2(arr + (j + 1) * asize, atarr, asize, barr + (j + 1) * bsize, btarr, bsize);
491: }
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /* See PetscInsertionSort_Private */
496: static inline PetscErrorCode PetscBinaryInsertionSort_Private(char *arr, char *tarr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt start, PetscInt right)
497: {
498: PetscInt i = start == left ? start + 1 : start;
500: PetscFunctionBegin;
501: for (; i <= right; ++i) {
502: PetscInt l = left, r = i;
503: Petsc_memcpy(tarr, arr + (i * size), size);
504: do {
505: const PetscInt m = l + ((r - l) >> 1);
506: if ((*cmp)(tarr, arr + (m * size), ctx) < 0) {
507: r = m;
508: } else {
509: l = m + 1;
510: }
511: } while (l < r);
512: Petsc_memmove(arr + ((l + 1) * size), arr + (l * size), (i - l) * size);
513: Petsc_memcpy(arr + (l * size), tarr, size);
514: }
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static inline PetscErrorCode PetscBinaryInsertionSortWithArray_Private(char *arr, char *atarr, size_t asize, char *barr, char *btarr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscInt left, PetscInt start, PetscInt right)
519: {
520: PetscInt i = start == left ? start + 1 : start;
522: PetscFunctionBegin;
523: for (; i <= right; ++i) {
524: PetscInt l = left, r = i;
525: Petsc_memcpy2(atarr, arr + (i * asize), asize, btarr, barr + (i * bsize), bsize);
526: do {
527: const PetscInt m = l + ((r - l) >> 1);
528: if ((*cmp)(atarr, arr + (m * asize), ctx) < 0) {
529: r = m;
530: } else {
531: l = m + 1;
532: }
533: } while (l < r);
534: Petsc_memmove2(arr + ((l + 1) * asize), arr + (l * asize), (i - l) * asize, barr + ((l + 1) * bsize), barr + (l * bsize), (i - l) * bsize);
535: Petsc_memcpy2(arr + (l * asize), atarr, asize, barr + (l * bsize), btarr, bsize);
536: }
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: typedef struct {
541: PetscInt size;
542: PetscInt start;
543: #if defined(__CRAY_AARCH64) /* segfaults with Cray compilers for aarch64 on a64FX */
544: } PetscTimSortStack;
545: #else
546: } PetscTimSortStack PETSC_ATTRIBUTEALIGNED(2 * sizeof(PetscInt));
547: #endif
549: typedef struct {
550: char *ptr PETSC_ATTRIBUTEALIGNED(PETSC_MEMALIGN);
551: size_t size;
552: size_t maxsize;
553: } PetscTimSortBuffer;
555: static inline PetscErrorCode PetscTimSortResizeBuffer_Private(PetscTimSortBuffer *buff, size_t newSize)
556: {
557: PetscFunctionBegin;
558: if (PetscLikely(newSize <= buff->size)) PetscFunctionReturn(PETSC_SUCCESS);
559: {
560: /* Can't be larger than n, there is merit to simply allocating buff to n to begin with */
561: size_t newMax = PetscMin(newSize * newSize, buff->maxsize);
562: PetscCall(PetscFree(buff->ptr));
563: PetscCall(PetscMalloc1(newMax, &buff->ptr));
564: buff->size = newMax;
565: }
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static inline PetscErrorCode PetscTimSortForceCollapse_Private(char *arr, size_t size, CompFunc cmp, PetscCtx ctx, PetscTimSortBuffer *buff, PetscTimSortStack *stack, PetscInt stacksize)
570: {
571: PetscFunctionBegin;
572: for (; stacksize; --stacksize) {
573: /* A = stack[i-1], B = stack[i], if A[-1] <= B[0] means sorted */
574: if ((*cmp)(arr + (stack[stacksize].start - 1) * size, arr + stack[stacksize].start * size, ctx) > 0) {
575: PetscInt l, m = stack[stacksize].start, r;
576: /* Search A for B[0] insertion */
577: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, stack[stacksize - 1].start, stack[stacksize].start - 1, (arr) + (stack[stacksize].start) * size, &l));
578: /* Search B for A[-1] insertion */
579: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, stack[stacksize].start, stack[stacksize].start + stack[stacksize].size - 1, (arr) + (stack[stacksize].start - 1) * size, &r));
580: if (m - l <= r - m) {
581: PetscCall(PetscTimSortResizeBuffer_Private(buff, (m - l + 1) * size));
582: PetscCall(PetscTimSortMergeLo_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
583: } else {
584: PetscCall(PetscTimSortResizeBuffer_Private(buff, (r - m + 1) * size));
585: PetscCall(PetscTimSortMergeHi_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
586: }
587: }
588: /* Update A with merge */
589: stack[stacksize - 1].size += stack[stacksize].size;
590: }
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: static inline PetscErrorCode PetscTimSortForceCollapseWithArray_Private(char *arr, size_t asize, char *barr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscTimSortBuffer *abuff, PetscTimSortBuffer *bbuff, PetscTimSortStack *stack, PetscInt stacksize)
595: {
596: PetscFunctionBegin;
597: for (; stacksize; --stacksize) {
598: /* A = stack[i-1], B = stack[i], if A[-1] <= B[0] means sorted */
599: if ((*cmp)(arr + (stack[stacksize].start - 1) * asize, arr + stack[stacksize].start * asize, ctx) > 0) {
600: PetscInt l, m = stack[stacksize].start, r;
601: /* Search A for B[0] insertion */
602: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, stack[stacksize - 1].start, stack[stacksize].start - 1, (arr) + (stack[stacksize].start) * asize, &l));
603: /* Search B for A[-1] insertion */
604: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, stack[stacksize].start, stack[stacksize].start + stack[stacksize].size - 1, (arr) + (stack[stacksize].start - 1) * asize, &r));
605: if (m - l <= r - m) {
606: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (m - l + 1) * asize));
607: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (m - l + 1) * bsize));
608: PetscCall(PetscTimSortMergeLoWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
609: } else {
610: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (r - m + 1) * asize));
611: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (r - m + 1) * bsize));
612: PetscCall(PetscTimSortMergeHiWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
613: }
614: }
615: /* Update A with merge */
616: stack[stacksize - 1].size += stack[stacksize].size;
617: }
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: static inline PetscErrorCode PetscTimSortMergeCollapse_Private(char *arr, size_t size, CompFunc cmp, PetscCtx ctx, PetscTimSortBuffer *buff, PetscTimSortStack *stack, PetscInt *stacksize)
622: {
623: PetscInt i = *stacksize;
625: PetscFunctionBegin;
626: while (i) {
627: PetscInt l, m, r, itemp = i;
629: if (i == 1) {
630: /* A = stack[i-1], B = stack[i] */
631: if (stack[i - 1].size < stack[i].size) {
632: /* if A[-1] <= B[0] then sorted */
633: if ((*cmp)(arr + (stack[i].start - 1) * size, arr + stack[i].start * size, ctx) > 0) {
634: m = stack[i].start;
635: /* Search A for B[0] insertion */
636: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, stack[i - 1].start, stack[i].start - 1, (arr) + stack[i].start * size, &l));
637: /* Search B for A[-1] insertion */
638: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, stack[i].start, stack[i].start + stack[i].size - 1, arr + (stack[i].start - 1) * size, &r));
639: if (m - l <= r - m) {
640: PetscCall(PetscTimSortResizeBuffer_Private(buff, (m - l + 1) * size));
641: PetscCall(PetscTimSortMergeLo_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
642: } else {
643: PetscCall(PetscTimSortResizeBuffer_Private(buff, (r - m + 1) * size));
644: PetscCall(PetscTimSortMergeHi_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
645: }
646: }
647: /* Update A with merge */
648: stack[i - 1].size += stack[i].size;
649: --i;
650: }
651: } else {
652: /* i > 2, i.e. C exists
653: A = stack[i-2], B = stack[i-1], C = stack[i]; */
654: if (stack[i - 2].size <= stack[i - 1].size + stack[i].size) {
655: if (stack[i - 2].size < stack[i].size) {
656: /* merge B into A, but if A[-1] <= B[0] then already sorted */
657: if ((*cmp)(arr + (stack[i - 1].start - 1) * size, arr + stack[i - 1].start * size, ctx) > 0) {
658: m = stack[i - 1].start;
659: /* Search A for B[0] insertion */
660: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, stack[i - 2].start, stack[i - 1].start - 1, (arr) + (stack[i - 1].start) * size, &l));
661: /* Search B for A[-1] insertion */
662: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, stack[i - 1].start, stack[i - 1].start + stack[i - 1].size - 1, (arr) + (stack[i - 1].start - 1) * size, &r));
663: if (m - l <= r - m) {
664: PetscCall(PetscTimSortResizeBuffer_Private(buff, (m - l + 1) * size));
665: PetscCall(PetscTimSortMergeLo_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
666: } else {
667: PetscCall(PetscTimSortResizeBuffer_Private(buff, (r - m + 1) * size));
668: PetscCall(PetscTimSortMergeHi_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
669: }
670: }
671: /* Update A with merge */
672: stack[i - 2].size += stack[i - 1].size;
673: /* Push C up the stack */
674: stack[i - 1].start = stack[i].start;
675: stack[i - 1].size = stack[i].size;
676: } else {
677: /* merge C into B */
678: mergeBC:
679: /* If B[-1] <= C[0] then... you know the drill */
680: if ((*cmp)(arr + (stack[i].start - 1) * size, arr + stack[i].start * size, ctx) > 0) {
681: m = stack[i].start;
682: /* Search B for C[0] insertion */
683: PetscCall(PetscGallopSearchLeft_Private(arr, size, cmp, ctx, stack[i - 1].start, stack[i].start - 1, arr + stack[i].start * size, &l));
684: /* Search C for B[-1] insertion */
685: PetscCall(PetscGallopSearchRight_Private(arr, size, cmp, ctx, stack[i].start, stack[i].start + stack[i].size - 1, (arr) + (stack[i].start - 1) * size, &r));
686: if (m - l <= r - m) {
687: PetscCall(PetscTimSortResizeBuffer_Private(buff, (m - l + 1) * size));
688: PetscCall(PetscTimSortMergeLo_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
689: } else {
690: PetscCall(PetscTimSortResizeBuffer_Private(buff, (r - m + 1) * size));
691: PetscCall(PetscTimSortMergeHi_Private(arr, buff->ptr, size, cmp, ctx, l, m, r));
692: }
693: }
694: /* Update B with merge */
695: stack[i - 1].size += stack[i].size;
696: }
697: --i;
698: } else if (stack[i - 1].size <= stack[i].size) {
699: /* merge C into B */
700: goto mergeBC;
701: }
702: }
703: if (itemp == i) break;
704: }
705: *stacksize = i;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static inline PetscErrorCode PetscTimSortMergeCollapseWithArray_Private(char *arr, size_t asize, char *barr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscTimSortBuffer *abuff, PetscTimSortBuffer *bbuff, PetscTimSortStack *stack, PetscInt *stacksize)
710: {
711: PetscInt i = *stacksize;
713: PetscFunctionBegin;
714: while (i) {
715: PetscInt l, m, r, itemp = i;
717: if (i == 1) {
718: /* A = stack[i-1], B = stack[i] */
719: if (stack[i - 1].size < stack[i].size) {
720: /* if A[-1] <= B[0] then sorted */
721: if ((*cmp)(arr + (stack[i].start - 1) * asize, arr + stack[i].start * asize, ctx) > 0) {
722: m = stack[i].start;
723: /* Search A for B[0] insertion */
724: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, stack[i - 1].start, stack[i].start - 1, (arr) + stack[i].start * asize, &l));
725: /* Search B for A[-1] insertion */
726: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, stack[i].start, stack[i].start + stack[i].size - 1, arr + (stack[i].start - 1) * asize, &r));
727: if (m - l <= r - m) {
728: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (m - l + 1) * asize));
729: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (m - l + 1) * bsize));
730: PetscCall(PetscTimSortMergeLoWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
731: } else {
732: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (r - m + 1) * asize));
733: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (r - m + 1) * bsize));
734: PetscCall(PetscTimSortMergeHiWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
735: }
736: }
737: /* Update A with merge */
738: stack[i - 1].size += stack[i].size;
739: --i;
740: }
741: } else {
742: /* i > 2, i.e. C exists
743: A = stack[i-2], B = stack[i-1], C = stack[i]; */
744: if (stack[i - 2].size <= stack[i - 1].size + stack[i].size) {
745: if (stack[i - 2].size < stack[i].size) {
746: /* merge B into A, but if A[-1] <= B[0] then already sorted */
747: if ((*cmp)(arr + (stack[i - 1].start - 1) * asize, arr + stack[i - 1].start * asize, ctx) > 0) {
748: m = stack[i - 1].start;
749: /* Search A for B[0] insertion */
750: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, stack[i - 2].start, stack[i - 1].start - 1, (arr) + (stack[i - 1].start) * asize, &l));
751: /* Search B for A[-1] insertion */
752: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, stack[i - 1].start, stack[i - 1].start + stack[i - 1].size - 1, (arr) + (stack[i - 1].start - 1) * asize, &r));
753: if (m - l <= r - m) {
754: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (m - l + 1) * asize));
755: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (m - l + 1) * bsize));
756: PetscCall(PetscTimSortMergeLoWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
757: } else {
758: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (r - m + 1) * asize));
759: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (r - m + 1) * bsize));
760: PetscCall(PetscTimSortMergeHiWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
761: }
762: }
763: /* Update A with merge */
764: stack[i - 2].size += stack[i - 1].size;
765: /* Push C up the stack */
766: stack[i - 1].start = stack[i].start;
767: stack[i - 1].size = stack[i].size;
768: } else {
769: /* merge C into B */
770: mergeBC:
771: /* If B[-1] <= C[0] then... you know the drill */
772: if ((*cmp)(arr + (stack[i].start - 1) * asize, arr + stack[i].start * asize, ctx) > 0) {
773: m = stack[i].start;
774: /* Search B for C[0] insertion */
775: PetscCall(PetscGallopSearchLeft_Private(arr, asize, cmp, ctx, stack[i - 1].start, stack[i].start - 1, arr + stack[i].start * asize, &l));
776: /* Search C for B[-1] insertion */
777: PetscCall(PetscGallopSearchRight_Private(arr, asize, cmp, ctx, stack[i].start, stack[i].start + stack[i].size - 1, (arr) + (stack[i].start - 1) * asize, &r));
778: if (m - l <= r - m) {
779: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (m - l + 1) * asize));
780: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (m - l + 1) * bsize));
781: PetscCall(PetscTimSortMergeLoWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
782: } else {
783: PetscCall(PetscTimSortResizeBuffer_Private(abuff, (r - m + 1) * asize));
784: PetscCall(PetscTimSortResizeBuffer_Private(bbuff, (r - m + 1) * bsize));
785: PetscCall(PetscTimSortMergeHiWithArray_Private(arr, abuff->ptr, asize, barr, bbuff->ptr, bsize, cmp, ctx, l, m, r));
786: }
787: }
788: /* Update B with merge */
789: stack[i - 1].size += stack[i].size;
790: }
791: --i;
792: } else if (stack[i - 1].size <= stack[i].size) {
793: /* merge C into B */
794: goto mergeBC;
795: }
796: }
797: if (itemp == i) break;
798: }
799: *stacksize = i;
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: /* March sequentially through the array building up a "run" of weakly increasing or strictly decreasing contiguous
804: elements. Decreasing runs are reversed by swapping. If the run is less than minrun, artificially extend it via either
805: binary insertion sort or regulat insertion sort */
806: static inline PetscErrorCode PetscTimSortBuildRun_Private(char *arr, char *tarr, size_t size, CompFunc cmp, PetscCtx ctx, PetscInt n, PetscInt minrun, PetscInt runstart, PetscInt *runend)
807: {
808: const PetscInt re = PetscMin(runstart + minrun, n - 1);
809: PetscInt ri = runstart;
811: PetscFunctionBegin;
812: if (PetscUnlikely(runstart == n - 1)) {
813: *runend = runstart;
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
816: /* guess whether run is ascending or descending and tally up the longest consecutive run. essentially a coinflip for random data */
817: if ((*cmp)(arr + (ri + 1) * size, arr + ri * size, ctx) < 0) {
818: ++ri;
819: while (ri < n - 1) {
820: if ((*cmp)(arr + (ri + 1) * size, arr + ri * size, ctx) >= 0) break;
821: ++ri;
822: }
823: {
824: PetscInt lo = runstart, hi = ri;
825: for (; lo < hi; ++lo, --hi) COPYSWAPPY(arr + lo * size, arr + hi * size, tarr, size);
826: }
827: } else {
828: ++ri;
829: while (ri < n - 1) {
830: if ((*cmp)(arr + (ri + 1) * size, arr + ri * size, ctx) < 0) break;
831: ++ri;
832: }
833: }
834: if (ri < re) {
835: /* the attempt failed, this section likely contains random data. If ri got close to minrun (within 50%) then we try
836: binary search */
837: if (ri - runstart <= minrun >> 1) {
838: ++MIN_GALLOP_GLOBAL; /* didn't get close hedge our bets against random data */
839: PetscCall(PetscInsertionSort_Private(arr, tarr, size, cmp, ctx, runstart, ri, re));
840: } else {
841: PetscCall(PetscBinaryInsertionSort_Private(arr, tarr, size, cmp, ctx, runstart, ri, re));
842: }
843: *runend = re;
844: } else *runend = ri;
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static inline PetscErrorCode PetscTimSortBuildRunWithArray_Private(char *arr, char *atarr, size_t asize, char *barr, char *btarr, size_t bsize, CompFunc cmp, PetscCtx ctx, PetscInt n, PetscInt minrun, PetscInt runstart, PetscInt *runend)
849: {
850: const PetscInt re = PetscMin(runstart + minrun, n - 1);
851: PetscInt ri = runstart;
853: PetscFunctionBegin;
854: if (PetscUnlikely(runstart == n - 1)) {
855: *runend = runstart;
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
858: /* guess whether run is ascending or descending and tally up the longest consecutive run. essentially a coinflip for random data */
859: if ((*cmp)(arr + (ri + 1) * asize, arr + (ri * asize), ctx) < 0) {
860: ++ri;
861: while (ri < n - 1) {
862: if ((*cmp)(arr + (ri + 1) * asize, arr + ri * asize, ctx) >= 0) break;
863: ++ri;
864: }
865: {
866: PetscInt lo = runstart, hi = ri;
867: for (; lo < hi; ++lo, --hi) COPYSWAPPY2(arr + lo * asize, arr + hi * asize, asize, barr + lo * bsize, barr + hi * bsize, bsize, atarr);
868: }
869: } else {
870: ++ri;
871: while (ri < n - 1) {
872: if ((*cmp)(arr + (ri + 1) * asize, arr + ri * asize, ctx) < 0) break;
873: ++ri;
874: }
875: }
876: if (ri < re) {
877: /* the attempt failed, this section likely contains random data. If ri got close to minrun (within 50%) then we try
878: binary search */
879: if (ri - runstart <= minrun >> 1) {
880: ++MIN_GALLOP_GLOBAL; /* didn't get close hedge our bets against random data */
881: PetscCall(PetscInsertionSortWithArray_Private(arr, atarr, asize, barr, btarr, bsize, cmp, ctx, runstart, ri, re));
882: } else {
883: PetscCall(PetscBinaryInsertionSortWithArray_Private(arr, atarr, asize, barr, btarr, bsize, cmp, ctx, runstart, ri, re));
884: }
885: *runend = re;
886: } else *runend = ri;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: /*@
891: PetscTimSort - Sorts an array in place in increasing order using Tim Peters <https://bugs.python.org/file4451/timsort.txt> adaptive sorting algorithm.
893: Not Collective, No Fortran Support
895: Input Parameters:
896: + n - number of values
897: . arr - array to be sorted
898: . size - size in bytes of the datatype held in arr
899: . cmp - function pointer to comparison function
900: - ctx - optional context to be passed to comparison function, NULL if not needed
902: Output Parameter:
903: . arr - sorted array
905: Level: developer
907: Notes:
908: Timsort makes the assumption that input data is already likely partially ordered, or that it contains contiguous
909: sections (termed 'runs') where the data is locally ordered (but not necessarily globally ordered). It therefore aims to
910: select slices of the array in such a way that resulting mergesorts operate on near perfectly length-balanced arrays. To
911: do so it repeatedly triggers attempts throughout to merge adjacent runs.
913: Should one run continuously "win" a comparison the algorithm begins the "gallop" phase. It will aggressively search
914: the "winner" for the location of the "losers" next entry (and vice versa) to copy all preceding elements into place in
915: bulk. However if the data is truly unordered (as is the case with random data) the immense gains possible from these
916: searches are expected __not__ to repay their costs. While adjacent arrays are almost all nearly the same size, they
917: likely all contain similar data.
919: Example Usage:
920: The comparison function must follow the `qsort()` comparison function paradigm, returning the sign of the difference
921: between its arguments. If left < right : return -1, if left == right : return 0, if left > right : return 1. The user
922: may also
923: change or reverse the order of the sort by flipping the above. Note that stability of the sort is only guaranteed if
924: the comparison function forms a valid trigraph. For example when sorting an array of type "my_type" in increasing
925: order
927: .vb
928: int increasing_comparison_function(const void *left, const void *right, PetscCtx ctx) {
929: my_type l = *(my_type *) left, r = *(my_type *) right;
930: return (l < r) ? -1 : (l > r);
931: }
932: .ve
934: Note the context is unused here but you may use it to pass and subsequently access whatever information required
935: inside the comparison function. The context pointer will unaltered except for any changes made inside the comparison function.
936: Then pass the function
937: .vb
938: PetscTimSort(n, arr, sizeof(arr[0]), increasing_comparison_function, ctx)
939: .ve
941: Fortran Notes:
942: To use this from Fortran you must write a comparison subroutine with 4 arguments which accepts `left`, `right`, `ctx` and
943: returns `result`. For example
944: .vb
945: subroutine CompareIntegers(left,right,ctx,result)
946: implicit none
948: PetscInt,intent(in) :: left, right
949: type(UserCtx) :: ctx
950: integer,intent(out) :: result
952: if (left < right) then
953: result = -1
954: else if (left == right) then
955: result = 0
956: else
957: result = 1
958: end if
959: return
960: end subroutine CompareIntegers
961: .ve
963: .seealso: `PetscTimSortWithArray()`, `PetscIntSortSemiOrdered()`, `PetscRealSortSemiOrdered()`, `PetscMPIIntSortSemiOrdered()`
964: @*/
965: PetscErrorCode PetscTimSort(PetscInt n, void *arr, size_t size, int (*cmp)(const void *, const void *, void *), PetscCtx ctx)
966: {
967: PetscInt stacksize = 0, minrun, runstart = 0, runend = 0;
968: PetscTimSortStack runstack[128];
969: PetscTimSortBuffer buff;
970: /* stacksize = log_phi(n) = log_2(n)/log_2(phi), so 128 is enough for ~5.614e26 elements.
971: It is so unlikely that this limit is reached that this is __never__ checked for */
973: PetscFunctionBegin;
974: /* Compute minrun. Minrun should be (32, 65) such that N/minrun
975: is a power of 2 or one plus a power of 2 */
976: {
977: PetscInt t = n, r = 0;
978: /* r becomes 1 if the least significant bits contain at least one off bit */
979: while (t >= 64) {
980: r |= t & 1;
981: t >>= 1;
982: }
983: minrun = t + r;
984: }
985: if (PetscDefined(USE_DEBUG)) {
986: PetscCall(PetscInfo(NULL, "minrun = %" PetscInt_FMT "\n", minrun));
987: if (n < 64) PetscCall(PetscInfo(NULL, "n %" PetscInt_FMT " < 64, consider using PetscSortInt() instead\n", n));
988: else PetscCheck(!(minrun < 32) && !(minrun > 65), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Calculated minrun %" PetscInt_FMT " not in range (32,65)", minrun);
989: }
990: PetscCall(PetscMalloc1((size_t)minrun * size, &buff.ptr));
991: buff.size = (size_t)minrun * size;
992: buff.maxsize = (size_t)n * size;
993: MIN_GALLOP_GLOBAL = MIN_GALLOP_CONST_GLOBAL;
994: while (runstart < n) {
995: /* Check if additional entries are at least partially ordered and build natural run */
996: PetscCall(PetscTimSortBuildRun_Private((char *)arr, buff.ptr, size, cmp, ctx, n, minrun, runstart, &runend));
997: runstack[stacksize].start = runstart;
998: runstack[stacksize].size = runend - runstart + 1;
999: PetscCall(PetscTimSortMergeCollapse_Private((char *)arr, size, cmp, ctx, &buff, runstack, &stacksize));
1000: ++stacksize;
1001: runstart = runend + 1;
1002: }
1003: /* Have been inside while, so discard last stacksize++ */
1004: --stacksize;
1005: PetscCall(PetscTimSortForceCollapse_Private((char *)arr, size, cmp, ctx, &buff, runstack, stacksize));
1006: PetscCall(PetscFree(buff.ptr));
1007: MIN_GALLOP_GLOBAL = MIN_GALLOP_CONST_GLOBAL;
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@
1012: PetscTimSortWithArray - Sorts an array in place in increasing order using Tim Peters <https://bugs.python.org/file4451/timsort.txt> adaptive sorting algorithm and
1013: reorders a second array to match the first. The arrays need not be the same type.
1015: Not Collective, No Fortran Support
1017: Input Parameters:
1018: + n - number of values
1019: . asize - size in bytes of the datatype held in arr
1020: . bsize - size in bytes of the datatype held in barr
1021: . cmp - function pointer to comparison function
1022: - ctx - optional context to be passed to comparison function, NULL if not needed
1024: Input/Output Parameters:
1025: + arr - array to be sorted, on output it is sorted
1026: - barr - array to be reordered, on output it is reordered
1028: Level: developer
1030: Notes:
1031: The arrays need not be of the same type, however barr MUST contain at least as many elements as arr and the two CANNOT
1032: overlap.
1034: Timsort makes the assumption that input data is already likely partially ordered, or that it contains contiguous
1035: sections (termed 'runs') where the data is locally ordered (but not necessarily globally ordered). It therefore aims
1036: to select slices of the array in such a way that resulting mergesorts operate on near perfectly length-balanced
1037: arrays. To do so it repeatedly triggers attempts throughout to merge adjacent runs.
1039: Should one run continuously "win" a comparison the algorithm begins the "gallop" phase. It will aggressively search
1040: the "winner" for the location of the "losers" next entry (and vice versa) to copy all preceding elements into place in
1041: bulk. However if the data is truly unordered (as is the case with random data) the immense gains possible from these
1042: searches are expected __not__ to repay their costs. While adjacent arrays are almost all nearly the same size, they
1043: likely all contain similar data.
1045: Example Usage:
1046: The comparison function must follow the `qsort()` comparison function paradigm, returning the sign of the difference
1047: between its arguments. If left < right \: return -1, if left == right \: return 0, if left > right \: return 1. The user
1048: may also change or reverse the order of the sort by flipping the above. Note that stability of the sort is only
1049: guaranteed if the comparison function forms a valid trigraph. For example when sorting an array of type "my_type" in
1050: increasing order
1052: .vb
1053: int increasing_comparison_function(const void *left, const void *right, PetscCtx ctx) {
1054: my_type l = *(my_type *) left, r = *(my_type *) right;
1055: return (l < r) ? -1 : (l > r);
1056: }
1057: .ve
1059: Note the context is unused here but you may use it to pass and subsequently access whatever information required
1060: inside the comparison function. The context pointer will unaltered except for any changes made inside the comparison function.
1061: Then pass the function
1062: .vb
1063: PetscTimSortWithArray(n, arr, sizeof(arr[0]), barr, sizeof(barr[0]), increasing_comparison_function, ctx)
1064: .ve
1066: Fortran Notes:
1067: To use this from Fortran you must write a comparison subroutine with 4 arguments which accepts `left`, `right`, `ctx` and
1068: returns `result`. For example
1069: .vb
1070: subroutine CompareIntegers(left,right,ctx,result)
1071: implicit none
1073: PetscInt,intent(in) :: left, right
1074: type(UserCtx) :: ctx
1075: integer,intent(out) :: result
1077: if (left < right) then
1078: result = -1
1079: else if (left == right) then
1080: result = 0
1081: else
1082: result = 1
1083: end if
1084: return
1085: end subroutine CompareIntegers
1086: .ve
1088: .seealso: `PetscTimSort()`, `PetscIntSortSemiOrderedWithArray()`, `PetscRealSortSemiOrderedWithArrayInt()`, `PetscMPIIntSortSemiOrderedWithArray()`
1089: @*/
1090: PetscErrorCode PetscTimSortWithArray(PetscInt n, void *arr, size_t asize, void *barr, size_t bsize, int (*cmp)(const void *, const void *, void *), PetscCtx ctx)
1091: {
1092: PetscInt stacksize = 0, minrun, runstart = 0, runend = 0;
1093: PetscTimSortStack runstack[128];
1094: PetscTimSortBuffer abuff, bbuff;
1095: /* stacksize = log_phi(n) = log_2(n)/log_2(phi), so 128 is enough for ~5.614e26 elements.
1096: It is so unlikely that this limit is reached that this is __never__ checked for */
1098: PetscFunctionBegin;
1099: /* Compute minrun. Minrun should be (32, 65) such that N/minrun
1100: is a power of 2 or one plus a power of 2 */
1101: {
1102: PetscInt t = n, r = 0;
1103: /* r becomes 1 if the least significant bits contain at least one off bit */
1104: while (t >= 64) {
1105: r |= t & 1;
1106: t >>= 1;
1107: }
1108: minrun = t + r;
1109: }
1110: if (PetscDefined(USE_DEBUG)) {
1111: PetscCall(PetscInfo(NULL, "minrun = %" PetscInt_FMT "\n", minrun));
1112: PetscCheck(n < 64 || (minrun >= 32 && minrun <= 65), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Calculated minrun %" PetscInt_FMT " not in range (32,65)", minrun);
1113: }
1114: PetscCall(PetscMalloc1((size_t)minrun * asize, &abuff.ptr));
1115: abuff.size = (size_t)minrun * asize;
1116: abuff.maxsize = (size_t)n * asize;
1117: PetscCall(PetscMalloc1((size_t)minrun * bsize, &bbuff.ptr));
1118: bbuff.size = (size_t)minrun * bsize;
1119: bbuff.maxsize = (size_t)n * bsize;
1120: MIN_GALLOP_GLOBAL = MIN_GALLOP_CONST_GLOBAL;
1121: while (runstart < n) {
1122: /* Check if additional entries are at least partially ordered and build natural run */
1123: PetscCall(PetscTimSortBuildRunWithArray_Private((char *)arr, abuff.ptr, asize, (char *)barr, bbuff.ptr, bsize, cmp, ctx, n, minrun, runstart, &runend));
1124: runstack[stacksize].start = runstart;
1125: runstack[stacksize].size = runend - runstart + 1;
1126: PetscCall(PetscTimSortMergeCollapseWithArray_Private((char *)arr, asize, (char *)barr, bsize, cmp, ctx, &abuff, &bbuff, runstack, &stacksize));
1127: ++stacksize;
1128: runstart = runend + 1;
1129: }
1130: /* Have been inside while, so discard last stacksize++ */
1131: --stacksize;
1132: PetscCall(PetscTimSortForceCollapseWithArray_Private((char *)arr, asize, (char *)barr, bsize, cmp, ctx, &abuff, &bbuff, runstack, stacksize));
1133: PetscCall(PetscFree(abuff.ptr));
1134: PetscCall(PetscFree(bbuff.ptr));
1135: MIN_GALLOP_GLOBAL = MIN_GALLOP_CONST_GLOBAL;
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: /*@
1140: PetscIntSortSemiOrdered - Sorts an array of `PetscInt` in place in increasing order.
1142: Not Collective
1144: Input Parameters:
1145: + n - number of values
1146: - arr - array of integers
1148: Output Parameter:
1149: . arr - sorted array of integers
1151: Level: intermediate
1153: Notes:
1154: If the array is less than 64 entries long `PetscSortInt()` is automatically used.
1156: This function serves as an alternative to `PetscSortInt()`. While this function works for any array of integers it is
1157: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1158: recommended that the user benchmark their code to see which routine is fastest.
1160: .seealso: `PetscTimSort()`, `PetscSortInt()`, `PetscSortIntWithPermutation()`
1161: @*/
1162: PetscErrorCode PetscIntSortSemiOrdered(PetscInt n, PetscInt arr[])
1163: {
1164: PetscFunctionBegin;
1165: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1166: PetscAssertPointer(arr, 2);
1167: if (n < 64) {
1168: PetscCall(PetscSortInt(n, arr));
1169: } else {
1170: PetscCall(PetscTimSort(n, arr, sizeof(PetscInt), Compare_PetscInt_Private, NULL));
1171: }
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: /*@
1176: PetscIntSortSemiOrderedWithArray - Sorts an array of `PetscInt` in place in increasing order and reorders a second
1177: `PetscInt` array to match the first.
1179: Not Collective
1181: Input Parameter:
1182: . n - number of values
1184: Input/Output Parameters:
1185: + arr1 - array of integers to be sorted, modified on output
1186: - arr2 - array of integers to be reordered, modified on output
1188: Level: intermediate
1190: Notes:
1191: The arrays CANNOT overlap.
1193: This function serves as an alternative to `PetscSortIntWithArray()`. While this function works for any array of integers it is
1194: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1195: recommended that the user benchmark their code to see which routine is fastest.
1197: .seealso: `PetscTimSortWithArray()`, `PetscSortIntWithArray()`, `PetscSortIntWithPermutation()`
1198: @*/
1199: PetscErrorCode PetscIntSortSemiOrderedWithArray(PetscInt n, PetscInt arr1[], PetscInt arr2[])
1200: {
1201: PetscFunctionBegin;
1202: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1203: PetscAssertPointer(arr1, 2);
1204: PetscAssertPointer(arr2, 3);
1205: /* cannot export out to PetscIntSortWithArray here since it isn't stable */
1206: PetscCall(PetscTimSortWithArray(n, arr1, sizeof(PetscInt), arr2, sizeof(PetscInt), Compare_PetscInt_Private, NULL));
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@
1211: PetscMPIIntSortSemiOrdered - Sorts an array of `PetscMPIInt` in place in increasing order.
1213: Not Collective
1215: Input Parameters:
1216: + n - number of values
1217: - arr - array of `PetscMPIInt`
1219: Output Parameter:
1220: . arr - sorted array of integers
1222: Level: intermediate
1224: Notes:
1225: If the array is less than 64 entries long `PetscSortMPIInt()` is automatically used.
1227: This function serves as an alternative to `PetscSortMPIInt()`. While this function works for any array of `PetscMPIInt` it is
1228: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1229: recommended that the user benchmark their code to see which routine is fastest.
1231: .seealso: `PetscTimSort()`, `PetscSortMPIInt()`
1232: @*/
1233: PetscErrorCode PetscMPIIntSortSemiOrdered(PetscInt n, PetscMPIInt arr[])
1234: {
1235: PetscFunctionBegin;
1236: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1237: PetscAssertPointer(arr, 2);
1238: if (n < 64) {
1239: PetscCall(PetscSortMPIInt(n, arr));
1240: } else {
1241: PetscCall(PetscTimSort(n, arr, sizeof(PetscMPIInt), Compare_PetscMPIInt_Private, NULL));
1242: }
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: /*@
1247: PetscMPIIntSortSemiOrderedWithArray - Sorts an array of `PetscMPIInt` in place in increasing order and reorders a second `PetscMPIInt`
1248: array to match the first.
1250: Not Collective
1252: Input Parameter:
1253: . n - number of values
1255: Input/Output Parameters:
1256: + arr1 - array of integers to be sorted, modified on output
1257: - arr2 - array of integers to be reordered, modified on output
1259: Level: intermediate
1261: Notes:
1262: The arrays CANNOT overlap.
1264: This function serves as an alternative to `PetscSortMPIIntWithArray()`. While this function works for any array of integers it is
1265: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1266: recommended that the user benchmark their code to see which routine is fastest.
1268: .seealso: `PetscTimSortWithArray()`, `PetscSortMPIIntWithArray()`, `PetscSortIntWithPermutation()`
1269: @*/
1270: PetscErrorCode PetscMPIIntSortSemiOrderedWithArray(PetscInt n, PetscMPIInt arr1[], PetscMPIInt arr2[])
1271: {
1272: PetscFunctionBegin;
1273: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1274: PetscAssertPointer(arr1, 2);
1275: PetscAssertPointer(arr2, 3);
1276: /* cannot export out to PetscMPIIntSortWithArray here since it isn't stable */
1277: PetscCall(PetscTimSortWithArray(n, arr1, sizeof(PetscMPIInt), arr2, sizeof(PetscMPIInt), Compare_PetscMPIInt_Private, NULL));
1278: PetscFunctionReturn(PETSC_SUCCESS);
1279: }
1281: /*@
1282: PetscRealSortSemiOrdered - Sorts an array of `PetscReal` in place in increasing order.
1284: Not Collective
1286: Input Parameters:
1287: + n - number of values
1288: - arr - array of `PetscReal`
1290: Output Parameter:
1291: . arr - sorted array of integers
1293: Level: intermediate
1295: Notes:
1296: If the array is less than 64 entries long `PetscSortReal()` is automatically used.
1298: This function serves as an alternative to `PetscSortReal()`. While this function works for any array of `PetscReal` it is
1299: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1300: recommended that the user benchmark their code to see which routine is fastest.
1302: .seealso: `PetscTimSort()`, `PetscSortReal()`, `PetscSortRealWithPermutation()`
1303: @*/
1304: PetscErrorCode PetscRealSortSemiOrdered(PetscInt n, PetscReal arr[])
1305: {
1306: PetscFunctionBegin;
1307: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1308: PetscAssertPointer(arr, 2);
1309: if (n < 64) {
1310: PetscCall(PetscSortReal(n, arr));
1311: } else {
1312: PetscCall(PetscTimSort(n, arr, sizeof(PetscReal), Compare_PetscReal_Private, NULL));
1313: }
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }
1317: /*@
1318: PetscRealSortSemiOrderedWithArrayInt - Sorts an array of `PetscReal` in place in increasing order and reorders a second
1319: array of `PetscInt` to match the first.
1321: Not Collective
1323: Input Parameter:
1324: . n - number of values
1326: Input/Output Parameters:
1327: + arr1 - array of `PetscReal` to be sorted, modified on output
1328: - arr2 - array of `PetscInt` to be reordered, modified on output
1330: Level: intermediate
1332: Notes:
1333: This function serves as an alternative to `PetscSortRealWithArray()`. While this function works for any array of `PetscReal` it is
1334: significantly faster if the array is not totally random. There are exceptions to this and so it is __highly__
1335: recommended that the user benchmark their code to see which routine is fastest.
1337: .seealso: `PetscTimSortWithArray()`, `PetscSortRealWithArrayInt()`, `PetscSortRealWithPermutation()`
1338: @*/
1339: PetscErrorCode PetscRealSortSemiOrderedWithArrayInt(PetscInt n, PetscReal arr1[], PetscInt arr2[])
1340: {
1341: PetscFunctionBegin;
1342: if (n <= 1) PetscFunctionReturn(PETSC_SUCCESS);
1343: PetscAssertPointer(arr1, 2);
1344: PetscAssertPointer(arr2, 3);
1345: /* cannot export out to PetscRealSortWithArrayInt here since it isn't stable */
1346: PetscCall(PetscTimSortWithArray(n, arr1, sizeof(PetscReal), arr2, sizeof(PetscInt), Compare_PetscReal_Private, NULL));
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }