Actual source code: sortd.c
petsc-3.12.5 2020-03-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: /* 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: .seealso: PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
50: @*/
51: PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
52: {
53: PetscInt j,k;
54: PetscReal tmp,vk;
58: if (n<8) {
59: for (k=0; k<n; k++) {
60: vk = v[k];
61: for (j=k+1; j<n; j++) {
62: if (vk > v[j]) {
63: SWAP(v[k],v[j],tmp);
64: vk = v[k];
65: }
66: }
67: }
68: } else PetscSortReal_Private(v,n-1);
69: return(0);
70: }
72: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
74: /* modified from PetscSortIntWithArray_Private */
75: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
76: {
78: PetscInt i,last,itmp;
79: PetscReal rvl,rtmp;
82: if (right <= 1) {
83: if (right == 1) {
84: if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
85: }
86: return(0);
87: }
88: SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
89: rvl = v[0];
90: last = 0;
91: for (i=1; i<=right; i++) {
92: if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
93: }
94: SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
95: PetscSortRealWithArrayInt_Private(v,V,last-1);
96: PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
97: return(0);
98: }
99: /*@
100: PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
101: changes a second PetscInt array to match the sorted first array.
103: Not Collective
105: Input Parameters:
106: + n - number of values
107: . i - array of integers
108: - I - second array of integers
110: Level: intermediate
112: .seealso: PetscSortReal()
113: @*/
114: PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
115: {
117: PetscInt j,k,itmp;
118: PetscReal rk,rtmp;
123: if (n<8) {
124: for (k=0; k<n; k++) {
125: rk = r[k];
126: for (j=k+1; j<n; j++) {
127: if (rk > r[j]) {
128: SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
129: rk = r[k];
130: }
131: }
132: }
133: } else {
134: PetscSortRealWithArrayInt_Private(r,Ii,n-1);
135: }
136: return(0);
137: }
139: /*@
140: PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
142: Not Collective
144: Input Parameters:
145: + key - the value to locate
146: . n - number of values in the array
147: . ii - array of values
148: - eps - tolerance used to compare
150: Output Parameter:
151: . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
153: Level: intermediate
155: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
156: @*/
157: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
158: {
159: PetscInt lo = 0,hi = n;
163: if (!n) {*loc = -1; return(0);}
165: while (hi - lo > 1) {
166: PetscInt mid = lo + (hi - lo)/2;
167: if (key < t[mid]) hi = mid;
168: else lo = mid;
169: }
170: *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
171: return(0);
172: }
174: /*@
175: PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
177: Not Collective
179: Input Parameters:
180: + n - number of values
181: - v - array of doubles
183: Output Parameter:
184: . n - number of non-redundant values
186: Level: intermediate
188: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
189: @*/
190: PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
191: {
193: PetscInt i,s = 0,N = *n, b = 0;
196: PetscSortReal(N,v);
197: for (i=0; i<N-1; i++) {
198: if (v[b+s+1] != v[b]) {
199: v[b+1] = v[b+s+1]; b++;
200: } else s++;
201: }
202: *n = N - s;
203: return(0);
204: }
206: /*@
207: PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
209: Not Collective
211: Input Parameters:
212: + ncut - splitig index
213: . n - number of values to sort
214: . a - array of values
215: - idx - index for array a
217: Output Parameters:
218: + a - permuted array of values such that its elements satisfy:
219: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
220: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
221: - idx - permuted index of array a
223: Level: intermediate
225: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
226: @*/
227: PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
228: {
229: PetscInt i,mid,last,itmp,j,first;
230: PetscScalar d,tmp;
231: PetscReal abskey;
234: first = 0;
235: last = n-1;
236: if (ncut < first || ncut > last) return(0);
238: while (1) {
239: mid = first;
240: d = a[mid];
241: abskey = PetscAbsScalar(d);
242: i = last;
243: for (j = first + 1; j <= i; ++j) {
244: d = a[j];
245: if (PetscAbsScalar(d) >= abskey) {
246: ++mid;
247: /* interchange */
248: tmp = a[mid]; itmp = idx[mid];
249: a[mid] = a[j]; idx[mid] = idx[j];
250: a[j] = tmp; idx[j] = itmp;
251: }
252: }
254: /* interchange */
255: tmp = a[mid]; itmp = idx[mid];
256: a[mid] = a[first]; idx[mid] = idx[first];
257: a[first] = tmp; idx[first] = itmp;
259: /* test for while loop */
260: if (mid == ncut) break;
261: else if (mid > ncut) last = mid - 1;
262: else first = mid + 1;
263: }
264: return(0);
265: }
267: /*@
268: PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
270: Not Collective
272: Input Parameters:
273: + ncut - splitig index
274: . n - number of values to sort
275: . a - array of values in PetscReal
276: - idx - index for array a
278: Output Parameters:
279: + a - permuted array of real values such that its elements satisfy:
280: abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
281: abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
282: - idx - permuted index of array a
284: Level: intermediate
286: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
287: @*/
288: PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
289: {
290: PetscInt i,mid,last,itmp,j,first;
291: PetscReal d,tmp;
292: PetscReal abskey;
295: first = 0;
296: last = n-1;
297: if (ncut < first || ncut > last) return(0);
299: while (1) {
300: mid = first;
301: d = a[mid];
302: abskey = PetscAbsReal(d);
303: i = last;
304: for (j = first + 1; j <= i; ++j) {
305: d = a[j];
306: if (PetscAbsReal(d) >= abskey) {
307: ++mid;
308: /* interchange */
309: tmp = a[mid]; itmp = idx[mid];
310: a[mid] = a[j]; idx[mid] = idx[j];
311: a[j] = tmp; idx[j] = itmp;
312: }
313: }
315: /* interchange */
316: tmp = a[mid]; itmp = idx[mid];
317: a[mid] = a[first]; idx[mid] = idx[first];
318: a[first] = tmp; idx[first] = itmp;
320: /* test for while loop */
321: if (mid == ncut) break;
322: else if (mid > ncut) last = mid - 1;
323: else first = mid + 1;
324: }
325: return(0);
326: }