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: {
 32:   PetscSorted(n,X,*sorted);
 33:   return(0);
 34: }

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

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

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

 64:    Not Collective

 66:    Input Parameters:
 67: +  n  - number of values
 68: -  v  - array of doubles

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

 75:    Level: intermediate

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

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

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

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

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

131:    Not Collective

133:    Input Parameters:
134: +  n  - number of values
135: .  i  - array of integers
136: -  I - second array of integers

138:    Level: intermediate

140: .seealso: PetscSortReal()
141: @*/
142: PetscErrorCode  PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
143: {
145:   PetscInt       j,k,itmp;
146:   PetscReal      rk,rtmp;

151:   if (n<8) {
152:     for (k=0; k<n; k++) {
153:       rk = r[k];
154:       for (j=k+1; j<n; j++) {
155:         if (rk > r[j]) {
156:           SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
157:           rk = r[k];
158:         }
159:       }
160:     }
161:   } else {
162:     PetscSortRealWithArrayInt_Private(r,Ii,n-1);
163:   }
164:   return(0);
165: }

167: /*@
168:   PetscFindReal - Finds a PetscReal in a sorted array of PetscReals

170:    Not Collective

172:    Input Parameters:
173: +  key - the value to locate
174: .  n   - number of values in the array
175: .  ii  - array of values
176: -  eps - tolerance used to compare

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

181:    Level: intermediate

183: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
184: @*/
185: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
186: {
187:   PetscInt lo = 0,hi = n;

191:   if (!n) {*loc = -1; return(0);}
194:   while (hi - lo > 1) {
195:     PetscInt mid = lo + (hi - lo)/2;
196:     if (key < t[mid]) hi = mid;
197:     else              lo = mid;
198:   }
199:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
200:   return(0);
201: }

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

206:    Not Collective

208:    Input Parameters:
209: +  n  - number of values
210: -  v  - array of doubles

212:    Output Parameter:
213: .  n - number of non-redundant values

215:    Level: intermediate

217: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
218: @*/
219: PetscErrorCode  PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
220: {
222:   PetscInt       i,s = 0,N = *n, b = 0;

225:   PetscSortReal(N,v);
226:   for (i=0; i<N-1; i++) {
227:     if (v[b+s+1] != v[b]) {
228:       v[b+1] = v[b+s+1]; b++;
229:     } else s++;
230:   }
231:   *n = N - s;
232:   return(0);
233: }

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

238:    Not Collective

240:    Input Parameters:
241: +  ncut  - splitig index
242: -  n     - number of values to sort

244:    Input/Output Parameters:
245: +  a     - array of values, on output the values are permuted such that its elements satisfy:
246:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
247:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
248: -  idx   - index for array a, on output permuted accordingly

250:    Level: intermediate

252: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
253: @*/
254: PetscErrorCode  PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
255: {
256:   PetscInt    i,mid,last,itmp,j,first;
257:   PetscScalar d,tmp;
258:   PetscReal   abskey;

261:   first = 0;
262:   last  = n-1;
263:   if (ncut < first || ncut > last) return(0);

265:   while (1) {
266:     mid    = first;
267:     d      = a[mid];
268:     abskey = PetscAbsScalar(d);
269:     i      = last;
270:     for (j = first + 1; j <= i; ++j) {
271:       d = a[j];
272:       if (PetscAbsScalar(d) >= abskey) {
273:         ++mid;
274:         /* interchange */
275:         tmp = a[mid];  itmp = idx[mid];
276:         a[mid] = a[j]; idx[mid] = idx[j];
277:         a[j] = tmp;    idx[j] = itmp;
278:       }
279:     }

281:     /* interchange */
282:     tmp = a[mid];      itmp = idx[mid];
283:     a[mid] = a[first]; idx[mid] = idx[first];
284:     a[first] = tmp;    idx[first] = itmp;

286:     /* test for while loop */
287:     if (mid == ncut) break;
288:     else if (mid > ncut) last = mid - 1;
289:     else first = mid + 1;
290:   }
291:   return(0);
292: }

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

297:    Not Collective

299:    Input Parameters:
300: +  ncut  - splitig index
301: -  n     - number of values to sort

303:    Input/Output Parameters:
304: +  a     - array of values, on output the values are permuted such that its elements satisfy:
305:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
306:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
307: -  idx   - index for array a, on output permuted accordingly

309:    Level: intermediate

311: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
312: @*/
313: PetscErrorCode  PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
314: {
315:   PetscInt  i,mid,last,itmp,j,first;
316:   PetscReal d,tmp;
317:   PetscReal abskey;

320:   first = 0;
321:   last  = n-1;
322:   if (ncut < first || ncut > last) return(0);

324:   while (1) {
325:     mid    = first;
326:     d      = a[mid];
327:     abskey = PetscAbsReal(d);
328:     i      = last;
329:     for (j = first + 1; j <= i; ++j) {
330:       d = a[j];
331:       if (PetscAbsReal(d) >= abskey) {
332:         ++mid;
333:         /* interchange */
334:         tmp = a[mid];  itmp = idx[mid];
335:         a[mid] = a[j]; idx[mid] = idx[j];
336:         a[j] = tmp;    idx[j] = itmp;
337:       }
338:     }

340:     /* interchange */
341:     tmp = a[mid];      itmp = idx[mid];
342:     a[mid] = a[first]; idx[mid] = idx[first];
343:     a[first] = tmp;    idx[first] = itmp;

345:     /* test for while loop */
346:     if (mid == ncut) break;
347:     else if (mid > ncut) last = mid - 1;
348:     else first = mid + 1;
349:   }
350:   return(0);
351: }