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