Logo AND Algorithmique Numérique Distribuée

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