Actual source code: sortd.c

petsc-master 2017-04-24
Report Typos and Errors

  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.

  7:  */
  8:  #include <petscsys.h>

 10: #define SWAP(a,b,t) {t=a;a=b;b=t;}

 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;

 19:   if (right <= 1) {
 20:     if (right == 1) {
 21:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
 22:     }
 23:     return(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:   return(0);
 35: }

 37: /*@
 38:    PetscSortReal - Sorts an array of doubles in place in increasing order.

 40:    Not Collective

 42:    Input Parameters:
 43: +  n  - number of values
 44: -  v  - array of doubles

 46:    Level: intermediate

 48:    Concepts: sorting^doubles

 50: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
 51: @*/
 52: PetscErrorCode  PetscSortReal(PetscInt n,PetscReal v[])
 53: {
 54:   PetscInt  j,k;
 55:   PetscReal tmp,vk;

 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:   return(0);
 70: }


 73: /*@
 74:    PetscSortRemoveDupsReal - Sorts an array of doubles in place in increasing order removes all duplicate entries

 76:    Not Collective

 78:    Input Parameters:
 79: +  n  - number of values
 80: -  v  - array of doubles

 82:    Output Parameter:
 83: .  n - number of non-redundant values

 85:    Level: intermediate

 87:    Concepts: sorting^doubles

 89: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
 90: @*/
 91: PetscErrorCode  PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
 92: {
 94:   PetscInt       i,s = 0,N = *n, b = 0;

 97:   PetscSortReal(N,v);
 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:   return(0);
105: }

107: /*@
108:    PetscSortSplit - Quick-sort split of an array of PetscScalars in place.

110:    Not Collective

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

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

124:    Level: intermediate

126:    Concepts: sorting^doubles

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;

137:   first = 0;
138:   last  = n-1;
139:   if (ncut < first || ncut > last) return(0);

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:     }

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;

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:   return(0);
168: }

170: /*@
171:    PetscSortSplitReal - Quick-sort split of an array of PetscReals in place.

173:    Not Collective

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

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

187:    Level: intermediate

189:    Concepts: sorting^doubles

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;

200:   first = 0;
201:   last  = n-1;
202:   if (ncut < first || ncut > last) return(0);

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:     }

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;

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:   return(0);
231: }