Actual source code: sortd.c

petsc-master 2020-07-05
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:    Level: intermediate

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

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

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

 97: /* modified from PetscSortIntWithArray_Private */
 98: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
 99: {
101:   PetscInt       i,last,itmp;
102:   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: {
140:   PetscInt       j,k,itmp;
141:   PetscReal      rk,rtmp;

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

162: /*@
163:   PetscFindReal - Finds a PetscReal in a sorted array of PetscReals

165:    Not Collective

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

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

176:    Level: intermediate

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

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

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

201:    Not Collective

203:    Input Parameters:
204: +  n  - number of values
205: -  v  - array of doubles

207:    Output Parameter:
208: .  n - number of non-redundant values

210:    Level: intermediate

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

220:   PetscSortReal(N,v);
221:   for (i=0; i<N-1; i++) {
222:     if (v[b+s+1] != v[b]) {
223:       v[b+1] = v[b+s+1]; b++;
224:     } else s++;
225:   }
226:   *n = N - s;
227:   return(0);
228: }

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

233:    Not Collective

235:    Input Parameters:
236: +  ncut  - splitig index
237: .  n     - number of values to sort
238: .  a     - array of values
239: -  idx   - index for array a

241:    Output Parameters:
242: +  a     - permuted array of values such that its elements satisfy:
243:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
244:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
245: -  idx   - permuted index of array a

247:    Level: intermediate

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

258:   first = 0;
259:   last  = n-1;
260:   if (ncut < first || ncut > last) return(0);

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

278:     /* interchange */
279:     tmp = a[mid];      itmp = idx[mid];
280:     a[mid] = a[first]; idx[mid] = idx[first];
281:     a[first] = tmp;    idx[first] = itmp;

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

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

294:    Not Collective

296:    Input Parameters:
297: +  ncut  - splitig index
298: .  n     - number of values to sort
299: .  a     - array of values in PetscReal
300: -  idx   - index for array a

302:    Output Parameters:
303: +  a     - permuted array of real values such that its elements satisfy:
304:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
305:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
306: -  idx   - permuted index of array a

308:    Level: intermediate

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

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

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

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

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