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(PetscInt 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, PetscInt right)
 42: {
 43:   PetscInt  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(PetscInt n, PetscReal v[])
 87: {
 88:   PetscInt  j, k;
 89:   PetscReal tmp, vk;

 91:   PetscFunctionBegin;
 92:   PetscAssertPointer(v, 2);
 93:   if (n < 8) {
 94:     for (k = 0; k < n; k++) {
 95:       vk = v[k];
 96:       for (j = k + 1; j < n; j++) {
 97:         if (vk > v[j]) {
 98:           SWAP(v[k], v[j], tmp);
 99:           vk = v[k];
100:         }
101:       }
102:     }
103:   } else PetscCall(PetscSortReal_Private(v, n - 1));
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: #define SWAP2ri(a, b, c, d, rt, it) \
108:   do { \
109:     rt = a; \
110:     a  = b; \
111:     b  = rt; \
112:     it = c; \
113:     c  = d; \
114:     d  = it; \
115:   } while (0)

117: /* modified from PetscSortIntWithArray_Private */
118: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscInt right)
119: {
120:   PetscInt  i, last, itmp;
121:   PetscReal rvl, rtmp;

123:   PetscFunctionBegin;
124:   if (right <= 1) {
125:     if (right == 1) {
126:       if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp);
127:     }
128:     PetscFunctionReturn(PETSC_SUCCESS);
129:   }
130:   SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp);
131:   rvl  = v[0];
132:   last = 0;
133:   for (i = 1; i <= right; i++) {
134:     if (v[i] < rvl) {
135:       last++;
136:       SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp);
137:     }
138:   }
139:   SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp);
140:   PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1));
141:   PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1)));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }
144: /*@
145:   PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order;
146:   changes a second `PetscInt` array to match the sorted first array.

148:   Not Collective

150:   Input Parameters:
151: + n  - number of values
152: . Ii - array of integers
153: - r  - second array of integers

155:   Level: intermediate

157: .seealso: `PetscSortReal()`
158: @*/
159: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n, PetscReal r[], PetscInt Ii[])
160: {
161:   PetscInt  j, k, itmp;
162:   PetscReal rk, rtmp;

164:   PetscFunctionBegin;
165:   PetscAssertPointer(r, 2);
166:   PetscAssertPointer(Ii, 3);
167:   if (n < 8) {
168:     for (k = 0; k < n; k++) {
169:       rk = r[k];
170:       for (j = k + 1; j < n; j++) {
171:         if (rk > r[j]) {
172:           SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
173:           rk = r[k];
174:         }
175:       }
176:     }
177:   } else {
178:     PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1));
179:   }
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*@
184:   PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s

186:   Not Collective

188:   Input Parameters:
189: + key - the value to locate
190: . n   - number of values in the array
191: . t   - array of values
192: - eps - tolerance used to compare

194:   Output Parameter:
195: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go

197:   Level: intermediate

199: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
200: @*/
201: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
202: {
203:   PetscInt lo = 0, hi = n;

205:   PetscFunctionBegin;
206:   PetscAssertPointer(loc, 5);
207:   if (!n) {
208:     *loc = -1;
209:     PetscFunctionReturn(PETSC_SUCCESS);
210:   }
211:   PetscAssertPointer(t, 3);
212:   PetscCheckSorted(n, t);
213:   while (hi - lo > 1) {
214:     PetscInt mid = lo + (hi - lo) / 2;
215:     if (key < t[mid]) hi = mid;
216:     else lo = mid;
217:   }
218:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*@
223:   PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries

225:   Not Collective

227:   Input Parameters:
228: + n - number of values
229: - v - array of doubles

231:   Output Parameter:
232: . n - number of non-redundant values

234:   Level: intermediate

236: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
237: @*/
238: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
239: {
240:   PetscInt i, s = 0, N = *n, b = 0;

242:   PetscFunctionBegin;
243:   PetscCall(PetscSortReal(N, v));
244:   for (i = 0; i < N - 1; i++) {
245:     if (v[b + s + 1] != v[b]) {
246:       v[b + 1] = v[b + s + 1];
247:       b++;
248:     } else s++;
249:   }
250:   *n = N - s;
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

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

257:   Not Collective

259:   Input Parameters:
260: + ncut - splitting index
261: - n    - number of values to sort

263:   Input/Output Parameters:
264: + a   - array of values, on output the values are permuted such that its elements satisfy:
265:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
266:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
267: - idx - index for array a, on output permuted accordingly

269:   Level: intermediate

271: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
272: @*/
273: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
274: {
275:   PetscInt    i, mid, last, itmp, j, first;
276:   PetscScalar d, tmp;
277:   PetscReal   abskey;

279:   PetscFunctionBegin;
280:   first = 0;
281:   last  = n - 1;
282:   if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);

284:   while (1) {
285:     mid    = first;
286:     d      = a[mid];
287:     abskey = PetscAbsScalar(d);
288:     i      = last;
289:     for (j = first + 1; j <= i; ++j) {
290:       d = a[j];
291:       if (PetscAbsScalar(d) >= abskey) {
292:         ++mid;
293:         /* interchange */
294:         tmp      = a[mid];
295:         itmp     = idx[mid];
296:         a[mid]   = a[j];
297:         idx[mid] = idx[j];
298:         a[j]     = tmp;
299:         idx[j]   = itmp;
300:       }
301:     }

303:     /* interchange */
304:     tmp        = a[mid];
305:     itmp       = idx[mid];
306:     a[mid]     = a[first];
307:     idx[mid]   = idx[first];
308:     a[first]   = tmp;
309:     idx[first] = itmp;

311:     /* test for while loop */
312:     if (mid == ncut) break;
313:     else if (mid > ncut) last = mid - 1;
314:     else first = mid + 1;
315:   }
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

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

322:   Not Collective

324:   Input Parameters:
325: + ncut - splitting index
326: - n    - number of values to sort

328:   Input/Output Parameters:
329: + a   - array of values, on output the values are permuted such that its elements satisfy:
330:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
331:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
332: - idx - index for array a, on output permuted accordingly

334:   Level: intermediate

336: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
337: @*/
338: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
339: {
340:   PetscInt  i, mid, last, itmp, j, first;
341:   PetscReal d, tmp;
342:   PetscReal abskey;

344:   PetscFunctionBegin;
345:   first = 0;
346:   last  = n - 1;
347:   if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS);

349:   while (1) {
350:     mid    = first;
351:     d      = a[mid];
352:     abskey = PetscAbsReal(d);
353:     i      = last;
354:     for (j = first + 1; j <= i; ++j) {
355:       d = a[j];
356:       if (PetscAbsReal(d) >= abskey) {
357:         ++mid;
358:         /* interchange */
359:         tmp      = a[mid];
360:         itmp     = idx[mid];
361:         a[mid]   = a[j];
362:         idx[mid] = idx[j];
363:         a[j]     = tmp;
364:         idx[j]   = itmp;
365:       }
366:     }

368:     /* interchange */
369:     tmp        = a[mid];
370:     itmp       = idx[mid];
371:     a[mid]     = a[first];
372:     idx[mid]   = idx[first];
373:     a[first]   = tmp;
374:     idx[first] = itmp;

376:     /* test for while loop */
377:     if (mid == ncut) break;
378:     else if (mid > ncut) last = mid - 1;
379:     else first = mid + 1;
380:   }
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }