Actual source code: sortd.c

petsc-3.12.0 2019-09-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: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
 14: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
 15: {
 16:   PetscInt  i,last;
 17:   PetscReal vl,tmp;

 20:   if (right <= 1) {
 21:     if (right == 1) {
 22:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
 23:     }
 24:     return(0);
 25:   }
 26:   SWAP(v[0],v[right/2],tmp);
 27:   vl   = v[0];
 28:   last = 0;
 29:   for (i=1; i<=right; i++) {
 30:     if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
 31:   }
 32:   SWAP(v[0],v[last],tmp);
 33:   PetscSortReal_Private(v,last-1);
 34:   PetscSortReal_Private(v+last+1,right-(last+1));
 35:   return(0);
 36: }

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

 41:    Not Collective

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

 47:    Level: intermediate

 49: .seealso: PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
 50: @*/
 51: PetscErrorCode  PetscSortReal(PetscInt n,PetscReal v[])
 52: {
 53:   PetscInt  j,k;
 54:   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: }

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

 74: /* modified from PetscSortIntWithArray_Private */
 75: static PetscErrorCode PetscSortRealWithArrayInt_Private(PetscReal *v,PetscInt *V,PetscInt right)
 76: {
 78:   PetscInt       i,last,itmp;
 79:   PetscReal      rvl,rtmp;

 82:   if (right <= 1) {
 83:     if (right == 1) {
 84:       if (v[0] > v[1]) SWAP2ri(v[0],v[1],V[0],V[1],rtmp,itmp);
 85:     }
 86:     return(0);
 87:   }
 88:   SWAP2ri(v[0],v[right/2],V[0],V[right/2],rtmp,itmp);
 89:   rvl  = v[0];
 90:   last = 0;
 91:   for (i=1; i<=right; i++) {
 92:     if (v[i] < rvl) {last++; SWAP2ri(v[last],v[i],V[last],V[i],rtmp,itmp);}
 93:   }
 94:   SWAP2ri(v[0],v[last],V[0],V[last],rtmp,itmp);
 95:   PetscSortRealWithArrayInt_Private(v,V,last-1);
 96:   PetscSortRealWithArrayInt_Private(v+last+1,V+last+1,right-(last+1));
 97:   return(0);
 98: }
 99: /*@
100:    PetscSortRealWithArrayInt - Sorts an array of PetscReal in place in increasing order;
101:        changes a second PetscInt array to match the sorted first array.

103:    Not Collective

105:    Input Parameters:
106: +  n  - number of values
107: .  i  - array of integers
108: -  I - second array of integers

110:    Level: intermediate

112: .seealso: PetscSortReal()
113: @*/
114: PetscErrorCode  PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
115: {
117:   PetscInt       j,k,itmp;
118:   PetscReal      rk,rtmp;

123:   if (n<8) {
124:     for (k=0; k<n; k++) {
125:       rk = r[k];
126:       for (j=k+1; j<n; j++) {
127:         if (rk > r[j]) {
128:           SWAP2ri(r[k],r[j],Ii[k],Ii[j],rtmp,itmp);
129:           rk = r[k];
130:         }
131:       }
132:     }
133:   } else {
134:     PetscSortRealWithArrayInt_Private(r,Ii,n-1);
135:   }
136:   return(0);
137: }

139: /*@
140:   PetscFindReal - Finds a PetscReal in a sorted array of PetscReals

142:    Not Collective

144:    Input Parameters:
145: +  key - the value to locate
146: .  n   - number of values in the array
147: .  ii  - array of values
148: -  eps - tolerance used to compare

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

153:    Level: intermediate

155: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
156: @*/
157: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
158: {
159:   PetscInt lo = 0,hi = n;

163:   if (!n) {*loc = -1; return(0);}
165:   while (hi - lo > 1) {
166:     PetscInt mid = lo + (hi - lo)/2;
167:     if (key < t[mid]) hi = mid;
168:     else              lo = mid;
169:   }
170:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
171:   return(0);
172: }

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

