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