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