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