fixed for VC++ windows
[charm.git] / src / ck-ldb / tm_bucket.c
1 #ifdef _WIN32
2 #include <windows.h>
3 #include <winbase.h>
4 #define random() rand()
5 #define srandom(x)  srand(x)
6 #endif
7 #include <stdio.h>
8 #include <float.h>
9 #include <math.h>
10 #include <assert.h>
11 #include "tm_tree.h"
12 #include "tm_bucket.h"
13 #include "tm_timings.h" 
14
15 #undef DEBUG
16 bucket_list_t global_bl; 
17
18 int tab_cmp(const void* x1,const void* x2){ 
19   int *e1,*e2,i1,i2,j1,j2;
20   double **tab;
21   bucket_list_t bl;
22
23   bl=global_bl;
24
25   e1=((int *)x1);
26   e2=((int *)x2);
27
28   
29   tab=bl->tab;
30
31   i1=e1[0];
32   j1=e1[1];
33   i2=e2[0];
34   j2=e2[1];
35
36   return tab[i1][j1]>tab[i2][j2]?-1:1;
37 }
38
39
40 int old_bucket_id(int i,int j,bucket_list_t bucket_list){
41   double *pivot,val;
42   int n,sup,inf,p;
43   pivot=bucket_list->pivot;
44   n=bucket_list->nb_buckets;  
45   val=bucket_list->tab[i][j];
46   
47   inf=-1;
48   sup=n;
49
50   while(sup-inf>1){
51     p=(sup+inf)/2;
52     //printf("%f [%d,%d,%d]=%f\n",val,inf,p,sup,pivot[p]);
53     if(val<pivot[p]){
54       inf=p;
55       if(inf==sup)
56         inf--;
57     }else{
58       sup=p;
59       if(sup==inf)
60         sup++;
61     }
62   }
63   //exit(-1);
64   return sup;
65 }
66
67
68 int bucket_id(int i,int j,bucket_list_t bucket_list){
69   double *pivot_tree,val;
70   int p,k;
71   pivot_tree=bucket_list->pivot_tree;
72   val=bucket_list->tab[i][j];
73
74
75
76   p=1;
77   for(k=0;k<bucket_list->max_depth;k++){
78     if(val>pivot_tree[p])
79       p=p*2;
80     else
81       p=p*2+1;
82   }
83
84   return  (int)pivot_tree[p];
85 }
86   
87
88
89 void  display_bucket(bucket_t *b){
90   printf("\tb.bucket=%p\n",b->bucket);
91   printf("\tb.bucket_len=%d\n",(int)b->bucket_len);
92   printf("\tb.nb_elem=%d\n",(int)b->nb_elem);
93
94 }
95
96 void check_bucket(bucket_t *b,double **tab,double inf, double sup,int N){
97   int k,i,j;
98   for(k=0;k<b->nb_elem;k++){
99     i=b->bucket[k].i;
100     j=b->bucket[k].j;
101     if((tab[i][j]<inf)||(tab[i][j]>sup)){
102       printf("[%d] (%d,%d):%f not in [%f,%f]\n",k,i,j,tab[i][j],inf,sup);
103       exit(-1);
104     }
105   }
106 }
107
108 void  display_bucket_list(bucket_list_t bucket_list){
109   int i;
110   double inf,sup;
111
112   for(i=0;i<bucket_list->nb_buckets-1;i++){
113     printf("pivot[%d]=%f\n",i,bucket_list->pivot[i]);
114   }
115   printf("\n");
116
117   for(i=0;i<bucket_list->nb_buckets;i++){
118     inf=bucket_list->pivot[i];
119     sup=bucket_list->pivot[i-1];
120     if(i==0)
121       sup=DBL_MAX;
122     if(i==bucket_list->nb_buckets-1)
123       inf=0;
124     printf("Bucket %d:\n",i);
125     display_bucket(bucket_list->bucket_tab[i]);
126     printf("\n");
127     check_bucket(bucket_list->bucket_tab[i],bucket_list->tab,inf,sup,bucket_list->N);
128   }
129   
130 }
131
132 void add_to_bucket(int id,int i,int j,bucket_list_t bucket_list){
133   bucket_t *bucket;
134   int N,n,size;
135
136   bucket=bucket_list->bucket_tab[id];
137   //display_bucket(bucket);
138   
139   if(bucket->bucket_len==bucket->nb_elem){
140     N=bucket_list->N;
141     n=bucket_list->nb_buckets;  
142     size=N*N/n;
143     //display_bucket(bucket);
144     bucket->bucket=(coord*)realloc(bucket->bucket,sizeof(coord)*(size+bucket->bucket_len));
145     bucket->bucket_len+=size;
146 #ifdef DEBUG
147     printf("malloc/realloc: %d\n",id);
148     printf("(%d,%d)\n",i,j);
149     display_bucket(bucket);
150     printf("\n");
151 #endif
152   }
153   
154  bucket->bucket[bucket->nb_elem].i=i;
155  bucket->bucket[bucket->nb_elem].j=j;
156  bucket->nb_elem++;
157
158   //printf("\n");
159   //exit(-1);
160 }
161
162 void dfs(int i,int inf,int sup,double *pivot,double *pivot_tree,int depth,int max_depth){
163   int p;
164   if(depth==max_depth)
165     return;
166
167   p=(inf+sup)/2;
168   pivot_tree[i]=pivot[p-1];
169
170   dfs(2*i,inf,p-1,pivot,pivot_tree,depth+1,max_depth);
171   dfs(2*i+1,p+1,sup,pivot,pivot_tree,depth+1,max_depth);
172 }
173
174 void  built_pivot_tree(bucket_list_t bucket_list){
175   double *pivot_tree,*pivot;
176   int n,i,k;
177   pivot=bucket_list->pivot;
178   n=bucket_list->nb_buckets;
179   pivot_tree=(double*)malloc(sizeof(double)*2*n);
180   bucket_list->max_depth=(int)log2(n);
181
182   dfs(1,1,n-1,pivot,pivot_tree,0,bucket_list->max_depth);
183
184   k=0;
185   for(i=n;i<2*n;i++)
186     pivot_tree[i]=k++;
187
188   bucket_list->pivot_tree=pivot_tree;  
189   /*
190   for(i=0;i<2*n;i++)
191     printf("%d:%f\t",i,pivot_tree[i]);
192   printf("\n");
193   */
194 }
195
196 void fill_buckets(bucket_list_t bucket_list){
197   int N,i,j,id;
198
199   N=bucket_list->N;
200
201   for(i=0;i<N;i++){
202     for(j=i+1;j<N;j++){
203       id=bucket_id(i,j,bucket_list);
204       add_to_bucket(id,i,j,bucket_list);
205     }
206   }
207   
208   
209 }
210
211 int is_power_of_2(int val){
212   int n=1;
213   do{
214     if(n==val)
215       return 1;
216     n<<=1;
217   }while(n>0);
218   return 0;
219 }
220
221
222 void partial_sort(bucket_list_t *bl,double **tab,int N,int nb_buckets){
223   int *sample;
224   int i,j,k,n;
225   double *pivot;
226   bucket_list_t bucket_list;
227
228
229   if(!is_power_of_2(nb_buckets)){
230     fprintf(stderr,"Error! Paramater nb_buckets is: %d and should be a power of 2\n",nb_buckets);
231     exit(-1);
232   }
233
234
235   bucket_list=(bucket_list_t)malloc(sizeof(_bucket_list_t));
236
237   bucket_list->tab=tab;
238   bucket_list->N=N;
239
240
241   n=pow(nb_buckets,2);
242   n=N;
243   sample=(int*)malloc(2*sizeof(int)*n);
244   
245   for(k=0;k<n;k++){
246     i=random()%(N-2)+1;
247     if(i==N-2)
248       j=N-1;
249     else
250       j=random()%(N-i-2)+i+1;
251     assert(i!=j);
252     assert(i<j);
253     assert(i<N);
254     assert(j<N);
255     sample[2*k]=i;
256     sample[2*k+1]=j;
257   }
258   
259   global_bl=bucket_list;
260   qsort(sample,n,2*sizeof(int),tab_cmp);
261
262   /*  for(k=0;k<N;k++){
263     i=sample[k]/N;
264     j=sample[k]%N;
265     printf("%f\n",tab[i][j]);
266   }
267   */
268   pivot=(double*)malloc(sizeof(double)*nb_buckets-1);
269   int id=2;
270   for(k=1;k<nb_buckets;k++){
271     i=sample[2*(id-1)];
272     j=sample[2*(id-1)+1];
273     id*=2;
274     
275
276     /*    i=sample[k*N/nb_buckets]/N;
277           j=sample[k*N/nb_buckets]%N;*/
278     pivot[k-1]=tab[i][j];
279     //printf("pivot[%d]=%f\n",k-1,tab[i][j]);
280   }
281
282   bucket_list->pivot=pivot;
283   bucket_list->nb_buckets=nb_buckets;
284   built_pivot_tree(bucket_list);
285   
286   bucket_list->bucket_tab=(bucket_t**)malloc(nb_buckets*sizeof(bucket_t*));
287   for(i=0;i<nb_buckets;i++){
288     bucket_list->bucket_tab[i]=(bucket_t*)calloc(1,sizeof(bucket_t));
289   }
290
291   fill_buckets(bucket_list);
292   
293   //display_bucket_list(bucket_list);
294
295   bucket_list->cur_bucket=0;
296   bucket_list->bucket_indice=0;
297   
298   free(sample);
299
300   *bl=bucket_list;
301 }
302
303 void next_bucket_elem(bucket_list_t bucket_list,int *i,int *j){
304   int N;
305   bucket_t *bucket=bucket_list->bucket_tab[bucket_list->cur_bucket];
306
307
308   while(bucket->nb_elem<=bucket_list->bucket_indice){
309     bucket_list->bucket_indice=0;
310     bucket=bucket_list->bucket_tab[bucket_list->cur_bucket++];
311
312     printf("### From bucket %d to bucket %d\n",bucket_list->cur_bucket-1,bucket_list->cur_bucket);
313   }
314
315   if(!bucket->sorted){
316     global_bl=bucket_list;
317     qsort(bucket->bucket,bucket->nb_elem,2*sizeof(int),tab_cmp);
318     bucket->sorted=1;
319   }
320
321
322   
323   N=bucket_list->N;
324
325   *i=bucket->bucket[bucket_list->bucket_indice].i;
326   *j=bucket->bucket[bucket_list->bucket_indice].j;
327   bucket_list->bucket_indice++;
328 }
329
330
331 int add_edge_3(double **tab,tree_t *tab_node, tree_t *parent,int i,int j,int N,int *nb_groups){
332   //printf("%d <-> %d ?\n",tab_node[i].id,tab_node[j].id);
333
334   if((!tab_node[i].parent) && (!tab_node[j].parent)){
335     if(parent){
336       parent->child[0]=&tab_node[i];
337       parent->child[1]=&tab_node[j];
338       tab_node[i].parent=parent;
339       tab_node[j].parent=parent;
340 #ifdef DEBUG
341       printf("%d: %d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id);
342 #endif
343       return 1;
344     } 
345     return 0;
346   }
347   
348   if(tab_node[i].parent && (!tab_node[j].parent)){
349     parent=tab_node[i].parent;
350     if(!parent->child[2]){
351       parent->child[2]=&tab_node[j];
352       tab_node[j].parent=parent;
353 #ifdef DEBUG
354       printf("%d: %d-%d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id,parent->child[2]->id);
355 #endif
356       (*nb_groups)++;
357     }
358     return 0;
359   }
360
361   if(tab_node[j].parent && (!tab_node[i].parent)){
362     parent=tab_node[j].parent;
363     if(!parent->child[2]){
364       parent->child[2]=&tab_node[i];
365       tab_node[i].parent=parent;
366 #ifdef DEBUG
367       printf("%d: %d-%d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id,parent->child[2]->id);
368 #endif      
369       (*nb_groups)++;
370     }
371     return 0;
372   }
373
374   return 0;
375 }
376
377 int try_add_edge(double **tab,tree_t *tab_node, tree_t *parent,int arity,int i,int j,int N,int *nb_groups){
378
379   assert(i!=j);
380
381   
382   switch(arity){
383   case 2:
384     if(tab_node[i].parent)
385       return 0;
386     if(tab_node[j].parent)
387       return 0;
388
389     parent->child[0]=&tab_node[i];
390     parent->child[1]=&tab_node[j];
391     tab_node[i].parent=parent;
392     tab_node[j].parent=parent;
393     
394     (*nb_groups)++;
395
396     return 1;
397   case 3:
398     return add_edge_3(tab,tab_node,parent,i,j,N,nb_groups);
399   default:
400     fprintf(stderr,"Cannot handle arity %d\n",parent->arity);
401     exit(-1);
402   }
403 }
404
405
406 void free_bucket(bucket_t *bucket){
407   free(bucket->bucket);
408   free(bucket);
409
410 }
411
412 void free_tab_bucket(bucket_t **bucket_tab,int N){
413   int i;
414   for(i=0;i<N;i++){
415     free_bucket(bucket_tab[i]);
416   }
417   free(bucket_tab);
418 }
419
420
421 void free_bucket_list(bucket_list_t bucket_list){
422
423   // Do not free the tab field it is used elsewhere
424
425   free_tab_bucket(bucket_list->bucket_tab,bucket_list->nb_buckets);
426   free(bucket_list->pivot);
427   free(bucket_list->pivot_tree);
428   free(bucket_list);
429 }
430
431 void bucket_grouping(double **tab,tree_t *tab_node, tree_t *new_tab_node, int arity,int N, int M,long int k){
432   bucket_list_t bucket_list;
433   double duration;
434   int l,i,j,nb_groups;
435   double val=0;
436  
437   TIC;
438   partial_sort(&bucket_list,tab,N,8);
439   duration=TOC;
440   printf("Partial sorting=%fs\n",duration);  
441  
442   TIC;
443   l=0;
444   i=0;
445   nb_groups=0;
446   while(l<M){
447     next_bucket_elem(bucket_list,&i,&j);
448     if(try_add_edge(tab,tab_node,&new_tab_node[l],arity,i,j,N,&nb_groups)){
449       l++;
450     }
451   }
452
453 #ifdef DEBUG
454   printf("l=%d,nb_groups=%d\n",l,nb_groups);
455 #endif
456
457   while(nb_groups<M){
458     next_bucket_elem(bucket_list,&i,&j);
459     try_add_edge(tab,tab_node,NULL,arity,i,j,N,&nb_groups);
460   }
461
462 #ifdef DEBUG
463   printf("l=%d,nb_groups=%d\n",l,nb_groups);
464 #endif
465
466   for(l=0;l<M;l++){
467     update_val(tab,&new_tab_node[l],N);      
468     val+=new_tab_node[l].val;
469   }
470
471
472       
473
474
475   duration=TOC;
476   printf("Grouping =%fs\n",duration);  
477
478   printf("Bucket: %d, indice:%d\n",bucket_list->cur_bucket,bucket_list->bucket_indice);
479
480   printf("val=%f\n",val);
481   free_bucket_list(bucket_list);
482
483   //  exit(-1);
484
485   //  display_grouping(new_tab_node,M,arity,val);
486   
487
488 }
489