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