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: }