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 - initial number of values
235: - v - array of values

237:   Output Parameters:
238: + n - number of non-redundant values
239: - v - array of non-redundant values

241:   Level: intermediate

243: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
244: @*/
245: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
246: {
247:   PetscInt i, s = 0, N = *n, b = 0;

249:   PetscFunctionBegin;
250:   PetscCall(PetscSortReal(N, v));
251:   for (i = 0; i < N - 1; i++) {
252:     if (v[b + s + 1] != v[b]) {
253:       v[b + 1] = v[b + s + 1];
254:       b++;
255:     } else s++;
256:   }
257:   *n = N - s;
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: /*@
262:   PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place.

264:   Not Collective

266:   Input Parameters:
267: + ncut - splitting index
268: - n    - number of values to sort

270:   Input/Output Parameters:
271: + a   - array of values, on output the values are permuted such that its elements satisfy:
272:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
273:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
274: - idx - index for array a, on output permuted accordingly

276:   Level: intermediate

278: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
279: @*/
280: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
281: {
282:   PetscInt    i, mid, last, itmp, j, first;
283:   PetscScalar d, tmp;
284:   PetscReal   abskey;

286:   PetscFunctionBegin;
287:   first = 0;
288:   last  = n - 1;
289:   if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);

291:   while (1) {
292:     mid    = first;
293:     d      = a[mid];
294:     abskey = PetscAbsScalar(d);
295:     i      = last;
296:     for (j = first + 1; j <= i; ++j) {
297:       d = a[j];
298:       if (PetscAbsScalar(d) >= abskey) {
299:         ++mid;
300:         /* interchange */
301:         tmp      = a[mid];
302:         itmp     = idx[mid];
303:         a[mid]   = a[j];
304:         idx[mid] = idx[j];
305:         a[j]     = tmp;
306:         idx[j]   = itmp;
307:       }
308:     }

310:     /* interchange */
311:     tmp        = a[mid];
312:     itmp       = idx[mid];
313:     a[mid]     = a[first];
314:     idx[mid]   = idx[first];
315:     a[first]   = tmp;
316:     idx[first] = itmp;

318:     /* test for while loop */
319:     if (mid == ncut) break;
320:     else if (mid > ncut) last = mid - 1;
321:     else first = mid + 1;
322:   }
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: /*@
327:   PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place.

329:   Not Collective

331:   Input Parameters:
332: + ncut - splitting index
333: - n    - number of values to sort

335:   Input/Output Parameters:
336: + a   - array of values, on output the values are permuted such that its elements satisfy:
337:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
338:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
339: - idx - index for array a, on output permuted accordingly

341:   Level: intermediate

343: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
344: @*/
345: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
346: {
347:   PetscInt  i, mid, last, itmp, j, first;
348:   PetscReal d, tmp;
349:   PetscReal abskey;

351:   PetscFunctionBegin;
352:   first = 0;
353:   last  = n - 1;
354:   if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);

356:   while (1) {
357:     mid    = first;
358:     d      = a[mid];
359:     abskey = PetscAbsReal(d);
360:     i      = last;
361:     for (j = first + 1; j <= i; ++j) {
362:       d = a[j];
363:       if (PetscAbsReal(d) >= abskey) {
364:         ++mid;
365:         /* interchange */
366:         tmp      = a[mid];
367:         itmp     = idx[mid];
368:         a[mid]   = a[j];
369:         idx[mid] = idx[j];
370:         a[j]     = tmp;
371:         idx[j]   = itmp;
372:       }
373:     }

375:     /* interchange */
376:     tmp        = a[mid];
377:     itmp       = idx[mid];
378:     a[mid]     = a[first];
379:     idx[mid]   = idx[first];
380:     a[first]   = tmp;
381:     idx[first] = itmp;

383:     /* test for while loop */
384:     if (mid == ncut) break;
385:     else if (mid > ncut) last = mid - 1;
386:     else first = mid + 1;
387:   }
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }