Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Corrected a major bugus.
[simgrid.git] / src / surf / sdp.c
1 /*      $Id$     */
2
3 /* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved.     */
4
5 /* This program is free software; you can redistribute it and/or modify it
6  * under the terms of the license (GNU LGPL) which comes with this package. */
7
8
9 #include "xbt/log.h"
10 #include "xbt/sysdep.h"
11 #include "xbt/mallocator.h"
12 #include "maxmin_private.h"
13
14 #include <stdlib.h>
15 #ifndef MATH
16 #include <math.h>
17 #endif
18
19 /*
20  * SDP specific variables.
21  */
22 #include <declarations.h>
23
24
25 static void create_cross_link(struct constraintmatrix *myconstraints, int k);
26
27 static void addentry(struct constraintmatrix *constraints, 
28               struct blockmatrix *, int matno, int blkno, int indexi, int indexj, double ent, int blocksize);
29
30
31 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
32                                 "Logging specific to SURF (sdp)");
33 /*
34 ########################################################################
35 ######################## Simple Proportionnal fairness #################
36 ########################################################################
37 ####### Original problem ##########
38 #
39 #  Max x_1*x_2*...*x_n
40 #    (A.x)_1 <= b_1
41 #    (A.x)_2 <= b_2
42 #     ...
43 #    (A.x)_m <= b_m
44 #    x>=0
45 #
46 #  We assume in the following that n=2^K
47 #
48 ####### Standard SDP form #########
49 #
50 # A SDP can be written in the following standard form:
51
52 #   (P) min c1*x1+c2*x2+...+cm*xm
53 #       st  F1*x1+F2*x2+...+Fm*xm-F0=X
54 #                                 X >= 0
55 #
56 # Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
57 # constraint X>0 means that X must be positive semidefinite positive.
58 #
59 ########## Equivalent SDP #########
60 #
61 #  Variables:
62 #
63 #    y(k,i) for k in [0,K] and i in [1,2^k]
64 #
65 #  Structure :
66 #                             y(0,1)
67 #                y(1,1)                      y(1,2)
68 #        y(2,1)        y(2,2)        y(2,3)        y(2,4)
69 #    y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
70 #    -------------------------------------------------------
71 #     x_1     x_2    x_3    x_4    x_5    x_6    x_7    x_8
72 #  
73 #
74 #  Structure Constraints:
75 #
76 #    forall k in [0,K-1], and i in [1,2^k]
77 #      [ y(k,  2i-1)   y(k-1, i) ]
78 #      [ y(k-1, i  )   y(k,  2i) ]
79 #
80 #    2^K-1
81 #
82 #  Capacity Constraints:
83 #     forall l in [1,m]
84 #        -(A.y(K,*))_l >= - b_l
85 #
86 #  Positivity Constraints:
87 #    forall k in [0,K], and i in [1,2^k]
88 #        y(k,i) >= 0
89 #
90 #  Minimize -y(0,1)
91 */
92
93 //typedef struct lmm_system {
94 //  int modified;
95 //  s_xbt_swag_t variable_set;  /* a list of lmm_variable_t */
96 //  s_xbt_swag_t constraint_set;        /* a list of lmm_constraint_t */
97 //  s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */
98 //  s_xbt_swag_t saturated_variable_set;        /* a list of lmm_variable_t */
99 //  s_xbt_swag_t saturated_constraint_set;      /* a list of lmm_constraint_t_t */
100 //  xbt_mallocator_t variable_mallocator;
101 //} s_lmm_system_t;
102
103 #define get_y(a,b) (pow(2,a)+b-1)
104
105 void sdp_solve(lmm_system_t sys)
106 {
107   lmm_variable_t var = NULL;
108
109   /*
110    * SDP mapping variables.
111    */
112   int K=0;
113   int nb_var = 0;
114   int nb_cnsts = 0;
115   int flows = 0;
116   int links = 0;
117   int nb_cnsts_capacity=0;
118   int nb_cnsts_struct=0;
119   int nb_cnsts_positivy=0;
120   int block_num=0;
121   int block_size; 
122   int total_block_size=0;
123   int *isdiag=NULL;
124   FILE *sdpout = fopen("SDPA-printf.tmp","w");
125   int blocksz = 0;
126   double *tempdiag = NULL;
127   int matno=0;
128   int i=0;
129   int j=0;
130   int k=0;
131   
132   /* 
133    * CSDP library specific variables.
134    */
135   struct blockmatrix C;
136   struct blockmatrix X,Z;
137   double *y;
138   double pobj, dobj;
139   double *a;
140   struct constraintmatrix *constraints;  
141  
142   /*
143    * Classic maxmin variables.
144    */
145   lmm_constraint_t cnst = NULL;
146   lmm_element_t elem = NULL;
147   xbt_swag_t cnst_list = NULL;
148   xbt_swag_t var_list = NULL;
149   xbt_swag_t elem_list = NULL;
150
151   if ( !(sys->modified))
152     return;
153
154   /* 
155    * Initialize the var list variable with only the active variables. 
156    * Associate an index in the swag variables.
157    */
158   i = 1;
159   var_list = &(sys->variable_set);
160   DEBUG1("Variable set : %d", xbt_swag_size(var_list));
161
162   xbt_swag_foreach(var, var_list) {
163     var->value = 0.0;
164     if(var->weight) var->index = i++;
165   }
166
167   cnst_list=&(sys->saturated_constraint_set);  
168   DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
169
170   /*
171    * Those fields are the top level description of the platform furnished in the xml file.
172    */
173   flows = i-1;
174   links = xbt_swag_size(&(sys->active_constraint_set));
175
176   /* 
177    * This  number is found based on the tree structure explained on top.
178    */
179   double tmp_k;
180
181   tmp_k = (double) log((double)flows)/log(2.0);
182   K = (int) ceil(tmp_k);
183
184
185   /* 
186    * The number of variables in the SDP program. 
187    */
188   nb_var = get_y(K, pow(2,K));
189   fprintf(sdpout,"%d\n", nb_var);
190  
191   /*
192    * Find the size of each group of constraints.
193    */
194   nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
195   nb_cnsts_struct = (int)pow(2,K) - 1;
196   nb_cnsts_positivy = (int)pow(2,K);
197
198   /*
199    * The total number of constraints.
200    */
201   nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
202   fprintf(sdpout,"%d\n", nb_cnsts);
203
204   /*
205    * Keep track of which blocks have off diagonal entries. 
206    */
207   isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
208   for (i=1; i<=nb_cnsts; i++)
209     isdiag[i]=1;
210
211   C.nblocks = nb_cnsts;
212   C.blocks  = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
213   constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
214
215   for(i = 1; i <= nb_var; i++){
216     constraints[i].blocks=NULL; 
217   }
218   
219   a = (double *)calloc(nb_var+1, sizeof(double)); 
220
221   /*
222    * Block sizes.
223    */
224   block_num=1;
225   block_size=0;
226
227   /*
228    * For each constraint block do.
229    */
230   for(i = 1; i <= nb_cnsts; i++){
231     
232     /*
233      * Structured blocks are size 2 and all others are size 1.
234      */
235     if(i <= nb_cnsts_struct){
236       total_block_size += block_size = 2;
237       fprintf(sdpout,"2 ");
238     }else{
239       total_block_size += block_size = 1;
240       fprintf(sdpout,"1 ");
241     }
242     
243     /*
244      * All blocks are matrices.
245      */
246     C.blocks[block_num].blockcategory = MATRIX;
247     C.blocks[block_num].blocksize = block_size;
248     C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
249  
250     block_num++;
251   }
252
253   fprintf(sdpout,"\n");
254
255
256   /*
257    * Creates de objective function array.
258    */
259   a = (double *)calloc((nb_var+1), sizeof(double));
260   
261   for(i = 1; i <= nb_var; i++){
262     if(get_y(0,1)==i){
263       fprintf(sdpout,"-1 ");
264       a[i]=-1;
265     }else{
266       fprintf(sdpout,"0 ");
267       a[i]=0;
268     }
269   }
270   fprintf(sdpout,"\n");
271
272
273   /*
274    * Structure contraint blocks.
275    */
276   block_num = 1;
277   matno = 1;
278   for(k = 1; k <= K; k++){
279     for(i = 1; i <= pow(2,k-1); i++){
280       matno=get_y(k,2*i-1);
281       fprintf(sdpout,"%d %d 1 1 1\n", matno  , block_num);
282       addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
283
284       matno=get_y(k,2*i);
285       fprintf(sdpout,"%d %d 2 2 1\n", matno  , block_num);
286       addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
287
288       matno=get_y(k-1,i);
289       fprintf(sdpout,"%d %d 1 2 1\n", matno  , block_num);
290       addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);      
291       
292       matno=get_y(k-1,i);
293       fprintf(sdpout,"%d %d 2 1 1\n", matno  , block_num);
294       addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
295       
296       isdiag[block_num] = 0;
297       block_num++;
298     }
299   }
300
301  
302   /*
303    * Capacity constraint block.
304    */
305   xbt_swag_foreach(cnst, cnst_list) {
306
307     fprintf(sdpout,"0 %d 1 1 %d\n", block_num,  (int) - (cnst->bound));    
308     addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);    
309
310
311     elem_list = &(cnst->element_set);
312     xbt_swag_foreach(elem, elem_list) {
313       if(elem->variable->weight <=0) break;
314       fprintf(sdpout,"%d %d 1 1 %d\n", elem->variable->index, block_num, (int) - (elem->variable->value)); 
315       addentry(constraints, &C, elem->variable->index, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
316  
317     }
318     block_num++;
319   }
320
321   
322   /*
323    * Positivy constraint blocks.
324    */
325   for(i = 1; i <= pow(2,K); i++){
326     matno=get_y(K, i);
327     fprintf(sdpout,"%d %d 1 1 1\n", matno, block_num);
328     addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
329     block_num++;
330   }
331
332
333   /*
334    * At this point, we'll stop to recognize whether any of the blocks
335    * are "hidden LP blocks"  and correct the block type if needed.
336    */
337   for (i=1; i<=nb_cnsts; i++){
338     if ((C.blocks[i].blockcategory != DIAG) && 
339         (isdiag[i]==1) && (C.blocks[i].blocksize > 1)){
340       /*
341        * We have a hidden diagonal block!
342        */
343       
344       printf("Block %d is actually diagonal.\n",i);
345       
346       blocksz=C.blocks[i].blocksize;
347       tempdiag=(double *)calloc((blocksz+1), sizeof(double));
348       for (j=1; j<=blocksz; j++)
349         tempdiag[j]=C.blocks[i].data.mat[ijtok(j,j,blocksz)];
350       free(C.blocks[i].data.mat);
351       C.blocks[i].data.vec=tempdiag;
352       C.blocks[i].blockcategory=DIAG;
353     };
354   };
355   
356   
357   /*
358    * Next, setup issparse and NULL out all nextbyblock pointers.
359    */
360   struct sparseblock *p=NULL;
361   for (i=1; i<=k; i++) {
362     p=constraints[i].blocks;
363     while (p != NULL){
364           /*
365            * First, set issparse.
366            */
367           if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15)){
368             p->issparse=0;
369           }else{
370             p->issparse=1;
371           };
372           
373           if (C.blocks[p->blocknum].blockcategory == DIAG)
374             p->issparse=1;
375           
376           /*
377            * Setup the cross links.
378            */
379           
380           p->nextbyblock=NULL;
381           p=p->next;
382     };
383   };
384
385
386   /*
387    * Create cross link reference.
388    */
389   create_cross_link(constraints, nb_var);
390   
391  
392   /*
393    * Debuging print problem in SDPA format.
394    */
395   printf("Printing SDPA...\n");
396   if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
397     char *tmp=strdup("SDPA.tmp");
398     write_prob(tmp, nb_cnsts, nb_var, C, a, constraints);  
399     //int write_prob(char *fname, int n, int k, struct blockmatrix C, double *a, struct constraintmatrix *constraints);
400     free(tmp);
401   }
402
403   /*
404    * Initialize parameters.
405    */
406   printf("Initializing solution...\n");
407   initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);  
408   
409
410
411   /*
412    * Call the solver.
413    */
414   printf("Calling the solver...\n");
415   int ret = easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
416
417   switch(ret){
418   case 0:
419   case 1: printf("SUCCESS The problem is primal infeasible\n");
420           break;
421
422   case 2: printf("SUCCESS The problem is dual infeasible\n");
423           break;
424  
425   case 3: printf("Partial SUCCESS A solution has been found, but full accuracy was not achieved. One or more of primal infeasibility, dual infeasibility, or relative duality gap are larger than their tolerances, but by a factor of less than 1000.\n");
426           break;
427
428   case 4: printf("Failure. Maximum number of iterations reached.");
429           break;
430
431   case 5: printf("Failure. Stuck at edge of primal feasibility.");
432           break;
433
434   }
435
436
437   /*
438    * Write out the solution if necessary.
439    */
440   xbt_swag_foreach(cnst, cnst_list) {
441
442     elem_list = &(cnst->element_set);
443     xbt_swag_foreach(elem, elem_list) {
444       if(elem->variable->weight <=0) break;
445       
446       i = (int)get_y(K, elem->variable->index);
447       elem->variable->value = y[i]; 
448       
449     }
450   }
451
452
453   /*
454    * Free up memory.
455    */
456   free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
457
458   fclose(sdpout);
459   free(isdiag);
460   sys->modified = 0;
461
462   if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
463     lmm_print(sys);
464   }
465
466 }
467
468
469 /*
470  * Create the cross_link reference in order to have a byblock list.
471  */
472 void create_cross_link(struct constraintmatrix *myconstraints, int k){
473
474   int i, j;
475   int blk;
476   struct sparseblock *p;
477   struct sparseblock *q;
478
479   struct sparseblock *prev;
480
481   /*
482    * Now, cross link.
483    */
484   prev=NULL;
485   for (i=1; i<=k; i++){
486     p=myconstraints[i].blocks;
487     while (p != NULL){
488       if (p->nextbyblock == NULL){
489         blk=p->blocknum;
490         
491         /*
492          * link in the remaining blocks.
493          */
494         for (j=i+1; j<=k; j++){
495           q=myconstraints[j].blocks;
496                   
497           while (q != NULL){
498             if (q->blocknum == p->blocknum){
499               if (p->nextbyblock == NULL){
500                 p->nextbyblock=q;
501                 q->nextbyblock=NULL;
502                 prev=q;
503               }
504               else{
505                 prev->nextbyblock=q;
506                 q->nextbyblock=NULL;
507                 prev=q;
508               }
509               break;
510             }
511             q=q->next;
512           }
513         }
514       }
515       p=p->next;
516     }
517   }
518 }
519
520
521
522  
523 void addentry(struct constraintmatrix *constraints, 
524               struct blockmatrix *C, 
525               int matno,
526               int blkno,
527               int indexi,
528               int indexj,
529               double ent, 
530               int blocksize)
531 {
532   struct sparseblock *p;
533   struct sparseblock *p_sav;
534
535   p=constraints[matno].blocks;
536   
537   if (matno != 0.0) {
538     if (p == NULL){
539       /*
540        * We haven't yet allocated any blocks.
541        */
542       p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
543       
544       //two entries because this library ignores indices starting in zerox
545       p->constraintnum=matno;
546       p->blocknum=blkno;
547       p->numentries=1;
548       p->next=NULL;
549
550       p->entries=calloc(p->numentries+1, sizeof(double));
551       p->iindices=calloc(p->numentries+1, sizeof(int));
552       p->jindices=calloc(p->numentries+1, sizeof(int));
553    
554       p->entries[p->numentries]=ent;
555       p->iindices[p->numentries]=indexi;
556       p->jindices[p->numentries]=indexj;
557
558       p->blocksize=blocksize;
559
560       constraints[matno].blocks=p;
561     } else {
562       /*
563        * We have some existing blocks.  See whether this block is already
564        * in the chain.
565        */
566       while (p != NULL){
567         if (p->blocknum == blkno){
568           /*
569            * Found the right block.
570            */
571           p->constraintnum=matno;
572           p->blocknum=blkno;
573           p->numentries=p->numentries+1;
574
575           p->entries = realloc(p->entries,  (p->numentries+1) * sizeof(double) );
576           p->iindices = realloc(p->iindices, (p->numentries+1) * sizeof(int) );
577           p->jindices = realloc(p->jindices, (p->numentries+1) * sizeof(int) );
578    
579           p->entries[p->numentries]=ent;
580           p->iindices[p->numentries]=indexi;
581           p->jindices[p->numentries]=indexj;
582         
583           return;
584         }
585         p_sav=p;
586         p=p->next;
587       }
588
589       /*
590        * If we get here, we have a non-empty structure but not the right block
591        * inside hence create a new structure.
592        */
593       
594       p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
595         
596       //two entries because this library ignores indices starting in zerox
597       p->constraintnum=matno;
598       p->blocknum=blkno;
599       p->numentries=1;
600       p->next=NULL;
601
602       p->entries=calloc(p->numentries+1, sizeof(double));
603       p->iindices=calloc(p->numentries+1, sizeof(int));
604       p->jindices=calloc(p->numentries+1, sizeof(int));
605    
606       p->entries[p->numentries]=ent;
607       p->iindices[p->numentries]=indexi;
608       p->jindices[p->numentries]=indexj;
609
610       p->blocksize=blocksize;
611     
612       p_sav->next=p;
613     }
614   } else {
615     if (ent != 0.0){
616         int blksz=C->blocks[blkno].blocksize;
617         if (C->blocks[blkno].blockcategory == DIAG){
618             C->blocks[blkno].data.vec[indexi]=ent;
619         }else{
620           C->blocks[blkno].data.mat[ijtok(indexi,indexj,blksz)]=ent;
621           C->blocks[blkno].data.mat[ijtok(indexj,indexi,blksz)]=ent;
622         };
623     };
624     
625   }
626 }
627