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 Output Parameters: 238 + n - number of non-redundant values 239 - v - array of non-redundant values 240 241 Level: intermediate 242 243 .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()` 244 @*/ 245 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[]) 246 { 247 PetscInt i, s = 0, N = *n, b = 0; 248 249 PetscFunctionBegin; 250 PetscCall(PetscSortReal(N, v)); 251 for (i = 0; i < N - 1; i++) { 252 if (v[b + s + 1] != v[b]) { 253 v[b + 1] = v[b + s + 1]; 254 b++; 255 } else s++; 256 } 257 *n = N - s; 258 PetscFunctionReturn(PETSC_SUCCESS); 259 } 260 261 /*@ 262 PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place. 263 264 Not Collective 265 266 Input Parameters: 267 + ncut - splitting index 268 - n - number of values to sort 269 270 Input/Output Parameters: 271 + a - array of values, on output the values are permuted such that its elements satisfy: 272 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 273 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 274 - idx - index for array a, on output permuted accordingly 275 276 Level: intermediate 277 278 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()` 279 @*/ 280 PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[]) 281 { 282 PetscInt i, mid, last, itmp, j, first; 283 PetscScalar d, tmp; 284 PetscReal abskey; 285 286 PetscFunctionBegin; 287 first = 0; 288 last = n - 1; 289 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS); 290 291 while (1) { 292 mid = first; 293 d = a[mid]; 294 abskey = PetscAbsScalar(d); 295 i = last; 296 for (j = first + 1; j <= i; ++j) { 297 d = a[j]; 298 if (PetscAbsScalar(d) >= abskey) { 299 ++mid; 300 /* interchange */ 301 tmp = a[mid]; 302 itmp = idx[mid]; 303 a[mid] = a[j]; 304 idx[mid] = idx[j]; 305 a[j] = tmp; 306 idx[j] = itmp; 307 } 308 } 309 310 /* interchange */ 311 tmp = a[mid]; 312 itmp = idx[mid]; 313 a[mid] = a[first]; 314 idx[mid] = idx[first]; 315 a[first] = tmp; 316 idx[first] = itmp; 317 318 /* test for while loop */ 319 if (mid == ncut) break; 320 else if (mid > ncut) last = mid - 1; 321 else first = mid + 1; 322 } 323 PetscFunctionReturn(PETSC_SUCCESS); 324 } 325 326 /*@ 327 PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place. 328 329 Not Collective 330 331 Input Parameters: 332 + ncut - splitting index 333 - n - number of values to sort 334 335 Input/Output Parameters: 336 + a - array of values, on output the values are permuted such that its elements satisfy: 337 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 338 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 339 - idx - index for array a, on output permuted accordingly 340 341 Level: intermediate 342 343 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()` 344 @*/ 345 PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[]) 346 { 347 PetscInt i, mid, last, itmp, j, first; 348 PetscReal d, tmp; 349 PetscReal abskey; 350 351 PetscFunctionBegin; 352 first = 0; 353 last = n - 1; 354 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS); 355 356 while (1) { 357 mid = first; 358 d = a[mid]; 359 abskey = PetscAbsReal(d); 360 i = last; 361 for (j = first + 1; j <= i; ++j) { 362 d = a[j]; 363 if (PetscAbsReal(d) >= abskey) { 364 ++mid; 365 /* interchange */ 366 tmp = a[mid]; 367 itmp = idx[mid]; 368 a[mid] = a[j]; 369 idx[mid] = idx[j]; 370 a[j] = tmp; 371 idx[j] = itmp; 372 } 373 } 374 375 /* interchange */ 376 tmp = a[mid]; 377 itmp = idx[mid]; 378 a[mid] = a[first]; 379 idx[mid] = idx[first]; 380 a[first] = tmp; 381 idx[first] = itmp; 382 383 /* test for while loop */ 384 if (mid == ncut) break; 385 else if (mid > ncut) last = mid - 1; 386 else first = mid + 1; 387 } 388 PetscFunctionReturn(PETSC_SUCCESS); 389 } 390