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 abskey = (d = a[mid],PetscAbsScalar(d)); 152 i = last; 153 for (j = first + 1; j <= i; ++j) { 154 if ((d = a[j],PetscAbsScalar(d)) >= abskey) { 155 ++mid; 156 /* interchange */ 157 tmp = a[mid]; itmp = idx[mid]; 158 a[mid] = a[j]; idx[mid] = idx[j]; 159 a[j] = tmp; idx[j] = itmp; 160 } 161 } 162 163 /* interchange */ 164 tmp = a[mid]; itmp = idx[mid]; 165 a[mid] = a[first]; idx[mid] = idx[first]; 166 a[first] = tmp; idx[first] = itmp; 167 168 /* test for while loop */ 169 if (mid == ncut) break; 170 else if (mid > ncut) last = mid - 1; 171 else first = mid + 1; 172 } 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "PetscSortSplitReal" 178 /*@ 179 PetscSortSplitReal - Quick-sort split of an array of PetscReals in place. 180 181 Not Collective 182 183 Input Parameters: 184 + ncut - splitig index 185 . n - number of values to sort 186 . a - array of values in PetscReal 187 - idx - index for array a 188 189 Output Parameters: 190 + a - permuted array of real values such that its elements satisfy: 191 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 192 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 193 - idx - permuted index of array a 194 195 Level: intermediate 196 197 Concepts: sorting^doubles 198 199 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 200 @*/ 201 PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[]) 202 { 203 PetscInt i,mid,last,itmp,j,first; 204 PetscReal d,tmp; 205 PetscReal abskey; 206 207 PetscFunctionBegin; 208 first = 0; 209 last = n-1; 210 if (ncut < first || ncut > last) PetscFunctionReturn(0); 211 212 while (1) { 213 mid = first; 214 abskey = (d = a[mid],PetscAbsReal(d)); 215 i = last; 216 for (j = first + 1; j <= i; ++j) { 217 if ((d = a[j],PetscAbsReal(d)) >= abskey) { 218 ++mid; 219 /* interchange */ 220 tmp = a[mid]; itmp = idx[mid]; 221 a[mid] = a[j]; idx[mid] = idx[j]; 222 a[j] = tmp; idx[j] = itmp; 223 } 224 } 225 226 /* interchange */ 227 tmp = a[mid]; itmp = idx[mid]; 228 a[mid] = a[first]; idx[mid] = idx[first]; 229 a[first] = tmp; idx[first] = itmp; 230 231 /* test for while loop */ 232 if (mid == ncut) break; 233 else if (mid > ncut) last = mid - 1; 234 else first = mid + 1; 235 } 236 PetscFunctionReturn(0); 237 } 238 239