Actual source code: sortd.c


  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 comparison routines.

  7:  */
  8: #include <petscsys.h>
  9: #include <petsc/private/petscimpl.h>

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

 13: /*@
 14:    PetscSortedReal - Determines whether the array is sorted.

 16:    Not Collective

 18:    Input Parameters:
 19: +  n  - number of values
 20: -  X  - array of integers

 22:    Output Parameters:
 23: .  sorted - flag whether the array is sorted

 25:    Level: intermediate

 27: .seealso: PetscSortReal(), PetscSortedInt(), PetscSortedMPIInt()
 28: @*/
 29: PetscErrorCode  PetscSortedReal(PetscInt n,const PetscReal X[],PetscBool *sorted)
 30: {
 31:   PetscSorted(n,X,*sorted);
 32:   return 0;
 33: }

 35: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
 36: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
 37: {
 38:   PetscInt  i,last;
 39:   PetscReal vl,tmp;

 41:   if (right <= 1) {
 42:     if (right == 1) {
 43:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
 44:     }
 45:     return 0;
 46:   }
 47:   SWAP(v[0],v[right/2],tmp);
 48:   vl   = v[0];
 49:   last = 0;
 50:   for (i=1; i<=right; i++) {
 51:     if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
 52:   }
 53:   SWAP(v[0],v[last],tmp);
 54:   PetscSortReal_Private(v,last-1);
 55:   PetscSortReal_Private(v+last+1,right-(last+1));
 56:   return 0;
 57: }

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

 62:    Not Collective

 64:    Input Parameters:
 65: +  n  - number of values
 66: -  v  - array of doubles

 68:    Notes:
 69:    This function serves as an alternative to PetscRealSortSemiOrdered(), and may perform faster especially if the array
 70:    is completely random. There are exceptions to this and so it is __highly__ recommended that the user benchmark their
 71:    code to see which routine is fastest.

 73:    Level: intermediate

 75: .seealso: PetscRealSortSemiOrdered(), PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
 76: @*/
 77: PetscErrorCode  PetscSortReal(PetscInt n,PetscReal v[])
 78: {
 79:   PetscInt  j,k;
 80:   PetscReal tmp,vk;

 83:   if (n<8) {
 84:     for (k=0; k<n; k++) {
 85:       vk = v[k];
 86:       for (j=k+1; j<n; j++) {
 87:         if (vk > v[j]) {
 88:           SWAP(v[k],v[j],tmp);
 89:           vk = v[k];
 90:         }
 91:       }
 92:     }
 93:   } else PetscSortReal_Private(v,n-1);
 94:   return 0;
 95: }

 97: #define SWAP2ri(a,b,c,d,rt,it) {rt=a;a=b;b=rt;it=c;c=d;d=it;}

 99: /* modified from PetscSortIntWithArray_Private */
100: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
101: {
102:   PetscInt       i,last,itmp;
103:   PetscReal      rvl,rtmp;

105:   if (right <= 1) {
106:     if (right == 1) {
107:       if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
108:     }
109:     return 0;
110:   }
111:   SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
112:   rvl  = v[0];
113:   last = 0;
114:   for (i=1; i<=right; i++) {
115:     if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
116:   }
117:   SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
118:   PetscSortRealWithArrayInt_Private(v,V,last-1);
119:   PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
120:   return 0;
121: }
122: /*@
123:    PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
124:        changes a second PetscInt array to match the sorted first array.

126:    Not Collective

128:    Input Parameters:
129: +  n  - number of values
130: .  i  - array of integers
131: -  I - second array of integers

133:    Level: intermediate

135: .seealso: PetscSortReal()
136: @*/
137: PetscErrorCode  PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
138: {
139:   PetscInt       j,k,itmp;
140:   PetscReal      rk,rtmp;

144:   if (n<8) {
145:     for (k=0; k<n; k++) {
146:       rk = r[k];
147:       for (j=k+1; j<n; j++) {
148:         if (rk > r[j]) {
149:           SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
150:           rk = r[k];
151:         }
152:       }
153:     }
154:   } else {
155:     PetscSortRealWithArrayInt_Private(r,Ii,n-1);
156:   }
157:   return 0;
158: }

