Actual source code: sortd.c


  2: /*
  3:    This file contains routines for sorting doubles.  Values are sorted in place.
  4:    These are provided because the general sort routines incur a great deal
  5:    of overhead in calling the comparison routines.

  7:  */
  8: #include <petscsys.h>
  9: #include <petsc/private/petscimpl.h>

 11: #define SWAP(a, b, t) \
 12:   { \
 13:     t = a; \
 14:     a = b; \
 15:     b = t; \
 16:   }

 18: /*@
 19:    PetscSortedReal - Determines whether the array of `PetscReal` is sorted.

 21:    Not Collective

 23:    Input Parameters:
 24: +  n  - number of values
 25: -  X  - array of integers

 27:    Output Parameter:
 28: .  sorted - flag whether the array is sorted

 30:    Level: intermediate

 32: .seealso: `PetscSortReal()`, `PetscSortedInt()`, `PetscSortedMPIInt()`
 33: @*/
 34: PetscErrorCode PetscSortedReal(PetscInt n, const PetscReal X[], PetscBool *sorted)
 35: {
 36:   PetscFunctionBegin;
 37:   PetscSorted(n, X, *sorted);
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
 42: static PetscErrorCode PetscSortReal_Private(PetscReal *v, PetscInt right)
 43: {
 44:   PetscInt  i, last;
 45:   PetscReal vl, tmp;

 47:   PetscFunctionBegin;
 48:   if (right <= 1) {
 49:     if (right == 1) {
 50:       if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
 51:     }
 52:     PetscFunctionReturn(PETSC_SUCCESS);
 53:   }
 54:   SWAP(v[0], v[right / 2], tmp);
 55:   vl   = v[0];
 56:   last = 0;
 57:   for (i = 1; i <= right; i++) {
 58:     if (v[i] < vl) {
 59:       last++;
 60:       SWAP(v[last], v[i], tmp);
 61:     }
 62:   }
 63:   SWAP(v[0], v[last], tmp);
 64:   PetscCall(PetscSortReal_Private(v, last - 1));
 65:   PetscCall(PetscSortReal_Private(v + last + 1, right - (last + 1)));
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: /*@
 70:    PetscSortReal - Sorts an array of `PetscReal` in place in increasing order.

 72:    Not Collective

 74:    Input Parameters:
 75: +  n  - number of values
 76: -  v  - array of doubles

 78:    Level: intermediate

 80:    Note:
 81:    This function serves as an alternative to `PetscRealSortSemiOrdered()`, and may perform faster especially if the array
 82:    is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
 83:    code to see which routine is fastest.

 85: .seealso: `PetscRealSortSemiOrdered()`, `PetscSortInt()`, `PetscSortRealWithPermutation()`, `PetscSortRealWithArrayInt()`
 86: @*/
 87: PetscErrorCode PetscSortReal(PetscInt n, PetscReal v[])
 88: {
 89:   PetscInt  j, k;
 90:   PetscReal tmp, vk;

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

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

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

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

149:    Not Collective

151:    Input Parameters:
152: +  n  - number of values
153: .  i  - array of integers
154: -  I - second array of integers

156:    Level: intermediate

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

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

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

187:    Not Collective

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

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

198:    Level: intermediate

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

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

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

226:    Not Collective

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

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

235:    Level: intermediate

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

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

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

258:    Not Collective

260:    Input Parameters:
261: +  ncut  - splitig index
262: -  n     - number of values to sort

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

270:    Level: intermediate

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

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

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

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

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

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

323:    Not Collective

325:    Input Parameters:
326: +  ncut  - splitig index
327: -  n     - number of values to sort

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

335:    Level: intermediate

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

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

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

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

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