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