Actual source code: sortd.c
petsc-3.13.6 2020-09-29
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 comparision routines.
7: */
8: #include <petscsys.h>
9: #include <petsc/private/petscimpl.h>
11: #define SWAP(a,b,t) {t=a;a=b;b=t;}
13: /*@
14: PetscSortedReal - Determines whether the array is sorted.
16: Not Collective
18: Input Parameters:
19: + n - number of values
20: - X - array of integers
22: Output Parameters:
23: . sorted - flag whether the array is sorted
25: Level: intermediate
27: .seealso: PetscSortReal(), PetscSortedInt(), PetscSortedMPIInt()
28: @*/
29: PetscErrorCode PetscSortedReal(PetscInt n,const PetscReal X[],PetscBool *sorted)
30: {
32: PetscSorted(n,X,*sorted);
33: return(0);
34: }
36: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
37: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
38: {
39: PetscInt i,last;
40: PetscReal vl,tmp;
43: if (right <= 1) {
44: if (right == 1) {
45: if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
46: }
47: return(0);
48: }
49: SWAP(v[0],v[right/2],tmp);
50: vl = v[0];
51: last = 0;
52: for (i=1; i<=right; i++) {
53: if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
54: }
55: SWAP(v[0],v[last],tmp);
56: PetscSortReal_Private(v,last-1);
57: PetscSortReal_Private(v+last+1,right-(last+1));
58: return(0);
59: }
61: /*@
62: PetscSortReal - Sorts an array of doubles in place in increasing order.
64: Not Collective
66: Input Parameters:
67: + n - number of values
68: - v - array of doubles
70: Level: intermediate
72: .seealso: PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
73: @*/
74: PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
75: {
76: PetscInt j,k;
77: PetscReal tmp,vk;
81: if (n<8) {
82: for (k=0; k<n; k++) {
83: vk = v[k];
84: for (j=k+1; j<n; j++) {
85: if (vk > v[j]) {
86: SWAP(v[k],v[j],tmp);
87: vk = v[k];
88: }
89: }
90: }
91: } else PetscSortReal_Private(v,n-1);
92: return(0);
93: }
95: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
97: /* modified from PetscSortIntWithArray_Private */
98: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
99: {
101: PetscInt i,last,itmp;
102: PetscReal rvl,rtmp;
105: if (right <= 1) {
106: if (right == 1) {
107: if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
108: }
109: return(0);
110: }
111: SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
112: rvl = v[0];
113: last = 0;
114: for (i=1; i<=right; i++) {
115: if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
116: }
117: SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
118: PetscSortRealWithArrayInt_Private(v,V,last-1);
119: PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
120: return(0);
121: }
122: /*@
123: PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
124: changes a second PetscInt array to match the sorted first array.
126: Not Collective
128: Input Parameters:
129: + n - number of values
130: . i - array of integers
131: - I - second array of integers
133: Level: intermediate
135: .seealso: PetscSortReal()
136: @*/
137: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
138: {
140: PetscInt j,k,itmp;
141: PetscReal rk,rtmp;
146: if (n<8) {
147: for (k=0; k<n; k++) {
148: rk = r[k];
149: for (j=k+1; j<n; j++) {
150: if (rk > r[j]) {
151: SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
152: rk = r[k];
153: }
154: }
155: }
156: } else {
157: PetscSortRealWithArrayInt_Private(r,Ii,n-1);
158: }
159: return(0);
160: }
162: /*@
163: PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
165: Not Collective
167: Input Parameters:
168: + key - the value to locate
169: . n - number of values in the array
170: . ii - array of values
171: - eps - tolerance used to compare
173: Output Parameter:
174: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
176: Level: intermediate
178: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
179: @*/
180: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
181: {
182: PetscInt lo = 0,hi = n;
186: if (!n) {*loc = -1; return(0);}
189: while (hi - lo > 1) {
190: PetscInt mid = lo + (hi - lo)/2;
191: if (key < t[mid]) hi = mid;
192: else lo = mid;
193: }
194: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
195: return(0);
196: }
198: /*@
199: PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
201: Not Collective
203: Input Parameters:
204: + n - number of values
205: - v - array of doubles
207: Output Parameter:
208: . n - number of non-redundant values
210: Level: intermediate
212: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
213: @*/
214: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
215: {
217: PetscInt i,s = 0,N = *n, b = 0;
220: PetscSortReal(N,v);
221: for (i=0; i<N-1; i++) {
222: if (v[b+s+1] != v[b]) {
223: v[b+1] = v[b+s+1]; b++;
224: } else s++;
225: }
226: *n = N - s;
227: return(0);
228: }
230: /*@
231: PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
233: Not Collective
235: Input Parameters:
236: + ncut - splitig index
237: . n - number of values to sort
238: . a - array of values
239: - idx - index for array a
241: Output Parameters:
242: + a - permuted array of values such that its elements satisfy:
243: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
244: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
245: - idx - permuted index of array a
247: Level: intermediate
249: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
250: @*/
251: PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
252: {
253: PetscInt i,mid,last,itmp,j,first;
254: PetscScalar d,tmp;
255: PetscReal abskey;
258: first = 0;
259: last = n-1;
260: if (ncut < first || ncut > last) return(0);
262: while (1) {
263: mid = first;
264: d = a[mid];
265: abskey = PetscAbsScalar(d);
266: i = last;
267: for (j = first + 1; j <= i; ++j) {
268: d = a[j];
269: if (PetscAbsScalar(d) >= abskey) {
270: ++mid;
271: /* interchange */
272: tmp = a[mid]; itmp = idx[mid];
273: a[mid] = a[j]; idx[mid] = idx[j];
274: a[j] = tmp; idx[j] = itmp;
275: }
276: }
278: /* interchange */
279: tmp = a[mid]; itmp = idx[mid];
280: a[mid] = a[first]; idx[mid] = idx[first];
281: a[first] = tmp; idx[first] = itmp;
283: /* test for while loop */
284: if (mid == ncut) break;
285: else if (mid > ncut) last = mid - 1;
286: else first = mid + 1;
287: }
288: return(0);
289: }
291: /*@
292: PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
294: Not Collective
296: Input Parameters:
297: + ncut - splitig index
298: . n - number of values to sort
299: . a - array of values in PetscReal
300: - idx - index for array a
302: Output Parameters:
303: + a - permuted array of real values such that its elements satisfy:
304: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
305: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
306: - idx - permuted index of array a
308: Level: intermediate
310: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
311: @*/
312: PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
313: {
314: PetscInt i,mid,last,itmp,j,first;
315: PetscReal d,tmp;
316: PetscReal abskey;
319: first = 0;
320: last = n-1;
321: if (ncut < first || ncut > last) return(0);
323: while (1) {
324: mid = first;
325: d = a[mid];
326: abskey = PetscAbsReal(d);
327: i = last;
328: for (j = first + 1; j <= i; ++j) {
329: d = a[j];
330: if (PetscAbsReal(d) >= abskey) {
331: ++mid;
332: /* interchange */
333: tmp = a[mid]; itmp = idx[mid];
334: a[mid] = a[j]; idx[mid] = idx[j];
335: a[j] = tmp; idx[j] = itmp;
336: }
337: }
339: /* interchange */
340: tmp = a[mid]; itmp = idx[mid];
341: a[mid] = a[first]; idx[mid] = idx[first];
342: a[first] = tmp; idx[first] = itmp;
344: /* test for while loop */
345: if (mid == ncut) break;
346: else if (mid > ncut) last = mid - 1;
347: else first = mid + 1;
348: }
349: return(0);
350: }