160: /*@
161:   PetscFindReal - Finds a PetscReal in a sorted array of PetscReals

163:    Not Collective

165:    Input Parameters:
166: +  key - the value to locate
167: .  n   - number of values in the array
168: .  ii  - array of values
169: -  eps - tolerance used to compare

171:    Output Parameter:
172: .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go

174:    Level: intermediate

176: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
177: @*/
178: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
179: {
180:   PetscInt lo = 0,hi = n;

183:   if (!n) {*loc = -1; return 0;}
186:   while (hi - lo > 1) {
187:     PetscInt mid = lo + (hi - lo)/2;
188:     if (key < t[mid]) hi = mid;
189:     else              lo = mid;
190:   }
191:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
192:   return 0;
193: }

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

198:    Not Collective

200:    Input Parameters:
201: +  n  - number of values
202: -  v  - array of doubles

204:    Output Parameter:
205: .  n - number of non-redundant values

207:    Level: intermediate

209: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
210: @*/
211: PetscErrorCode  PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
212: {
213:   PetscInt       i,s = 0,N = *n, b = 0;

215:   PetscSortReal(N,v);
216:   for (i=0; i<N-1; i++) {
217:     if (v[b+s+1] != v[b]) {
218:       v[b+1] = v[b+s+1]; b++;
219:     } else s++;
220:   }
221:   *n = N - s;
222:   return 0;
223: }

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

228:    Not Collective

230:    Input Parameters:
231: +  ncut  - splitig index
232: -  n     - number of values to sort

234:    Input/Output Parameters:
235: +  a     - array of values, on output the values are permuted such that its elements satisfy:
236:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
237:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
238: -  idx   - index for array a, on output permuted accordingly

240:    Level: intermediate

242: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
243: @*/
244: PetscErrorCode  PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
245: {
246:   PetscInt    i,mid,last,itmp,j,first;
247:   PetscScalar d,tmp;
248:   PetscReal   abskey;

250:   first = 0;
251:   last  = n-1;
252:   if (ncut < first || ncut > last) return 0;

254:   while (1) {
255:     mid    = first;
256:     d      = a[mid];
257:     abskey = PetscAbsScalar(d);
258:     i      = last;
259:     for (j = first + 1; j <= i; ++j) {
260:       d = a[j];
261:       if (PetscAbsScalar(d) >= abskey) {
262:         ++mid;
263:         /* interchange */
264:         tmp = a[mid];  itmp = idx[mid];
265:         a[mid] = a[j]; idx[mid] = idx[j];
266:         a[j] = tmp;    idx[j] = itmp;
267:       }
268:     }

270:     /* interchange */
271:     tmp = a[mid];      itmp = idx[mid];
272:     a[mid] = a[first]; idx[mid] = idx[first];
273:     a[first] = tmp;    idx[first] = itmp;

275:     /* test for while loop */
276:     if (mid == ncut) break;
277:     else if (mid > ncut) last = mid - 1;
278:     else first = mid + 1;
279:   }
280:   return 0;
281: }

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

286:    Not Collective

288:    Input Parameters:
289: +  ncut  - splitig index
290: -  n     - number of values to sort

292:    Input/Output Parameters:
293: +  a     - array of values, on output the values are permuted such that its elements satisfy:
294:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
295:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
296: -  idx   - index for array a, on output permuted accordingly

298:    Level: intermediate

300: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
301: @*/
302: PetscErrorCode  PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
303: {
304:   PetscInt  i,mid,last,itmp,j,first;
305:   PetscReal d,tmp;
306:   PetscReal abskey;

308:   first = 0;
309:   last  = n-1;
310:   if (ncut < first || ncut > last) return 0;

312:   while (1) {
313:     mid    = first;
314:     d      = a[mid];
315:     abskey = PetscAbsReal(d);
316:     i      = last;
317:     for (j = first + 1; j <= i; ++j) {
318:       d = a[j];
319:       if (PetscAbsReal(d) >= abskey) {
320:         ++mid;
321:         /* interchange */
322:         tmp = a[mid];  itmp = idx[mid];
323:         a[mid] = a[j]; idx[mid] = idx[j];
324:         a[j] = tmp;    idx[j] = itmp;
325:       }
326:     }

328:     /* interchange */
329:     tmp = a[mid];      itmp = idx[mid];
330:     a[mid] = a[first]; idx[mid] = idx[first];
331:     a[first] = tmp;    idx[first] = itmp;

333:     /* test for while loop */
334:     if (mid == ncut) break;
335:     else if (mid > ncut) last = mid - 1;
336:     else first = mid + 1;
337:   }
338:   return 0;
339: }