xref: /petsc/src/sys/utils/sortd.c (revision d97c5584da8f12a1e3bac82946f8ccd8395196fa)
1e5c89e4eSSatish Balay #define PETSC_DLL
2e5c89e4eSSatish Balay /*
3e5c89e4eSSatish Balay    This file contains routines for sorting doubles.  Values are sorted in place.
4e5c89e4eSSatish Balay    These are provided because the general sort routines incur a great deal
5e5c89e4eSSatish Balay    of overhead in calling the comparision routines.
6e5c89e4eSSatish Balay 
7e5c89e4eSSatish Balay  */
8e5c89e4eSSatish Balay #include "petsc.h"                /*I  "petsc.h"  I*/
9e5c89e4eSSatish Balay #include "petscsys.h"             /*I  "petscsys.h"    I*/
10e5c89e4eSSatish Balay 
11e5c89e4eSSatish Balay #define SWAP(a,b,t) {t=a;a=b;b=t;}
12e5c89e4eSSatish Balay 
13e5c89e4eSSatish Balay #undef __FUNCT__
14e5c89e4eSSatish Balay #define __FUNCT__ "PetscSortReal_Private"
15e5c89e4eSSatish Balay /* A simple version of quicksort; taken from Kernighan and Ritchie, page 87 */
16e5c89e4eSSatish Balay static PetscErrorCode PetscSortReal_Private(PetscReal *v,PetscInt right)
17e5c89e4eSSatish Balay {
18e5c89e4eSSatish Balay   PetscInt  i,last;
19e5c89e4eSSatish Balay   PetscReal vl,tmp;
20e5c89e4eSSatish Balay 
21e5c89e4eSSatish Balay   PetscFunctionBegin;
22e5c89e4eSSatish Balay   if (right <= 1) {
23e5c89e4eSSatish Balay     if (right == 1) {
24e5c89e4eSSatish Balay       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
25e5c89e4eSSatish Balay     }
26e5c89e4eSSatish Balay     PetscFunctionReturn(0);
27e5c89e4eSSatish Balay   }
28e5c89e4eSSatish Balay   SWAP(v[0],v[right/2],tmp);
29e5c89e4eSSatish Balay   vl   = v[0];
30e5c89e4eSSatish Balay   last = 0;
31e5c89e4eSSatish Balay   for (i=1; i<=right; i++) {
32e5c89e4eSSatish Balay     if (v[i] < vl) {last++; SWAP(v[last],v[i],tmp);}
33e5c89e4eSSatish Balay   }
34e5c89e4eSSatish Balay   SWAP(v[0],v[last],tmp);
35e5c89e4eSSatish Balay   PetscSortReal_Private(v,last-1);
36e5c89e4eSSatish Balay   PetscSortReal_Private(v+last+1,right-(last+1));
37e5c89e4eSSatish Balay   PetscFunctionReturn(0);
38e5c89e4eSSatish Balay }
39e5c89e4eSSatish Balay 
40e5c89e4eSSatish Balay #undef __FUNCT__
41e5c89e4eSSatish Balay #define __FUNCT__ "PetscSortReal"
42e5c89e4eSSatish Balay /*@
43e5c89e4eSSatish Balay    PetscSortReal - Sorts an array of doubles in place in increasing order.
44e5c89e4eSSatish Balay 
45e5c89e4eSSatish Balay    Not Collective
46e5c89e4eSSatish Balay 
47e5c89e4eSSatish Balay    Input Parameters:
48e5c89e4eSSatish Balay +  n  - number of values
49e5c89e4eSSatish Balay -  v  - array of doubles
50e5c89e4eSSatish Balay 
51e5c89e4eSSatish Balay    Level: intermediate
52e5c89e4eSSatish Balay 
53e5c89e4eSSatish Balay    Concepts: sorting^doubles
54e5c89e4eSSatish Balay 
55e5c89e4eSSatish Balay .seealso: PetscSortInt(), PetscSortRealWithPermutation()
56e5c89e4eSSatish Balay @*/
57e5c89e4eSSatish Balay PetscErrorCode PETSC_DLLEXPORT PetscSortReal(PetscInt n,PetscReal v[])
58e5c89e4eSSatish Balay {
59e5c89e4eSSatish Balay   PetscInt  j,k;
60e5c89e4eSSatish Balay   PetscReal tmp,vk;
61e5c89e4eSSatish Balay 
62e5c89e4eSSatish Balay   PetscFunctionBegin;
63e5c89e4eSSatish Balay   if (n<8) {
64e5c89e4eSSatish Balay     for (k=0; k<n; k++) {
65e5c89e4eSSatish Balay       vk = v[k];
66e5c89e4eSSatish Balay       for (j=k+1; j<n; j++) {
67e5c89e4eSSatish Balay 	if (vk > v[j]) {
68e5c89e4eSSatish Balay 	  SWAP(v[k],v[j],tmp);
69e5c89e4eSSatish Balay 	  vk = v[k];
70e5c89e4eSSatish Balay 	}
71e5c89e4eSSatish Balay       }
72e5c89e4eSSatish Balay     }
73e5c89e4eSSatish Balay   } else {
74e5c89e4eSSatish Balay     PetscSortReal_Private(v,n-1);
75e5c89e4eSSatish Balay   }
76e5c89e4eSSatish Balay   PetscFunctionReturn(0);
77e5c89e4eSSatish Balay }
78e5c89e4eSSatish Balay 
79*d97c5584SHong Zhang 
80*d97c5584SHong Zhang 
81*d97c5584SHong Zhang #undef __FUNCT__
82*d97c5584SHong Zhang #define __FUNCT__ "PetscSortSplit"
83*d97c5584SHong Zhang /*@
84*d97c5584SHong Zhang    PetscSortSplit - Quick-sort split of an array in place.
85*d97c5584SHong Zhang                     Modified from SPARSEKIT2qsplit()
86*d97c5584SHong Zhang 
87*d97c5584SHong Zhang    Not Collective
88*d97c5584SHong Zhang 
89*d97c5584SHong Zhang    Input Parameters:
90*d97c5584SHong Zhang +  ncut  - splitig index
91*d97c5584SHong Zhang .  n     - number of values to sort
92*d97c5584SHong Zhang .  a     - array of values
93*d97c5584SHong Zhang -  idx   - index for array a
94*d97c5584SHong Zhang 
95*d97c5584SHong Zhang    Output Parameters:
96*d97c5584SHong Zhang +  a     - permuted array of values such that its elements satisfy:
97*d97c5584SHong Zhang            abs(a[i]) >= abs(a[ncut-1]) for i < ncut and
98*d97c5584SHong Zhang            abs(a[i]) <= abs(a[ncut-1]) for i >= ncut
99*d97c5584SHong Zhang -  idx   - permuted index of array a
100*d97c5584SHong Zhang 
101*d97c5584SHong Zhang    Level: intermediate
102*d97c5584SHong Zhang 
103*d97c5584SHong Zhang    Concepts: sorting^doubles
104*d97c5584SHong Zhang 
105*d97c5584SHong Zhang .seealso: PetscSortInt(), PetscSortRealWithPermutation()
106*d97c5584SHong Zhang @*/
107*d97c5584SHong Zhang PetscErrorCode PETSC_DLLEXPORT PetscSortSplit(PetscInt ncut,PetscInt n,PetscScalar a[],PetscInt idx[])
108*d97c5584SHong Zhang {
109*d97c5584SHong Zhang   PetscInt    i,mid,last,itmp,j,first;
110*d97c5584SHong Zhang   PetscScalar d,tmp;
111*d97c5584SHong Zhang   PetscReal   abskey;
112*d97c5584SHong Zhang 
113*d97c5584SHong Zhang   PetscFunctionBegin;
114*d97c5584SHong Zhang   first = 0;
115*d97c5584SHong Zhang   last  = n-1;
116*d97c5584SHong Zhang   if (ncut < first || ncut > last) PetscFunctionReturn(0);
117*d97c5584SHong Zhang 
118*d97c5584SHong Zhang   while (1){
119*d97c5584SHong Zhang     mid = first;
120*d97c5584SHong Zhang     abskey = (d = a[mid],PetscAbsScalar(d));
121*d97c5584SHong Zhang     i = last;
122*d97c5584SHong Zhang     for (j = first + 1; j <= i; ++j) {
123*d97c5584SHong Zhang       if ((d = a[j],PetscAbsScalar(d)) >= abskey) {
124*d97c5584SHong Zhang         ++mid;
125*d97c5584SHong Zhang         /* interchange */
126*d97c5584SHong Zhang         tmp = a[mid];  itmp = idx[mid];
127*d97c5584SHong Zhang         a[mid] = a[j]; idx[mid] = idx[j];
128*d97c5584SHong Zhang         a[j] = tmp;    idx[j] = itmp;
129*d97c5584SHong Zhang       }
130*d97c5584SHong Zhang     }
131*d97c5584SHong Zhang 
132*d97c5584SHong Zhang     /* interchange */
133*d97c5584SHong Zhang     tmp = a[mid];      itmp = idx[mid];
134*d97c5584SHong Zhang     a[mid] = a[first]; idx[mid] = idx[first];
135*d97c5584SHong Zhang     a[first] = tmp;    idx[first] = itmp;
136*d97c5584SHong Zhang 
137*d97c5584SHong Zhang     /* test for while loop */
138*d97c5584SHong Zhang     if (mid == ncut) {
139*d97c5584SHong Zhang       break;
140*d97c5584SHong Zhang     } else if (mid > ncut){
141*d97c5584SHong Zhang       last = mid - 1;
142*d97c5584SHong Zhang     } else {
143*d97c5584SHong Zhang       first = mid + 1;
144*d97c5584SHong Zhang     }
145*d97c5584SHong Zhang   }
146*d97c5584SHong Zhang   PetscFunctionReturn(0);
147*d97c5584SHong Zhang }
148