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