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 { 73 PetscSortReal_Private(v,n-1); 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "PetscSortSplit" 80 /*@ 81 PetscSortSplit - Quick-sort split of an array of PetscScalars in place. 82 83 Not Collective 84 85 Input Parameters: 86 + ncut - splitig index 87 . n - number of values to sort 88 . a - array of values 89 - idx - index for array a 90 91 Output Parameters: 92 + a - permuted array of values such that its elements satisfy: 93 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 94 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 95 - idx - permuted index of array a 96 97 Level: intermediate 98 99 Concepts: sorting^doubles 100 101 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 102 @*/ 103 PetscErrorCode PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[]) 104 { 105 PetscInt i,mid,last,itmp,j,first; 106 PetscScalar d,tmp; 107 PetscReal abskey; 108 109 PetscFunctionBegin; 110 first = 0; 111 last = n-1; 112 if (ncut < first || ncut > last) PetscFunctionReturn(0); 113 114 while (1){ 115 mid = first; 116 abskey = (d = a[mid],PetscAbsScalar(d)); 117 i = last; 118 for (j = first + 1; j <= i; ++j) { 119 if ((d = a[j],PetscAbsScalar(d)) >= abskey) { 120 ++mid; 121 /* interchange */ 122 tmp = a[mid]; itmp = idx[mid]; 123 a[mid] = a[j]; idx[mid] = idx[j]; 124 a[j] = tmp; idx[j] = itmp; 125 } 126 } 127 128 /* interchange */ 129 tmp = a[mid]; itmp = idx[mid]; 130 a[mid] = a[first]; idx[mid] = idx[first]; 131 a[first] = tmp; idx[first] = itmp; 132 133 /* test for while loop */ 134 if (mid == ncut) { 135 break; 136 } else if (mid > ncut){ 137 last = mid - 1; 138 } else { 139 first = mid + 1; 140 } 141 } 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "PetscSortSplitReal" 147 /*@ 148 PetscSortSplitReal - Quick-sort split of an array of PetscReals in place. 149 150 Not Collective 151 152 Input Parameters: 153 + ncut - splitig index 154 . n - number of values to sort 155 . a - array of values in PetscReal 156 - idx - index for array a 157 158 Output Parameters: 159 + a - permuted array of real values such that its elements satisfy: 160 abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 161 abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 162 - idx - permuted index of array a 163 164 Level: intermediate 165 166 Concepts: sorting^doubles 167 168 .seealso: PetscSortInt(), PetscSortRealWithPermutation() 169 @*/ 170 PetscErrorCode PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[]) 171 { 172 PetscInt i,mid,last,itmp,j,first; 173 PetscReal d,tmp; 174 PetscReal abskey; 175 176 PetscFunctionBegin; 177 first = 0; 178 last = n-1; 179 if (ncut < first || ncut > last) PetscFunctionReturn(0); 180 181 while (1){ 182 mid = first; 183 abskey = (d = a[mid],PetscAbsReal(d)); 184 i = last; 185 for (j = first + 1; j <= i; ++j) { 186 if ((d = a[j],PetscAbsReal(d)) >= abskey) { 187 ++mid; 188 /* interchange */ 189 tmp = a[mid]; itmp = idx[mid]; 190 a[mid] = a[j]; idx[mid] = idx[j]; 191 a[j] = tmp; idx[j] = itmp; 192 } 193 } 194 195 /* interchange */ 196 tmp = a[mid]; itmp = idx[mid]; 197 a[mid] = a[first]; idx[mid] = idx[first]; 198 a[first] = tmp; idx[first] = itmp; 199 200 /* test for while loop */ 201 if (mid == ncut) { 202 break; 203 } else if (mid > ncut){ 204 last = mid - 1; 205 } else { 206 first = mid + 1; 207 } 208 } 209 PetscFunctionReturn(0); 210 } 211 212