Actual source code: sortd.c

petsc-master 2017-11-21
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:    Concepts: sorting^doubles

 51: .seealso: PetscSortInt(), PetscSortRealWithPermutation(), PetscSortRealWithArrayInt()
 52: @*/
 53: PetscErrorCode  PetscSortReal(PetscInt n,PetscReal v[])
 54: {
 55:   PetscInt  j,k;
 56:   PetscReal tmp,vk;

 60:   if (n<8) {
 61:     for (k=0; k<n; k++) {
 62:       vk = v[k];
 63:       for (j=k+1; j<n; j++) {
 64:         if (vk > v[j]) {
 65:           SWAP(v[k],v[j],tmp);
 66:           vk = v[k];
 67:         }
 68:       }
 69:     }
 70:   } else PetscSortReal_Private(v,n-1);
 71:   return(0);
 72: }

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

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

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

105:    Not Collective

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

112:    Level: intermediate

114:    Concepts: sorting^ints with array

116: .seealso: PetscSortReal()
117: @*/
118: PetscErrorCode  PetscSortRealWithArrayInt(PetscInt n,PetscReal r[],PetscInt Ii[])
119: {
121:   PetscInt       j,k,itmp;
122:   PetscReal      rk,rtmp;

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

143: /*@
144:   PetscFindReal - Finds a PetscReal in a sorted array of PetscReals

146:    Not Collective

148:    Input Parameters:
149: +  key - the value to locate
150: .  n   - number of values in the array
151: .  ii  - array of values
152: -  eps - tolerance used to compare

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

157:    Level: intermediate

159:    Concepts: sorting^ints

161: .seealso: PetscSortReal(), PetscSortRealWithArrayInt()
162: @*/
163: PetscErrorCode PetscFindReal(PetscReal key, PetscInt n, const PetscReal t[], PetscReal eps, PetscInt *loc)
164: {
165:   PetscInt lo = 0,hi = n;

169:   if (!n) {*loc = -1; return(0);}
171:   while (hi - lo > 1) {
172:     PetscInt mid = lo + (hi - lo)/2;
173:     if (key < t[mid]) hi = mid;
174:     else              lo = mid;
175:   }
176:   *loc = (PetscAbsReal(key - t[lo]) < eps) ? lo : -(lo + (key > t[lo]) + 1);
177:   return(0);
178: }

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

183:    Not Collective

185:    Input Parameters:
186: +  n  - number of values
187: -  v  - array of doubles

189:    Output Parameter:
190: .  n - number of non-redundant values

192:    Level: intermediate

194:    Concepts: sorting^doubles

196: .seealso: PetscSortReal(), PetscSortRemoveDupsInt()
197: @*/
198: PetscErrorCode  PetscSortRemoveDupsReal(PetscInt *n,PetscReal v[])
199: {
201:   PetscInt       i,s = 0,N = *n, b = 0;

204:   PetscSortReal(N,v);
205:   for (i=0; i<N-1; i++) {
206:     if (v[b+s+1] != v[b]) {
207:       v[b+1] = v[b+s+1]; b++;
208:     } else s++;
209:   }
210:   *n = N - s;
211:   return(0);
212: }

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

217:    Not Collective

219:    Input Parameters:
220: +  ncut  - splitig index
221: .  n     - number of values to sort
222: .  a     - array of values
223: -  idx   - index for array a

225:    Output Parameters:
226: +  a     - permuted array of values such that its elements satisfy:
227:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
228:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
229: -  idx   - permuted index of array a

231:    Level: intermediate

233:    Concepts: sorting^doubles

235: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
236: @*/
237: PetscErrorCode  PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
238: {
239:   PetscInt    i,mid,last,itmp,j,first;
240:   PetscScalar d,tmp;
241:   PetscReal   abskey;

244:   first = 0;
245:   last  = n-1;
246:   if (ncut < first || ncut > last) return(0);

248:   while (1) {
249:     mid    = first;
250:     d      = a[mid];
251:     abskey = PetscAbsScalar(d);
252:     i      = last;
253:     for (j = first + 1; j <= i; ++j) {
254:       d = a[j];
255:       if (PetscAbsScalar(d) >= abskey) {
256:         ++mid;
257:         /* interchange */
258:         tmp = a[mid];  itmp = idx[mid];
259:         a[mid] = a[j]; idx[mid] = idx[j];
260:         a[j] = tmp;    idx[j] = itmp;
261:       }
262:     }

264:     /* interchange */
265:     tmp = a[mid];      itmp = idx[mid];
266:     a[mid] = a[first]; idx[mid] = idx[first];
267:     a[first] = tmp;    idx[first] = itmp;

269:     /* test for while loop */
270:     if (mid == ncut) break;
271:     else if (mid > ncut) last = mid - 1;
272:     else first = mid + 1;
273:   }
274:   return(0);
275: }

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

280:    Not Collective

282:    Input Parameters:
283: +  ncut  - splitig index
284: .  n     - number of values to sort
285: .  a     - array of values in PetscReal
286: -  idx   - index for array a

288:    Output Parameters:
289: +  a     - permuted array of real values such that its elements satisfy:
290:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
291:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
292: -  idx   - permuted index of array a

294:    Level: intermediate

296:    Concepts: sorting^doubles

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

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

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

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

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