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 comparison 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 @*/ 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 */ 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__ recommended 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 @*/ 79 PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[]) 80 { 81 PetscInt j,k; 82 PetscReal tmp,vk; 83 84 PetscFunctionBegin; 85 PetscValidRealPointer(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 */ 103 static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right) 104 { 105 PetscInt i,last,itmp; 106 PetscReal rvl,rtmp; 107 108 PetscFunctionBegin; 109 if (right <= 1) { 110 if (right == 1) { 111 if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp); 112 } 113 PetscFunctionReturn(0); 114 } 115 SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp); 116 rvl = v[0]; 117 last = 0; 118 for (i=1; i<=right; i++) { 119 if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);} 120 } 121 SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp); 122 PetscCall(PetscSortRealWithArrayInt_Private(v,V,last-1)); 123 PetscCall(PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1))); 124 PetscFunctionReturn(0); 125 } 126 /*@ 127 PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order; 128 changes a second PetscInt array to match the sorted first array. 129 130 Not Collective 131 132 Input Parameters: 133 + n - number of values 134 . i - array of integers 135 - I - second array of integers 136 137 Level: intermediate 138 139 .seealso: `PetscSortReal()` 140 @*/ 141 PetscErrorCode PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[]) 142 { 143 PetscInt j,k,itmp; 144 PetscReal rk,rtmp; 145 146 PetscFunctionBegin; 147 PetscValidRealPointer(r,2); 148 PetscValidIntPointer(Ii,3); 149 if (n<8) { 150 for (k=0; k<n; k++) { 151 rk = r[k]; 152 for (j=k+1; j<n; j++) { 153 if (rk > r[j]) { 154 SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp); 155 rk = r[k]; 156 } 157 } 158 } 159 } else { 160 PetscCall(PetscSortRealWithArrayInt_Private(r,Ii,n-1)); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 /*@ 166 PetscFindReal - Finds a PetscReal in a sorted array of PetscReals 167 168 Not Collective 169 170 Input Parameters: 171 + key - the value to locate 172 . n - number of values in the array 173 . ii - array of values 174 - eps - tolerance used to compare 175 176 Output Parameter: 177 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 178 179 Level: intermediate 180 181 .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()` 182 @*/ 183 PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc) 184 { 185 PetscInt lo = 0,hi = n; 186 187 PetscFunctionBegin; 188 PetscValidIntPointer(loc,5); 189 if (!n) {*loc = -1; PetscFunctionReturn(0);} 190 PetscValidRealPointer(t,3); 191 PetscCheckSorted(n,t); 192 while (hi - lo > 1) { 193 PetscInt mid = lo + (hi - lo)/2; 194 if (key < t[mid]) hi = mid; 195 else lo = mid; 196 } 197 *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1); 198 PetscFunctionReturn(0); 199 } 200 201 /*@ 202 PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries 203 204 Not Collective 205 206 Input Parameters: 207 + n - number of values 208 - v - array of doubles 209 210 Output Parameter: 211 . n - number of non-redundant values 212 213 Level: intermediate 214 215 .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()` 216 @*/ 217 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[]) 218 { 219 PetscInt i,s = 0,N = *n, b = 0; 220 221 PetscFunctionBegin; 222 PetscCall(PetscSortReal(N,v)); 223 for (i=0; i<N-1; i++) { 224 if (v[b+s+1] != v[b]) { 225 v[b+1] = v[b+s+1]; b++; 226 } else s++; 227 } 228 *n = N - s; 229 PetscFunctionReturn(0); 230 } 231 232 /*@ 233 PetscSortSplit - Quick-sort split of an array of PetscScalars in place. 234 235 Not Collective 236 237 Input Parameters: 238 + ncut - splitig index 239 - n - number of values to sort 240 241 Input/Output Parameters: 242 + a - array of values, on output the values are permuted 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 - index for array a, on output permuted accordingly 246 247 Level: intermediate 248 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; 256 257 PetscFunctionBegin; 258 first = 0; 259 last = n-1; 260 if (ncut < first || ncut > last) PetscFunctionReturn(0); 261 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 } 277 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; 282 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 PetscFunctionReturn(0); 289 } 290 291 /*@ 292 PetscSortSplitReal - Quick-sort split of an array of PetscReals in place. 293 294 Not Collective 295 296 Input Parameters: 297 + ncut - splitig index 298 - n - number of values to sort 299 300 Input/Output Parameters: 301 + a - array of values, on output the values are permuted such that its elements satisfy: 302 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 303 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 304 - idx - index for array a, on output permuted accordingly 305 306 Level: intermediate 307 308 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()` 309 @*/ 310 PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[]) 311 { 312 PetscInt i,mid,last,itmp,j,first; 313 PetscReal d,tmp; 314 PetscReal abskey; 315 316 PetscFunctionBegin; 317 first = 0; 318 last = n-1; 319 if (ncut < first || ncut > last) PetscFunctionReturn(0); 320 321 while (1) { 322 mid = first; 323 d = a[mid]; 324 abskey = PetscAbsReal(d); 325 i = last; 326 for (j = first + 1; j <= i; ++j) { 327 d = a[j]; 328 if (PetscAbsReal(d) >= abskey) { 329 ++mid; 330 /* interchange */ 331 tmp = a[mid]; itmp = idx[mid]; 332 a[mid] = a[j]; idx[mid] = idx[j]; 333 a[j] = tmp; idx[j] = itmp; 334 } 335 } 336 337 /* interchange */ 338 tmp = a[mid]; itmp = idx[mid]; 339 a[mid] = a[first]; idx[mid] = idx[first]; 340 a[first] = tmp; idx[first] = itmp; 341 342 /* test for while loop */ 343 if (mid == ncut) break; 344 else if (mid > ncut) last = mid - 1; 345 else first = mid + 1; 346 } 347 PetscFunctionReturn(0); 348 } 349