1
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.
6
7 */
8 #include <petscsys.h> /*I "petscsys.h" I*/
9 #include <petsc/private/petscimpl.h>
10
11 #define SWAP(a,b,t) {t=a;a=b;b=t;}
12
13 /*@
14 PetscSortedReal - Determines whether the array is sorted.
15
16 Not Collective
17
18 Input Parameters:
19 + n - number of values
20 - X - array of integers
21
22 Output Parameters:
23 . sorted - flag whether the array is sorted
24
25 Level: intermediate
26
27 .seealso: PetscSortReal(), PetscSortedInt(), PetscSortedMPIInt()
28 @*/
PetscSortedReal(PetscInt n,const PetscReal X[],PetscBool * sorted)29 PetscErrorCode PetscSortedReal(PetscInt n,const PetscReal X[],PetscBool *sorted)
30 {
31 PetscFunctionBegin;
32 PetscSorted(n,X,*sorted);
33 PetscFunctionReturn(0);
34 }
35
36 /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
PetscSortReal_Private(PetscReal * v,PetscInt right)37 static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
38 {
39 PetscInt i,last;
40 PetscReal vl,tmp;
41
42 PetscFunctionBegin;
43 if (right <= 1) {
44 if (right == 1) {
45 if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
46 }
47 PetscFunctionReturn(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 PetscFunctionReturn(0);
59 }
60
61 /*@
62 PetscSortReal - Sorts an array of doubles in place in increasing order.
63
64 Not Collective
65
66 Input Parameters:
67 + n - number of values
68 - v - array of doubles
69
70 Notes:
71 This function serves as an alternative to PetscRealSortSemiOrdered(), and may perform faster especially if the array
72 is completely random. There are exceptions to this and so it is __highly__ recomended that the user benchmark their
73 code to see which routine is fastest.
74
75 Level: intermediate
76
77 .seealso: PetscRealSortSemiOrdered(), PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
78 @*/
PetscSortReal(PetscInt n,PetscReal v[])79 PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[])
80 {
81 PetscInt j,k;
82 PetscReal tmp,vk;
83
84 PetscFunctionBegin;
85 PetscValidPointer(v,2);
86 if (n<8) {
87 for (k=0; k<n; k++) {
88 vk = v[k];
89 for (j=k+1; j<n; j++) {
90 if (vk > v[j]) {
91 SWAP(v[k],v[j],tmp);
92 vk = v[k];
93 }
94 }
95 }
96 } else PetscSortReal_Private(v,n-1);
97 PetscFunctionReturn(0);
98 }
99
100 #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}
101
102 /* modified from PetscSortIntWithArray_Private */
PetscSortRealWithArrayInt_Private(PetscReal * v,PetscInt * V,PetscInt right)103 static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
104 {
105 PetscErrorCode ierr;
106 PetscInt i,last,itmp;
107 PetscReal rvl,rtmp;
108
109 PetscFunctionBegin;
110 if (right <= 1) {
111 if (right == 1) {
112 if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
113 }
114 PetscFunctionReturn(0);
115 }
116 SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
117 rvl = v[0];
118 last = 0;
119 for (i=1; i<=right; i++) {
120 if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
121 }
122 SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
123 ierr = PetscSortRealWithArrayInt_Private(v,V,last-1);CHKERRQ(ierr);
124 ierr = PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr);
125 PetscFunctionReturn(0);
126 }
127 /*@
128 PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
129 changes a second PetscInt array to match the sorted first array.
130
131 Not Collective
132
133 Input Parameters:
134 + n - number of values
135 . i - array of integers
136 - I - second array of integers
137
138 Level: intermediate
139
140 .seealso: PetscSortReal()
141 @*/
PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])142 PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
143 {
144 PetscErrorCode ierr;
145 PetscInt j,k,itmp;
146 PetscReal rk,rtmp;
147
148 PetscFunctionBegin;
149 PetscValidPointer(r,2);
150 PetscValidPointer(Ii,3);
151 if (n<8) {
152 for (k=0; k<n; k++) {
153 rk = r[k];
154 for (j=k+1; j<n; j++) {
155 if (rk > r[j]) {
156 SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
157 rk = r[k];
158 }
159 }
160 }
161 } else {
162 ierr = PetscSortRealWithArrayInt_Private(r,Ii,n-1);CHKERRQ(ierr);
163 }
164 PetscFunctionReturn(0);
165 }
166
167 /*@
168 PetscFindReal - Finds a PetscReal in a sorted array of PetscReals
169
170 Not Collective
171
172 Input Parameters:
173 + key - the value to locate
174 . n - number of values in the array
175 . ii - array of values
176 - eps - tolerance used to compare
177
178 Output Parameter:
179 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go
180
181 Level: intermediate
182
183 .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
184 @*/
PetscFindReal(PetscReal key,PetscInt n,const PetscReal t[],PetscReal eps,PetscInt * loc)185 PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
186 {
187 PetscInt lo = 0,hi = n;
188
189 PetscFunctionBegin;
190 PetscValidPointer(loc,4);
191 if (!n) {*loc = -1; PetscFunctionReturn(0);}
192 PetscValidPointer(t,3);
193 PetscCheckSorted(n,t);
194 while (hi - lo > 1) {
195 PetscInt mid = lo + (hi - lo)/2;
196 if (key < t[mid]) hi = mid;
197 else lo = mid;
198 }
199 *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
200 PetscFunctionReturn(0);
201 }
202
203 /*@
204 PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries
205
206 Not Collective
207
208 Input Parameters:
209 + n - number of values
210 - v - array of doubles
211
212 Output Parameter:
213 . n - number of non-redundant values
214
215 Level: intermediate
216
217 .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
218 @*/
PetscSortRemoveDupsReal(PetscInt * n,PetscReal v[])219 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
220 {
221 PetscErrorCode ierr;
222 PetscInt i,s = 0,N = *n, b = 0;
223
224 PetscFunctionBegin;
225 ierr = PetscSortReal(N,v);CHKERRQ(ierr);
226 for (i=0; i<N-1; i++) {
227 if (v[b+s+1] != v[b]) {
228 v[b+1] = v[b+s+1]; b++;
229 } else s++;
230 }
231 *n = N - s;
232 PetscFunctionReturn(0);
233 }
234
235 /*@
236 PetscSortSplit - Quick-sort split of an array of PetscScalars in place.
237
238 Not Collective
239
240 Input Parameters:
241 + ncut - splitig index
242 . n - number of values to sort
243 . a - array of values
244 - idx - index for array a
245
246 Output Parameters:
247 + a - permuted array of values such that its elements satisfy:
248 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
249 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
250 - idx - permuted index of array a
251
252 Level: intermediate
253
254 .seealso: PetscSortInt(), PetscSortRealWithPermutation()
255 @*/
PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])256 PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
257 {
258 PetscInt i,mid,last,itmp,j,first;
259 PetscScalar d,tmp;
260 PetscReal abskey;
261
262 PetscFunctionBegin;
263 first = 0;
264 last = n-1;
265 if (ncut < first || ncut > last) PetscFunctionReturn(0);
266
267 while (1) {
268 mid = first;
269 d = a[mid];
270 abskey = PetscAbsScalar(d);
271 i = last;
272 for (j = first + 1; j <= i; ++j) {
273 d = a[j];
274 if (PetscAbsScalar(d) >= abskey) {
275 ++mid;
276 /* interchange */
277 tmp = a[mid]; itmp = idx[mid];
278 a[mid] = a[j]; idx[mid] = idx[j];
279 a[j] = tmp; idx[j] = itmp;
280 }
281 }
282
283 /* interchange */
284 tmp = a[mid]; itmp = idx[mid];
285 a[mid] = a[first]; idx[mid] = idx[first];
286 a[first] = tmp; idx[first] = itmp;
287
288 /* test for while loop */
289 if (mid == ncut) break;
290 else if (mid > ncut) last = mid - 1;
291 else first = mid + 1;
292 }
293 PetscFunctionReturn(0);
294 }
295
296 /*@
297 PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.
298
299 Not Collective
300
301 Input Parameters:
302 + ncut - splitig index
303 . n - number of values to sort
304 . a - array of values in PetscReal
305 - idx - index for array a
306
307 Output Parameters:
308 + a - permuted array of real values such that its elements satisfy:
309 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
310 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
311 - idx - permuted index of array a
312
313 Level: intermediate
314
315 .seealso: PetscSortInt(), PetscSortRealWithPermutation()
316 @*/
PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])317 PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
318 {
319 PetscInt i,mid,last,itmp,j,first;
320 PetscReal d,tmp;
321 PetscReal abskey;
322
323 PetscFunctionBegin;
324 first = 0;
325 last = n-1;
326 if (ncut < first || ncut > last) PetscFunctionReturn(0);
327
328 while (1) {
329 mid = first;
330 d = a[mid];
331 abskey = PetscAbsReal(d);
332 i = last;
333 for (j = first + 1; j <= i; ++j) {
334 d = a[j];
335 if (PetscAbsReal(d) >= abskey) {
336 ++mid;
337 /* interchange */
338 tmp = a[mid]; itmp = idx[mid];
339 a[mid] = a[j]; idx[mid] = idx[j];
340 a[j] = tmp; idx[j] = itmp;
341 }
342 }
343
344 /* interchange */
345 tmp = a[mid]; itmp = idx[mid];
346 a[mid] = a[first]; idx[mid] = idx[first];
347 a[first] = tmp; idx[first] = itmp;
348
349 /* test for while loop */
350 if (mid == ncut) break;
351 else if (mid > ncut) last = mid - 1;
352 else first = mid + 1;
353 }
354 PetscFunctionReturn(0);
355 }
356