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 Parameters:
 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:   PetscSorted(n, X, *sorted);
 37:   return 0;
 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:   if (right <= 1) {
 47:     if (right == 1) {
 48:       if (v[0] > v[1]) SWAP(v[0], v[1], tmp);
 49:     }
 50:     return 0;
 51:   }
 52:   SWAP(v[0], v[right / 2], tmp);
 53:   vl   = v[0];
 54:   last = 0;
 55:   for (i = 1; i <= right; i++) {
 56:     if (v[i] < vl) {
 57:       last++;
 58:       SWAP(v[last], v[i], tmp);
 59:     }
 60:   }
 61:   SWAP(v[0], v[last], tmp);
 62:   PetscSortReal_Private(v, last - 1);
 63:   PetscSortReal_Private(v + last + 1, right - (last + 1));
 64:   return 0;
 65: }

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

 70:    Not Collective

 72:    Input Parameters:
 73: +  n  - number of values
 74: -  v  - array of doubles

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

 81:    Level: intermediate

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

 91:   if (n < 8) {
 92:     for (k = 0; k < n; k++) {
 93:       vk = v[k];
 94:       for (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 PetscSortReal_Private(v, n - 1);
102:   return 0;
103: }

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

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

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

145:    Not Collective

147:    Input Parameters:
148: +  n  - number of values
149: .  i  - array of integers
150: -  I - second array of integers

152:    Level: intermediate

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

163:   if (n < 8) {
164:     for (k = 0; k < n; k++) {
165:       rk = r[k];
166:       for (j = k + 1; j < n; j++) {
167:         if (rk > r[j]) {
168:           SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp);
169:           rk = r[k];
170:         }
171:       }
172:     }
173:   } else {
174:     PetscSortRealWithArrayInt_Private(r, Ii, n - 1);
175:   }
176:   return 0;
177: }

179: /*@
180:   PetscFindReal - Finds a PetscReal` in a sorted array of `PetscReal`s

182:    Not Collective

184:    Input Parameters:
185: +  key - the value to locate
186: .  n   - number of values in the array
187: .  ii  - array of values
188: -  eps - tolerance used to compare

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

193:    Level: intermediate

195: .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()`
196: @*/
197: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
198: {
199:   PetscInt lo = 0, hi = n;

202:   if (!n) {
203:     *loc = -1;
204:     return 0;
205:   }
208:   while (hi - lo > 1) {
209:     PetscInt mid = lo + (hi - lo) / 2;
210:     if (key < t[mid]) hi = mid;
211:     else lo = mid;
212:   }
213:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
214:   return 0;
215: }

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

220:    Not Collective

222:    Input Parameters:
223: +  n  - number of values
224: -  v  - array of doubles

226:    Output Parameter:
227: .  n - number of non-redundant values

229:    Level: intermediate

231: .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()`
232: @*/
233: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[])
234: {
235:   PetscInt i, s = 0, N = *n, b = 0;

237:   PetscSortReal(N, v);
238:   for (i = 0; i < N - 1; i++) {
239:     if (v[b + s + 1] != v[b]) {
240:       v[b + 1] = v[b + s + 1];
241:       b++;
242:     } else s++;
243:   }
244:   *n = N - s;
245:   return 0;
246: }

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

251:    Not Collective

253:    Input Parameters:
254: +  ncut  - splitig index
255: -  n     - number of values to sort

257:    Input/Output Parameters:
258: +  a     - array of values, on output the values are permuted such that its elements satisfy:
259:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
260:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
261: -  idx   - index for array a, on output permuted accordingly

263:    Level: intermediate

265: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
266: @*/
267: PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[])
268: {
269:   PetscInt    i, mid, last, itmp, j, first;
270:   PetscScalar d, tmp;
271:   PetscReal   abskey;

273:   first = 0;
274:   last  = n - 1;
275:   if (ncut < first || ncut > last) return 0;

277:   while (1) {
278:     mid    = first;
279:     d      = a[mid];
280:     abskey = PetscAbsScalar(d);
281:     i      = last;
282:     for (j = first + 1; j <= i; ++j) {
283:       d = a[j];
284:       if (PetscAbsScalar(d) >= abskey) {
285:         ++mid;
286:         /* interchange */
287:         tmp      = a[mid];
288:         itmp     = idx[mid];
289:         a[mid]   = a[j];
290:         idx[mid] = idx[j];
291:         a[j]     = tmp;
292:         idx[j]   = itmp;
293:       }
294:     }

296:     /* interchange */
297:     tmp        = a[mid];
298:     itmp       = idx[mid];
299:     a[mid]     = a[first];
300:     idx[mid]   = idx[first];
301:     a[first]   = tmp;
302:     idx[first] = itmp;

304:     /* test for while loop */
305:     if (mid == ncut) break;
306:     else if (mid > ncut) last = mid - 1;
307:     else first = mid + 1;
308:   }
309:   return 0;
310: }

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

315:    Not Collective

317:    Input Parameters:
318: +  ncut  - splitig index
319: -  n     - number of values to sort

321:    Input/Output Parameters:
322: +  a     - array of values, on output the values are permuted such that its elements satisfy:
323:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
324:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
325: -  idx   - index for array a, on output permuted accordingly

327:    Level: intermediate

329: .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()`
330: @*/
331: PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[])
332: {
333:   PetscInt  i, mid, last, itmp, j, first;
334:   PetscReal d, tmp;
335:   PetscReal abskey;

337:   first = 0;
338:   last  = n - 1;
339:   if (ncut < first || ncut > last) return 0;

341:   while (1) {
342:     mid    = first;
343:     d      = a[mid];
344:     abskey = PetscAbsReal(d);
345:     i      = last;
346:     for (j = first + 1; j <= i; ++j) {
347:       d = a[j];
348:       if (PetscAbsReal(d) >= abskey) {
349:         ++mid;
350:         /* interchange */
351:         tmp      = a[mid];
352:         itmp     = idx[mid];
353:         a[mid]   = a[j];
354:         idx[mid] = idx[j];
355:         a[j]     = tmp;
356:         idx[j]   = itmp;
357:       }
358:     }

360:     /* interchange */
361:     tmp        = a[mid];
362:     itmp       = idx[mid];
363:     a[mid]     = a[first];
364:     idx[mid]   = idx[first];
365:     a[first]   = tmp;
366:     idx[first] = itmp;

368:     /* test for while loop */
369:     if (mid == ncut) break;
370:     else if (mid > ncut) last = mid - 1;
371:     else first = mid + 1;
372:   }
373:   return 0;
374: }