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