Actual source code: sortd.c

petsc-master 2020-11-29
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>
  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__ recomended 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
243: .  a     - array of values
244: -  idx   - index for array a

246:    Output Parameters:
247: +  a     - permuted array of values such that its elements satisfy:
248:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
249:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
250: -  idx   - permuted index of array a

252:    Level: intermediate

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

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

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

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

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

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

299:    Not Collective

301:    Input Parameters:
302: +  ncut  - splitig index
303: .  n     - number of values to sort
304: .  a     - array of values in PetscReal
305: -  idx   - index for array a

307:    Output Parameters:
308: +  a     - permuted array of real values such that its elements satisfy:
309:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
310:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
311: -  idx   - permuted index of array a

313:    Level: intermediate

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

324:   first = 0;
325:   last  = n-1;
326:   if (ncut < first || ncut > last) return(0);

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

344:     /* interchange */
345:     tmp = a[mid];      itmp = idx[mid];
346:     a[mid] = a[first]; idx[mid] = idx[first];
347:     a[first] = tmp;    idx[first] = itmp;

349:     /* test for while loop */
350:     if (mid == ncut) break;
351:     else if (mid > ncut) last = mid - 1;
352:     else first = mid + 1;
353:   }
354:   return(0);
355: }