177:    Not Collective

179:    Input Parameters:
180: +  n  - number of values
181: -  v  - array of doubles

183:    Output Parameter:
184: .  n - number of non-redundant values

186:    Level: intermediate

188: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
189: @*/
190: PetscErrorCode  PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
191: {
193:   PetscInt       i,s = 0,N = *n, b = 0;

196:   PetscSortReal(N,v);
197:   for (i=0; i<N-1; i++) {
198:     if (v[b+s+1] != v[b]) {
199:       v[b+1] = v[b+s+1]; b++;
200:     } else s++;
201:   }
202:   *n = N - s;
203:   return(0);
204: }

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

209:    Not Collective

211:    Input Parameters:
212: +  ncut  - splitig index
213: .  n     - number of values to sort
214: .  a     - array of values
215: -  idx   - index for array a

217:    Output Parameters:
218: +  a     - permuted array of values such that its elements satisfy:
219:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
220:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
221: -  idx   - permuted index of array a

223:    Level: intermediate

225: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
226: @*/
227: PetscErrorCode  PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
228: {
229:   PetscInt    i,mid,last,itmp,j,first;
230:   PetscScalar d,tmp;
231:   PetscReal   abskey;

234:   first = 0;
235:   last  = n-1;
236:   if (ncut < first || ncut > last) return(0);

238:   while (1) {
239:     mid    = first;
240:     d      = a[mid];
241:     abskey = PetscAbsScalar(d);
242:     i      = last;
243:     for (j = first + 1; j <= i; ++j) {
244:       d = a[j];
245:       if (PetscAbsScalar(d) >= abskey) {
246:         ++mid;
247:         /* interchange */
248:         tmp = a[mid];  itmp = idx[mid];
249:         a[mid] = a[j]; idx[mid] = idx[j];
250:         a[j] = tmp;    idx[j] = itmp;
251:       }
252:     }

254:     /* interchange */
255:     tmp = a[mid];      itmp = idx[mid];
256:     a[mid] = a[first]; idx[mid] = idx[first];
257:     a[first] = tmp;    idx[first] = itmp;

259:     /* test for while loop */
260:     if (mid == ncut) break;
261:     else if (mid > ncut) last = mid - 1;
262:     else first = mid + 1;
263:   }
264:   return(0);
265: }

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

270:    Not Collective

272:    Input Parameters:
273: +  ncut  - splitig index
274: .  n     - number of values to sort
275: .  a     - array of values in PetscReal
276: -  idx   - index for array a

278:    Output Parameters:
279: +  a     - permuted array of real values such that its elements satisfy:
280:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
281:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
282: -  idx   - permuted index of array a

284:    Level: intermediate

286: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
287: @*/
288: PetscErrorCode  PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
289: {
290:   PetscInt  i,mid,last,itmp,j,first;
291:   PetscReal d,tmp;
292:   PetscReal abskey;

295:   first = 0;
296:   last  = n-1;
297:   if (ncut < first || ncut > last) return(0);

299:   while (1) {
300:     mid    = first;
301:     d      = a[mid];
302:     abskey = PetscAbsReal(d);
303:     i      = last;
304:     for (j = first + 1; j <= i; ++j) {
305:       d = a[j];
306:       if (PetscAbsReal(d) >= abskey) {
307:         ++mid;
308:         /* interchange */
309:         tmp = a[mid];  itmp = idx[mid];
310:         a[mid] = a[j]; idx[mid] = idx[j];
311:         a[j] = tmp;    idx[j] = itmp;
312:       }
313:     }

315:     /* interchange */
316:     tmp = a[mid];      itmp = idx[mid];
317:     a[mid] = a[first]; idx[mid] = idx[first];
318:     a[first] = tmp;    idx[first] = itmp;

320:     /* test for while loop */
321:     if (mid == ncut) break;
322:     else if (mid > ncut) last = mid - 1;
323:     else first = mid + 1;
324:   }
325:   return(0);
326: }