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(PetscInt 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, PetscInt right) 42 { 43 PetscInt 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(PetscInt n, PetscReal v[]) 87 { 88 PetscInt j, k; 89 PetscReal tmp, vk; 90 91 PetscFunctionBegin; 92 PetscAssertPointer(v, 2); 93 if (n < 8) { 94 for (k = 0; k < n; k++) { 95 vk = v[k]; 96 for (j = k + 1; j < n; j++) { 97 if (vk > v[j]) { 98 SWAP(v[k], v[j], tmp); 99 vk = v[k]; 100 } 101 } 102 } 103 } else PetscCall(PetscSortReal_Private(v, n - 1)); 104 PetscFunctionReturn(PETSC_SUCCESS); 105 } 106 107 #define SWAP2ri(a, b, c, d, rt, it) \ 108 do { \ 109 rt = a; \ 110 a = b; \ 111 b = rt; \ 112 it = c; \ 113 c = d; \ 114 d = it; \ 115 } while (0) 116 117 /* modified from PetscSortIntWithArray_Private */ 118 static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v, PetscInt *V, PetscInt right) 119 { 120 PetscInt i, last, itmp; 121 PetscReal rvl, rtmp; 122 123 PetscFunctionBegin; 124 if (right <= 1) { 125 if (right == 1) { 126 if (v[0] > v[1]) SWAP2ri(v[0], v[1], V[0], V[1], rtmp, itmp); 127 } 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 SWAP2ri(v[0], v[right / 2], V[0], V[right / 2], rtmp, itmp); 131 rvl = v[0]; 132 last = 0; 133 for (i = 1; i <= right; i++) { 134 if (v[i] < rvl) { 135 last++; 136 SWAP2ri(v[last], v[i], V[last], V[i], rtmp, itmp); 137 } 138 } 139 SWAP2ri(v[0], v[last], V[0], V[last], rtmp, itmp); 140 PetscCall(PetscSortRealWithArrayInt_Private(v, V, last - 1)); 141 PetscCall(PetscSortRealWithArrayInt_Private(v + last + 1, V + last + 1, right - (last + 1))); 142 PetscFunctionReturn(PETSC_SUCCESS); 143 } 144 /*@ 145 PetscSortRealWithArrayInt - Sorts an array of `PetscReal` in place in increasing order; 146 changes a second `PetscInt` array to match the sorted first array. 147 148 Not Collective 149 150 Input Parameters: 151 + n - number of values 152 . Ii - array of integers 153 - r - second array of integers 154 155 Level: intermediate 156 157 .seealso: `PetscSortReal()` 158 @*/ 159 PetscErrorCode PetscSortRealWithArrayInt(PetscInt n, PetscReal r[], PetscInt Ii[]) 160 { 161 PetscInt j, k, itmp; 162 PetscReal rk, rtmp; 163 164 PetscFunctionBegin; 165 PetscAssertPointer(r, 2); 166 PetscAssertPointer(Ii, 3); 167 if (n < 8) { 168 for (k = 0; k < n; k++) { 169 rk = r[k]; 170 for (j = k + 1; j < n; j++) { 171 if (rk > r[j]) { 172 SWAP2ri(r[k], r[j], Ii[k], Ii[j], rtmp, itmp); 173 rk = r[k]; 174 } 175 } 176 } 177 } else { 178 PetscCall(PetscSortRealWithArrayInt_Private(r, Ii, n - 1)); 179 } 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 /*@ 184 PetscFindReal - Finds a `PetscReal` in a sorted array of `PetscReal`s 185 186 Not Collective 187 188 Input Parameters: 189 + key - the value to locate 190 . n - number of values in the array 191 . t - array of values 192 - eps - tolerance used to compare 193 194 Output Parameter: 195 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 196 197 Level: intermediate 198 199 .seealso: `PetscSortReal()`, `PetscSortRealWithArrayInt()` 200 @*/ 201 PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc) 202 { 203 PetscInt lo = 0, hi = n; 204 205 PetscFunctionBegin; 206 PetscAssertPointer(loc, 5); 207 if (!n) { 208 *loc = -1; 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 PetscAssertPointer(t, 3); 212 while (hi - lo > 1) { 213 PetscInt mid = lo + (hi - lo) / 2; 214 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]); 215 if (key < t[mid]) hi = mid; 216 else lo = mid; 217 } 218 *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221 222 /*@ 223 PetscSortRemoveDupsReal - Sorts an array of `PetscReal` in place in increasing order and removes all duplicate entries 224 225 Not Collective 226 227 Input Parameters: 228 + n - number of values 229 - v - array of doubles 230 231 Output Parameter: 232 . n - number of non-redundant values 233 234 Level: intermediate 235 236 .seealso: `PetscSortReal()`, `PetscSortRemoveDupsInt()` 237 @*/ 238 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n, PetscReal v[]) 239 { 240 PetscInt i, s = 0, N = *n, b = 0; 241 242 PetscFunctionBegin; 243 PetscCall(PetscSortReal(N, v)); 244 for (i = 0; i < N - 1; i++) { 245 if (v[b + s + 1] != v[b]) { 246 v[b + 1] = v[b + s + 1]; 247 b++; 248 } else s++; 249 } 250 *n = N - s; 251 PetscFunctionReturn(PETSC_SUCCESS); 252 } 253 254 /*@ 255 PetscSortSplit - Quick-sort split of an array of `PetscScalar`s in place. 256 257 Not Collective 258 259 Input Parameters: 260 + ncut - splitting index 261 - n - number of values to sort 262 263 Input/Output Parameters: 264 + a - array of values, on output the values are permuted such that its elements satisfy: 265 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 266 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 267 - idx - index for array a, on output permuted accordingly 268 269 Level: intermediate 270 271 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()` 272 @*/ 273 PetscErrorCode PetscSortSplit(PetscInt ncut, PetscInt n, PetscScalar a[], PetscInt idx[]) 274 { 275 PetscInt i, mid, last, itmp, j, first; 276 PetscScalar d, tmp; 277 PetscReal abskey; 278 279 PetscFunctionBegin; 280 first = 0; 281 last = n - 1; 282 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS); 283 284 while (1) { 285 mid = first; 286 d = a[mid]; 287 abskey = PetscAbsScalar(d); 288 i = last; 289 for (j = first + 1; j <= i; ++j) { 290 d = a[j]; 291 if (PetscAbsScalar(d) >= abskey) { 292 ++mid; 293 /* interchange */ 294 tmp = a[mid]; 295 itmp = idx[mid]; 296 a[mid] = a[j]; 297 idx[mid] = idx[j]; 298 a[j] = tmp; 299 idx[j] = itmp; 300 } 301 } 302 303 /* interchange */ 304 tmp = a[mid]; 305 itmp = idx[mid]; 306 a[mid] = a[first]; 307 idx[mid] = idx[first]; 308 a[first] = tmp; 309 idx[first] = itmp; 310 311 /* test for while loop */ 312 if (mid == ncut) break; 313 else if (mid > ncut) last = mid - 1; 314 else first = mid + 1; 315 } 316 PetscFunctionReturn(PETSC_SUCCESS); 317 } 318 319 /*@ 320 PetscSortSplitReal - Quick-sort split of an array of `PetscReal`s in place. 321 322 Not Collective 323 324 Input Parameters: 325 + ncut - splitting index 326 - n - number of values to sort 327 328 Input/Output Parameters: 329 + a - array of values, on output the values are permuted such that its elements satisfy: 330 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 331 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 332 - idx - index for array a, on output permuted accordingly 333 334 Level: intermediate 335 336 .seealso: `PetscSortInt()`, `PetscSortRealWithPermutation()` 337 @*/ 338 PetscErrorCode PetscSortSplitReal(PetscInt ncut, PetscInt n, PetscReal a[], PetscInt idx[]) 339 { 340 PetscInt i, mid, last, itmp, j, first; 341 PetscReal d, tmp; 342 PetscReal abskey; 343 344 PetscFunctionBegin; 345 first = 0; 346 last = n - 1; 347 if (ncut < first || ncut > last) PetscFunctionReturn(PETSC_SUCCESS); 348 349 while (1) { 350 mid = first; 351 d = a[mid]; 352 abskey = PetscAbsReal(d); 353 i = last; 354 for (j = first + 1; j <= i; ++j) { 355 d = a[j]; 356 if (PetscAbsReal(d) >= abskey) { 357 ++mid; 358 /* interchange */ 359 tmp = a[mid]; 360 itmp = idx[mid]; 361 a[mid] = a[j]; 362 idx[mid] = idx[j]; 363 a[j] = tmp; 364 idx[j] = itmp; 365 } 366 } 367 368 /* interchange */ 369 tmp = a[mid]; 370 itmp = idx[mid]; 371 a[mid] = a[first]; 372 idx[mid] = idx[first]; 373 a[first] = tmp; 374 idx[first] = itmp; 375 376 /* test for while loop */ 377 if (mid == ncut) break; 378 else if (mid > ncut) last = mid - 1; 379 else first = mid + 1; 380 } 381 PetscFunctionReturn(PETSC_SUCCESS); 382 } 383