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