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 10 #define SWAP(a,b,t) {t=a;a=b;b=t;} 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "PetscSortReal_Private" 14 /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */ 15 static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right) 16 { 17 PetscInt i,last; 18 PetscReal vl,tmp; 19 20 PetscFunctionBegin; 21 if (right <= 1) { 22 if (right == 1) { 23 if (v[0] > v[1]) SWAP(v[0],v[1],tmp); 24 } 25 PetscFunctionReturn(0); 26 } 27 SWAP(v[0],v[right/2],tmp); 28 vl = v[0]; 29 last = 0; 30 for (i=1; i<=right; i++) { 31 if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);} 32 } 33 SWAP(v[0],v[last],tmp); 34 PetscSortReal_Private(v,last-1); 35 PetscSortReal_Private(v+last+1,right-(last+1)); 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "PetscSortReal" 41 /*@ 42 PetscSortReal - Sorts an array of doubles in place in increasing order. 43 44 Not Collective 45 46 Input Parameters: 47 + n - number of values 48 - v - array of doubles 49 50 Level: intermediate 51 52 Concepts: sorting^doubles 53 54 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 55 @*/ 56 PetscErrorCode PetscSortReal(PetscInt n,PetscReal v[]) 57 { 58 PetscInt j,k; 59 PetscReal tmp,vk; 60 61 PetscFunctionBegin; 62 if (n<8) { 63 for (k=0; k<n; k++) { 64 vk = v[k]; 65 for (j=k+1; j<n; j++) { 66 if (vk > v[j]) { 67 SWAP(v[k],v[j],tmp); 68 vk = v[k]; 69 } 70 } 71 } 72 } else PetscSortReal_Private(v,n-1); 73 PetscFunctionReturn(0); 74 } 75 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "PetscSortRemoveDupsReal" 79 /*@ 80 PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries 81 82 Not Collective 83 84 Input Parameters: 85 + n - number of values 86 - v - array of doubles 87 88 Output Parameter: 89 . n - number of non-redundant values 90 91 Level: intermediate 92 93 Concepts: sorting^doubles 94 95 .seealso: PetscSortReal(), PetscSortRemoveDupsInt() 96 @*/ 97 PetscErrorCode PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[]) 98 { 99 PetscErrorCode ierr; 100 PetscInt i,s = 0,N = *n, b = 0; 101 102 PetscFunctionBegin; 103 ierr = PetscSortReal(N,v);CHKERRQ(ierr); 104 for (i=0; i<N-1; i++) { 105 if (v[b+s+1] != v[b]) { 106 v[b+1] = v[b+s+1]; b++; 107 } else s++; 108 } 109 *n = N - s; 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "PetscSortSplit" 115 /*@ 116 PetscSortSplit - Quick-sort split of an array of PetscScalars in place. 117 118 Not Collective 119 120 Input Parameters: 121 + ncut - splitig index 122 . n - number of values to sort 123 . a - array of values 124 - idx - index for array a 125 126 Output Parameters: 127 + a - permuted array of values such that its elements satisfy: 128 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 129 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 130 - idx - permuted index of array a 131 132 Level: intermediate 133 134 Concepts: sorting^doubles 135 136 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 137 @*/ 138 PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[]) 139 { 140 PetscInt i,mid,last,itmp,j,first; 141 PetscScalar d,tmp; 142 PetscReal abskey; 143 144 PetscFunctionBegin; 145 first = 0; 146 last = n-1; 147 if (ncut < first || ncut > last) PetscFunctionReturn(0); 148 149 while (1) { 150 mid = first; 151 d = a[mid]; 152 abskey = PetscAbsScalar(d); 153 i = last; 154 for (j = first + 1; j <= i; ++j) { 155 d = a[j]; 156 if (PetscAbsScalar(d) >= abskey) { 157 ++mid; 158 /* interchange */ 159 tmp = a[mid]; itmp = idx[mid]; 160 a[mid] = a[j]; idx[mid] = idx[j]; 161 a[j] = tmp; idx[j] = itmp; 162 } 163 } 164 165 /* interchange */ 166 tmp = a[mid]; itmp = idx[mid]; 167 a[mid] = a[first]; idx[mid] = idx[first]; 168 a[first] = tmp; idx[first] = itmp; 169 170 /* test for while loop */ 171 if (mid == ncut) break; 172 else if (mid > ncut) last = mid - 1; 173 else first = mid + 1; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "PetscSortSplitReal" 180 /*@ 181 PetscSortSplitReal - Quick-sort split of an array of PetscReals in place. 182 183 Not Collective 184 185 Input Parameters: 186 + ncut - splitig index 187 . n - number of values to sort 188 . a - array of values in PetscReal 189 - idx - index for array a 190 191 Output Parameters: 192 + a - permuted array of real values such that its elements satisfy: 193 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 194 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 195 - idx - permuted index of array a 196 197 Level: intermediate 198 199 Concepts: sorting^doubles 200 201 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 202 @*/ 203 PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[]) 204 { 205 PetscInt i,mid,last,itmp,j,first; 206 PetscReal d,tmp; 207 PetscReal abskey; 208 209 PetscFunctionBegin; 210 first = 0; 211 last = n-1; 212 if (ncut < first || ncut > last) PetscFunctionReturn(0); 213 214 while (1) { 215 mid = first; 216 d = a[mid]; 217 abskey = PetscAbsReal(d); 218 i = last; 219 for (j = first + 1; j <= i; ++j) { 220 d = a[j]; 221 if (PetscAbsReal(d) >= abskey) { 222 ++mid; 223 /* interchange */ 224 tmp = a[mid]; itmp = idx[mid]; 225 a[mid] = a[j]; idx[mid] = idx[j]; 226 a[j] = tmp; idx[j] = itmp; 227 } 228 } 229 230 /* interchange */ 231 tmp = a[mid]; itmp = idx[mid]; 232 a[mid] = a[first]; idx[mid] = idx[first]; 233 a[first] = tmp; idx[first] = itmp; 234 235 /* test for while loop */ 236 if (mid == ncut) break; 237 else if (mid > ncut) last = mid - 1; 238 else first = mid + 1; 239 } 240 PetscFunctionReturn(0); 241 } 242 243