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>                /*I  "petscsys.h"  I*/

 10: #define SWAP(a,b,t) {t=a;a=b;b=t;}
 11: 
 14: /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
 15: static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
 16: {
 17:   PetscInt  i,last;
 18:   PetscReal vl,tmp;
 19: 
 21:   if (right <= 1) {
 22:     if (right == 1) {
 23:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
 24:     }
 25:     return(0);
 26:   }
 27:   SWAP(v[0],v[right/2],tmp);
 28:   vl   = v[0];
 29:   last = 0;
 30:   for (i=1; i<=right; i++) {
 31:     if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
 32:   }
 33:   SWAP(v[0],v[last],tmp);
 34:   PetscSortReal_Private(v,last-1);
 35:   PetscSortReal_Private(v+last+1,right-(last+1));
 36:   return(0);
 37: }

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

 44:    Not Collective

 46:    Input Parameters:
 47: +  n  - number of values
 48: -  v  - array of doubles

 50:    Level: intermediate

 52:    Concepts: sorting^doubles

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

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

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

 83:    Not Collective

 85:    Input Parameters:
 86: +  ncut  - splitig index
 87: .  n     - number of values to sort
 88: .  a     - array of values
 89: -  idx   - index for array a

 91:    Output Parameters:
 92: +  a     - permuted array of values such that its elements satisfy:
 93:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 
 94:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 
 95: -  idx   - permuted index of array a

 97:    Level: intermediate

 99:    Concepts: sorting^doubles

101: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
102: @*/
103: PetscErrorCode  PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
104: {
105:   PetscInt    i,mid,last,itmp,j,first;
106:   PetscScalar d,tmp;
107:   PetscReal   abskey;

110:   first = 0;
111:   last  = n-1;
112:   if (ncut < first || ncut > last) return(0);

114:   while (1){
115:     mid = first;
116:     abskey = (d = a[mid],PetscAbsScalar(d));
117:     i = last;
118:     for (j = first + 1; j <= i; ++j) {
119:       if ((d = a[j],PetscAbsScalar(d)) >= abskey) {
120:         ++mid;
121:         /* interchange */
122:         tmp = a[mid];  itmp = idx[mid];
123:         a[mid] = a[j]; idx[mid] = idx[j];
124:         a[j] = tmp;    idx[j] = itmp;
125:       }
126:     }
127: 
128:     /* interchange */
129:     tmp = a[mid];      itmp = idx[mid];
130:     a[mid] = a[first]; idx[mid] = idx[first];
131:     a[first] = tmp;    idx[first] = itmp;

133:     /* test for while loop */
134:     if (mid == ncut) {
135:       break;
136:     } else if (mid > ncut){
137:       last = mid - 1;
138:     } else {
139:       first = mid + 1;
140:     }
141:   }
142:   return(0);
143: }

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

150:    Not Collective

152:    Input Parameters:
153: +  ncut  - splitig index
154: .  n     - number of values to sort
155: .  a     - array of values in PetscReal
156: -  idx   - index for array a

158:    Output Parameters:
159: +  a     - permuted array of real values such that its elements satisfy:
160:            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and 
161:            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut 
162: -  idx   - permuted index of array a

164:    Level: intermediate

166:    Concepts: sorting^doubles

168: .seealso: PetscSortInt(), PetscSortRealWithPermutation()
169: @*/
170: PetscErrorCode  PetscSortSplitReal(PetscInt ncut,PetscInt n,PetscReal a[],PetscInt idx[])
171: {
172:   PetscInt    i,mid,last,itmp,j,first;
173:   PetscReal   d,tmp;
174:   PetscReal   abskey;

177:   first = 0;
178:   last  = n-1;
179:   if (ncut < first || ncut > last) return(0);

181:   while (1){
182:     mid = first;
183:     abskey = (d = a[mid],PetscAbsReal(d));
184:     i = last;
185:     for (j = first + 1; j <= i; ++j) {
186:       if ((d = a[j],PetscAbsReal(d)) >= abskey) {
187:         ++mid;
188:         /* interchange */
189:         tmp = a[mid];  itmp = idx[mid];
190:         a[mid] = a[j]; idx[mid] = idx[j];
191:         a[j] = tmp;    idx[j] = itmp;
192:       }
193:     }
194: 
195:     /* interchange */
196:     tmp = a[mid];      itmp = idx[mid];
197:     a[mid] = a[first]; idx[mid] = idx[first];
198:     a[first] = tmp;    idx[first] = itmp;

200:     /* test for while loop */
201:     if (mid == ncut) {
202:       break;
203:     } else if (mid > ncut){
204:       last = mid - 1;
205:     } else {
206:       first = mid + 1;
207:     }
208:   }
209:   return(0);
210: }