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