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 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: }