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