MIDAPACK - MIcrowave Data Analysis PACKage  1.1b
Parallel software tools for high performance CMB DA analysis
 All Data Structures Files Functions Variables Typedefs Groups Pages
csort.c
Go to the documentation of this file.
1 
8 #include <stdlib.h>
9 #include <string.h>
10 //int per page size
11 #define PAGE 1024
12 
18 void insertion_sort(int *indices, int count){
19  int i, j;
20  int tmp;
21  for(i=0; i<count-1 ; i++){
22  tmp = indices[i+1];
23  j=i;
24  while(j != -1 && tmp < indices[j]){
25  indices[j+1] = indices[j];
26  indices[j]=tmp;
27  j--;
28  }
29  }
30 }
31 
40 void quick_sort(int *indices, int left, int right){
41  int pivot;
42  int tmp, key;
43  int i,j;
44  if (left<right){
45  key=indices[left];
46  i=left+1;
47  j=right;
48  while(i<=j){
49  while((i<=right) && (indices[i]<=key)) i++;
50  while ((j>left) && (indices[j]>key)) j--;
51  if(i<j){
52  tmp=indices[i];
53  indices[i] = indices[j];
54  indices[j] = tmp;
55  i++;
56  j--;
57  }
58  }
59  tmp=indices[left];
60  indices[left] = indices[j];
61  indices[j] = tmp;
62  pivot = j;
63  quick_sort(indices, left, pivot-1);
64  quick_sort(indices, pivot+1, right);
65  }
66 }
67 
73 void bubble_sort(int *indices, int count){
74  int i, j, tmp;
75  for (i=(count-1); i>0; i--){
76  for (j = 1; j <= i; j++){
77  if (indices[j-1] > indices[j]){
78  tmp = indices[j-1];
79  indices[j-1] = indices[j];
80  indices[j] = tmp;
81  }
82  }
83  }
84 }
85 
95 int counting_sort(int *indices, int count){
96  int *buf;
97  int i, j, k;
98  int min, max;
99  min=indices[0];
100  max=indices[0];
101  for(i=1; i<count ; i++){
102  if(indices[i]>max){
103  max=indices[i];
104  }
105  else{
106  if(indices[i]<min)
107  min=indices[i];
108  }
109  }
110  buf = (int *) calloc(max-min+1, sizeof(int));
111  for(i=0; i<count ; i++){
112  buf[indices[i]-min]=1;
113  }
114  j=0;
115  for(i=0; i<(max-min+1) ; i++){
116  if(buf[i]==1){
117  indices[j]=min+i;
118  j++;
119  }
120  }
121  free(buf);
122  return j;
123 }
124 
129 void shell_sort(int *a,int n){
130  int j,i,k,m,mid;
131  for(m = n/2;m>0;m/=2){
132  for(j = m;j< n;j++){
133  for(i=j-m;i>=0;i-=m){
134  if(a[i+m]>=a[i])
135  break;
136  else{
137  mid = a[i];
138  a[i] = a[i+m];
139  a[i+m] = mid;
140  }
141  }
142  }
143  }
144 }
145 
146 
161 int ssort(int *indices, int count, int flag){
162  int i, n;
163  int *ptr_i, *ptr_o;
164  switch(flag){
165  case 0 :
166  quick_sort(indices, 0, count-1);
167  break;
168  case 1 :
169  bubble_sort(indices, count);
170  break;
171  case 2 :
172  insertion_sort(indices, count);
173  break;
174  case 3 :
175  n=counting_sort(indices, count);
176  return n;
177  case 4 :
178  shell_sort(indices, count);
179  break;
180  }
181  ptr_i = indices;
182  ptr_o = indices;
183  n=1;
184  for(i=0; i<count-1; i++){
185  ptr_i++;
186  if(*ptr_i != *ptr_o){
187  ptr_o++;
188  n++;
189  *ptr_o = *ptr_i;
190  }
191  }
192  return n;
193 }
194 
195 //optimized version is faster than the other implementation but there is a bug !!!
196 #ifdef W_OPENMP
197 int omp_psort_opt(int *A, int nA, int flag){
198  int i;
199  int *count, *disp;
200  int q, r;
201  int p, l;
202  int tid, nths;
203  int *buf;
204  int *ptr_i, *ptr_o;
205  int n, k, d;
206 
207  #pragma omp parallel shared(nths)
208  {//---fork---just to get the number of threads
209  nths = omp_get_num_threads();
210  }//---join---
211 
212  p = nA/PAGE; //number of full pages
213  q = p/nths; //full pages per thread
214  l = p%nths; //full pages left
215  r = nA%PAGE; //number of elements the last page not full
216 
217  count = (int *) malloc(nths *sizeof(int));
218  disp = (int *) malloc(nths *sizeof(int));
219 
220  for(i=0; i<nths; i++){
221  count[i] = q*PAGE;
222  if(i<l)
223  count[i] += PAGE;
224  if(i==l)
225  count[i] += r;
226  }
227 
228  disp[0] = 0;
229  for(i=0; i<nths-1; i++){
230  disp[i+1] = disp[i] + count[i];
231  }
232 
233  #pragma omp parallel private(tid, n, k, d, buf) shared(nths, A, nA, disp, count)
234  {//---fork---1st step, sort on local chunk
235  tid = omp_get_thread_num();
236 
237  buf = (int *) malloc(nA*sizeof(int));
238  //buf = (int *) malloc(count[tid]*sizeof(int));
239  memcpy(buf, A+disp[tid], count[tid]*sizeof(int));
240 
241  n = ssort(buf, count[tid], flag);
242  count[tid]=n;
243 
244  memcpy(A+disp[tid], buf, n*sizeof(int));
245 
246  nths = omp_get_num_threads();
247 
248 
249  k = nths;
250  d = 1;
251  while(k>1){
252  #pragma omp barrier
253  if(tid%(2*d)==0 && tid+d<nths){
254  set_or(A+disp[tid], count[tid] , A+disp[tid+d], count[tid+d], buf);
255  count[tid]=n;
256  memcpy(A+disp[tid], buf, n*sizeof(int));
257  d*=2;
258  k/=2;
259  }
260  }
261  free(buf);
262  }//---join---
263 
264  nA=count[0];
265 // printf("\nA :");
266 // for(i=0; i<nA; i++){
267 // printf(" %d",A[i]);
268 // }
269  free(count);
270  free(disp);
271  return nA;
272 }
273 
291 int omp_psort(int *A, int nA, int flag){
292  int i;
293  int *count, *disp;
294  int q, r;
295  int tid, nths;
296  int *buf;
297  int *ptr_i, *ptr_o;
298  int n, k, d;
299 
300  #pragma omp parallel private(tid) shared(nths)
301  {//---fork---just to get the number of threads
302  nths = omp_get_num_threads();
303  }//---join---
304 
305  q = nA/nths;
306  r = nA%nths;
307 
308  count = (int *) malloc(nths *sizeof(int));
309  disp = (int *) malloc(nths *sizeof(int));
310 
311  for(i=0; i<nths; i++){
312  if(i<r){
313  count[i] = q+1;
314  }
315  else{
316  count[i] = q;
317  }
318  }
319 
320  disp[0] = 0;
321  for(i=0; i<nths-1; i++){
322  disp[i+1] = disp[i] + count[i];
323  }
324 
325  #pragma omp parallel private(tid, n) shared(A, disp, count)
326  {//---fork---1st step, sort on local chunk
327  tid = omp_get_thread_num();
328  n = ssort(A+disp[tid], count[tid], flag);
329  count[tid]=n;
330  }//---join---
331 
332  buf = (int *) malloc(nA*sizeof(int));
333 
334  #pragma omp parallel private(tid, n, k, d) shared(nths, nA, A, disp, count, buf)
335  {//---fork---2nd step, gathering with a binary tree scheme
336  tid = omp_get_thread_num();
337  nths = omp_get_num_threads();
338  k = nths;
339  d = 1;
340  while(k>1){
341  if(tid%(2*d)==0 && tid+d<nths){
342 // printf("\nd %d, thread %d, count+ %d, disp %d",d , tid, count[tid], disp[tid]);
343  n = card_or(A+disp[tid], count[tid] , A+disp[tid+d], count[tid+d]);
344  set_or(A+disp[tid], count[tid] , A+disp[tid+d], count[tid+d], buf+disp[tid]);
345  count[tid]=n;
346  memcpy(A+disp[tid], buf+disp[tid], n*sizeof(int));
347  d*=2;
348  k/=2;
349  }
350  #pragma omp barrier
351  }
352  }//---join---
353 
354  nA=count[0];
355  free(buf);
356  free(count);
357  free(disp);
358  return nA;
359 }
360 #endif
361 
362 
364 int sorted(int *indices, int count){
365  int i=0;
366  while(i<count-2){
367  if(indices[i]>indices[i+1]){
368  return 1;
369  }
370  else{
371  i++;
372  }
373  }
374  return 0;
375 }
376 
377 int monotony(int *indices, int count){
378  int i=0;
379  while(i<count-2){
380  if(indices[i]>=indices[i+1]){
381  return 1;
382  }
383  else{
384  i++;
385  }
386  }
387  return 0;
388 }