Actual source code: sortd.c
petsc-3.9.4 2018-09-11
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: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
14: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
15: {
16: PetscInt i,last;
17: PetscReal vl,tmp;
20: if (right <= 1) {
21: if (right == 1) {
22: if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
23: }
24: return(0);
25: }
26: SWAP(v[0],v[right/2],tmp);
27: vl = v[0];
28: last = 0;
29: for (i=1; i<=right; i++) {
30: if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
31: }
32: SWAP(v[0],v[last],tmp);
33: PetscSortReal_Private(v,last-1);
34: PetscSortReal_Private(v+last+1,right-(last+1));
35: return(0);
36: }
38: /*@
39: PetscSortReal - Sorts an array of doubles in place in increasing order.
41: Not Collective
43: Input Parameters:
44: + n - number of values
45: - v - array of doubles
47: Level: intermediate
49: Concepts: sorting^doubles
51: .seealso: PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
52: @*/
53: PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
54: {
55: PetscInt j,k;
56: PetscReal tmp,vk;
60: if (n<8) {
61: for (k=0; k<n; k++) {
62: vk = v[k];
63: for (j=k+1; j<n; j++) {
64: if (vk > v[j]) {
65: SWAP(v[k],v[j],tmp);
66: vk = v[k];
67: }
68: }
69: }
70: } else PetscSortReal_Private(v,n-1);
71: return(0);
72: }
74: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
76: /* modified from PetscSortIntWithArray_Private */
77: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
78: {
80: PetscInt i,last,itmp;
81: PetscReal rvl,rtmp;
84: if (right <= 1) {
85: if (right == 1) {
86: if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
87: }
88: return(0);
89: }
90: SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
91: rvl = v[0];
92: last = 0;
93: for (i=1; i<=right; i++) {
94: if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
95: }
96: SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
97: PetscSortRealWithArrayInt_Private(v,V,last-1);
98: PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
99: return(0);
100: }
101: /*@
102: PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
103: changes a second PetscInt array to match the sorted first array.
105: Not Collective
107: Input Parameters:
108: + n - number of values
109: . i - array of integers
110: - I - second array of integers
112: Level: intermediate
114: Concepts: sorting^ints with array
116: .seealso: PetscSortReal()
117: @*/
118: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
119: {
121: PetscInt j,k,itmp;
122: PetscReal rk,rtmp;
127: if (n<8) {
128: for (k=0; k<n; k++) {
129: rk = r[k];
130: for (j=k+1; j<n; j++) {
131: if (rk > r[j]) {
132: SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
133: rk = r[k];
134: }
135: }
136: }
137: } else {
138: PetscSortRealWithArrayInt_Private(r,Ii,n-1);
139: }
140: return(0);
141: }
143: /*@
144: PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
146: Not Collective
148: Input Parameters:
149: + key - the value to locate
150: . n - number of values in the array
151: . ii - array of values
152: - eps - tolerance used to compare
154: Output Parameter:
155: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
157: Level: intermediate
159: Concepts: sorting^ints
161: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
162: @*/
163: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
164: {
165: PetscInt lo = 0,hi = n;
169: if (!n) {*loc = -1; return(0);}
171: while (hi - lo > 1) {
172: PetscInt mid = lo + (hi - lo)/2;
173: if (key < t[mid]) hi = mid;
174: else lo = mid;
175: }
176: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
177: return(0);
178: }
180: /*@
181: PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
183: Not Collective
185: Input Parameters:
186: + n - number of values
187: - v - array of doubles
189: Output Parameter:
190: . n - number of non-redundant values
192: Level: intermediate
194: Concepts: sorting^doubles
196: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
197: @*/
198: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
199: {
201: PetscInt i,s = 0,N = *n, b = 0;
204: PetscSortReal(N,v);
205: for (i=0; i<N-1; i++) {
206: if (v[b+s+1] != v[b]) {
207: v[b+1] = v[b+s+1]; b++;
208: } else s++;
209: }
210: *n = N - s;
211: return(0);
212: }
214: /*@
215: PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
217: Not Collective
219: Input Parameters:
220: + ncut - splitig index
221: . n - number of values to sort
222: . a - array of values
223: - idx - index for array a
225: Output Parameters:
226: + a - permuted array of values such that its elements satisfy:
227: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
228: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
229: - idx - permuted index of array a
231: Level: intermediate
233: Concepts: sorting^doubles
235: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
236: @*/
237: PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
238: {
239: PetscInt i,mid,last,itmp,j,first;
240: PetscScalar d,tmp;
241: PetscReal abskey;
244: first = 0;
245: last = n-1;
246: if (ncut < first || ncut > last) return(0);
248: while (1) {
249: mid = first;
250: d = a[mid];
251: abskey = PetscAbsScalar(d);
252: i = last;
253: for (j = first + 1; j <= i; ++j) {
254: d = a[j];
255: if (PetscAbsScalar(d) >= abskey) {
256: ++mid;
257: /* interchange */
258: tmp = a[mid]; itmp = idx[mid];
259: a[mid] = a[j]; idx[mid] = idx[j];
260: a[j] = tmp; idx[j] = itmp;
261: }
262: }
264: /* interchange */
265: tmp = a[mid]; itmp = idx[mid];
266: a[mid] = a[first]; idx[mid] = idx[first];
267: a[first] = tmp; idx[first] = itmp;
269: /* test for while loop */
270: if (mid == ncut) break;
271: else if (mid > ncut) last = mid - 1;
272: else first = mid + 1;
273: }
274: return(0);
275: }
277: /*@
278: PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
280: Not Collective
282: Input Parameters:
283: + ncut - splitig index
284: . n - number of values to sort
285: . a - array of values in PetscReal
286: - idx - index for array a
288: Output Parameters:
289: + a - permuted array of real values such that its elements satisfy:
290: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
291: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
292: - idx - permuted index of array a
294: Level: intermediate
296: Concepts: sorting^doubles
298: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
299: @*/
300: PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
301: {
302: PetscInt i,mid,last,itmp,j,first;
303: PetscReal d,tmp;
304: PetscReal abskey;
307: first = 0;
308: last = n-1;
309: if (ncut < first || ncut > last) return(0);
311: while (1) {
312: mid = first;
313: d = a[mid];
314: abskey = PetscAbsReal(d);
315: i = last;
316: for (j = first + 1; j <= i; ++j) {
317: d = a[j];
318: if (PetscAbsReal(d) >= abskey) {
319: ++mid;
320: /* interchange */
321: tmp = a[mid]; itmp = idx[mid];
322: a[mid] = a[j]; idx[mid] = idx[j];
323: a[j] = tmp; idx[j] = itmp;
324: }
325: }
327: /* interchange */
328: tmp = a[mid]; itmp = idx[mid];
329: a[mid] = a[first]; idx[mid] = idx[first];
330: a[first] = tmp; idx[first] = itmp;
332: /* test for while loop */
333: if (mid == ncut) break;
334: else if (mid > ncut) last = mid - 1;
335: else first = mid + 1;
336: }
337: return(0);
338: }