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;
 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) i++;
 171   }
 172   flows=i;
 173   DEBUG1("Variable set : %d", xbt_swag_size(var_list));
 174   xbt_swag_foreach(var, var_list) {
 175     var->value = 0.0;
 176     if(var->weight) var->index = i--;
 177   }
 178
 179   cnst_list=&(sys->active_constraint_set);
 180   DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
 181
 182   /*
 183    * Those fields are the top level description of the platform furnished in the xml file.
 184    */
 185   links = xbt_swag_size(&(sys->active_constraint_set));
 186
 187   /*
 188    * This  number is found based on the tree structure explained on top.
 189    */
 190   tmp_k = (double) log((double)flows)/log(2.0);
 191   K = (int) ceil(tmp_k);
 192
 193
 194   /*
 195    * The number of variables in the SDP program.
 196    */
 197   nb_var = get_y(K, pow(2,K));
 198   DEBUG1("Number of variables in the SDP program : %d", nb_var);
 199
 200
 201   /*
 202    * Find the size of each group of constraints.
 203    */
 204   nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
 205   nb_cnsts_struct = (int)pow(2,K) - 1;
 206   nb_cnsts_positivy = (int)pow(2,K);
 207
 208   /*
 209    * The total number of constraints.
 210    */
 211   nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
 212   CDEBUG1(surf_sdp_out,"Number of constraints = %d", nb_cnsts);
 213
 214   /*
 215    * Keep track of which blocks have off diagonal entries.
 216    */
 217   isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
 218   for (i=1; i<=nb_cnsts; i++)
 219     isdiag[i]=1;
 220
 221   C.nblocks = nb_cnsts;
 222   C.blocks  = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
 223   constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
 224
 225   for(i = 1; i <= nb_var; i++){
 226     constraints[i].blocks=NULL;
 227   }
 228
 229   a = (double *)calloc(nb_var+1, sizeof(double));
 230
 231   /*
 232    * Block sizes.
 233    */
 234   block_num=1;
 235   block_size=0;
 236
 237   /*
 238    * For each constraint block do.
 239    */
 240   for(i = 1; i <= nb_cnsts; i++){
 241
 242     /*
 243      * Structured blocks are size 2 and all others are size 1.
 244      */
 245     if(i <= nb_cnsts_struct){
 246       total_block_size += block_size = 2;
 247       DEBUG0("2 ");
 248     }else{
 249       total_block_size += block_size = 1;
 250       CDEBUG0(surf_sdp_out,"1 ");
 251     }
 252
 253     /*
 254      * All blocks are matrices.
 255      */
 256     C.blocks[block_num].blockcategory = MATRIX;
 257     C.blocks[block_num].blocksize = block_size;
 258     C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
 259
 260     block_num++;
 261   }
 262
 263   CDEBUG0(surf_sdp_out," ");
 264
 265
 266   /*
 267    * Creates de objective function array.
 268    */
 269   a = (double *)calloc((nb_var+1), sizeof(double));
 270
 271   for(i = 1; i <= nb_var; i++){
 272     if(get_y(0,1)==i){
 273       CDEBUG0(surf_sdp_out,"-1 ");
 274       a[i]=-1;
 275     }else{
 276       CDEBUG0(surf_sdp_out,"0 ");
 277       a[i]=0;
 278     }
 279   }
 280   CDEBUG0(surf_sdp_out,"\n");
 281
 282
 283   /*
 284    * Structure contraint blocks.
 285    */
 286   block_num = 1;
 287   matno = 1;
 288   for(k = 1; k <= K; k++){
 289     for(i = 1; i <= pow(2,k-1); i++){
 290       matno=get_y(k,2*i-1);
 291       CDEBUG2(surf_sdp_out,"%d %d 1 1 1\n", matno  , block_num);
 292       addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
 293
 294       matno=get_y(k,2*i);
 295       CDEBUG2(surf_sdp_out,"%d %d 2 2 1\n", matno  , block_num);
 296       addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
 297
 298       matno=get_y(k-1,i);
 299       CDEBUG2(surf_sdp_out,"%d %d 1 2 1\n", matno  , block_num);
 300       addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);
 301
 302       matno=get_y(k-1,i);
 303       CDEBUG2(surf_sdp_out,"%d %d 2 1 1\n", matno  , block_num);
 304       addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
 305
 306       isdiag[block_num] = 0;
 307       block_num++;
 308     }
 309   }
 310
 311
 312   /*
 313    * Capacity constraint block.
 314    */
 315   xbt_swag_foreach(cnst, cnst_list) {
 316
 317     CDEBUG2(surf_sdp_out,"0 %d 1 1 %f\n", block_num,  - (cnst->bound));
 318     addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);
 319
 320     elem_list = &(cnst->element_set);
 321     xbt_swag_foreach(elem, elem_list){
 322       if(elem->variable->weight <=0) break;
 323       matno=get_y(K,elem->variable->index);
 324       CDEBUG3(surf_sdp_out,"%d %d 1 1 %f\n", elem->variable->index, block_num, - (elem->value));
 325       addentry(constraints, &C, matno, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
 326
 327     }
 328     block_num++;
 329   }
 330
 331
 332   /*
 333    * Positivy constraint blocks.
 334    */
 335   for(i = 1; i <= pow(2,K); i++){
 336     matno=get_y(K, i);
 337     CDEBUG2(surf_sdp_out,"%d %d 1 1 1\n", matno, block_num);
 338     addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
 339     block_num++;
 340   }
 341
 342
 343   /*
 344    * At this point, we'll stop to recognize whether any of the blocks
 345    * are "hidden LP blocks"  and correct the block type if needed.
 346    */
 347   for (i=1; i<=nb_cnsts; i++){
 348     if ((C.blocks[i].blockcategory != DIAG) &&
 349         (isdiag[i]==1) && (C.blocks[i].blocksize > 1)){
 350       /*
 351        * We have a hidden diagonal block!
 352        */
 353
 354       //printf("Block %d is actually diagonal.\n",i);
 355
 356       blocksz=C.blocks[i].blocksize;
 357       tempdiag=(double *)calloc((blocksz+1), sizeof(double));
 358       for (j=1; j<=blocksz; j++)
 359         tempdiag[j]=C.blocks[i].data.mat[ijtok(j,j,blocksz)];
 360       free(C.blocks[i].data.mat);
 361       C.blocks[i].data.vec=tempdiag;
 362       C.blocks[i].blockcategory=DIAG;
 363     };
 364   };
 365
 366
 367   /*
 368    * Next, setup issparse and NULL out all nextbyblock pointers.
 369    */
 370   p=NULL;
 371   for (i=1; i<=k; i++) {
 372     p=constraints[i].blocks;
 373     while (p != NULL){
 374           /*
 375            * First, set issparse.
 376            */
 377           if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15)){
 378             p->issparse=0;
 379           }else{
 380             p->issparse=1;
 381           };
 382
 383           if (C.blocks[p->blocknum].blockcategory == DIAG)
 384             p->issparse=1;
 385
 386           /*
 387            * Setup the cross links.
 388            */
 389
 390           p->nextbyblock=NULL;
 391           p=p->next;
 392     };
 393   };
 394
 395
 396   /*
 397    * Create cross link reference.
 398    */
 399   create_cross_link(constraints, nb_var);
 400
 401
 402   /*
 403    * Debuging print problem in SDPA format.
 404    */
 405   if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
 406     DEBUG0("Printing SDPA...");
 407     tmp=strdup("SURF-PROPORTIONNAL.sdpa");
 408     write_prob(tmp,total_block_size,nb_var,C,a,constraints);
 409     free(tmp);
 410   }
 411
 412   /*
 413    * Initialize parameters.
 414    */
 415   DEBUG0("Initializing solution...");
 416   initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
 417
 418
 419
 420   /*
 421    * Call the solver.
 422    */
 423   DEBUG0("Calling the solver...");
 424   stdout_sav=stdout;
 425   stdout=fopen("/dev/null","w");
 426   ret = easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
 427   fclose(stdout);
 428   stdout=stdout_sav;
 429
 430   switch(ret){
 431   case 0:
 432   case 1: DEBUG0("SUCCESS The problem is primal infeasible");
 433           break;
 434
 435   case 2: DEBUG0("SUCCESS The problem is dual infeasible");
 436           break;
 437
 438   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.");
 439           break;
 440
 441   case 4: DEBUG0("Failure. Maximum number of iterations reached.");
 442           break;
 443
 444   case 5: DEBUG0("Failure. Stuck at edge of primal feasibility.");
 445           break;
 446
 447   }
 448
 449   if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
 450     tmp=strdup("SURF-PROPORTIONNAL.sol");
 451     write_sol(tmp,total_block_size, nb_var, X, y, Z);
 452     free(tmp);
 453   }
 454
 455   /*
 456    * Write out the solution if necessary.
 457    */
 458   xbt_swag_foreach(cnst, cnst_list) {
 459
 460     elem_list = &(cnst->element_set);
 461     xbt_swag_foreach(elem, elem_list) {
 462       if(elem->variable->weight <=0) break;
 463
 464       i = (int)get_y(K, elem->variable->index);
 465       elem->variable->value = y[i];
 466
 467     }
 468   }
 469
 470
 471   /*
 472    * Free up memory.
 473    */
 474   free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
 475
 476
 477  free(isdiag);
 478  free(tempdiag);
 479  free(tmp);
 480  free(a);
 481  free(y);
 482
 483  //  fclose(sdpout);
 484  free(isdiag);
 485  sys->modified = 0;
 486
 487  if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
 488    lmm_print(sys);
 489  }
 490
 491 }
 492
 493
 494 /*
 495  * Create the cross_link reference in order to have a byblock list.
 496  */
 497 void create_cross_link(struct constraintmatrix *myconstraints, int k){
 498
 499   int i, j;
 500   int blk;
 501   struct sparseblock *p;
 502   struct sparseblock *q;
 503
 504   struct sparseblock *prev;
 505
 506   /*
 507    * Now, cross link.
 508    */
 509   prev=NULL;
 510   for (i=1; i<=k; i++){
 511     p=myconstraints[i].blocks;
 512     while (p != NULL){
 513       if (p->nextbyblock == NULL){
 514         blk=p->blocknum;
 515
 516         /*
 517          * link in the remaining blocks.
 518          */
 519         for (j=i+1; j<=k; j++){
 520           q=myconstraints[j].blocks;
 521
 522           while (q != NULL){
 523             if (q->blocknum == p->blocknum){
 524               if (p->nextbyblock == NULL){
 525                 p->nextbyblock=q;
 526                 q->nextbyblock=NULL;
 527                 prev=q;
 528               }
 529               else{
 530                 prev->nextbyblock=q;
 531                 q->nextbyblock=NULL;
 532                 prev=q;
 533               }
 534               break;
 535             }
 536             q=q->next;
 537           }
 538         }
 539       }
 540       p=p->next;
 541     }
 542   }
 543 }
 544
 545
 546
 547
 548 void addentry(struct constraintmatrix *constraints,
 549               struct blockmatrix *C,
 550               int matno,
 551               int blkno,
 552               int indexi,
 553               int indexj,
 554               double ent,
 555               int blocksize)
 556 {
 557   struct sparseblock *p;
 558   struct sparseblock *p_sav;
 559
 560   p=constraints[matno].blocks;
 561
 562   if (matno != 0.0) {
 563     if (p == NULL){
 564       /*
 565        * We haven't yet allocated any blocks.
 566        */
 567       p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
 568
 569       //two entries because this library ignores indices starting in zerox
 570       p->constraintnum=matno;
 571       p->blocknum=blkno;
 572       p->numentries=1;
 573       p->next=NULL;
 574
 575       p->entries=calloc(p->numentries+1, sizeof(double));
 576       p->iindices=calloc(p->numentries+1, sizeof(int));
 577       p->jindices=calloc(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       p->blocksize=blocksize;
 584
 585       constraints[matno].blocks=p;
 586     } else {
 587       /*
 588        * We have some existing blocks.  See whether this block is already
 589        * in the chain.
 590        */
 591       while (p != NULL){
 592         if (p->blocknum == blkno){
 593           /*
 594            * Found the right block.
 595            */
 596           p->constraintnum=matno;
 597           p->blocknum=blkno;
 598           p->numentries=p->numentries+1;
 599
 600           p->entries = realloc(p->entries,  (p->numentries+1) * sizeof(double) );
 601           p->iindices = realloc(p->iindices, (p->numentries+1) * sizeof(int) );
 602           p->jindices = realloc(p->jindices, (p->numentries+1) * sizeof(int) );
 603
 604           p->entries[p->numentries]=ent;
 605           p->iindices[p->numentries]=indexi;
 606           p->jindices[p->numentries]=indexj;
 607
 608           return;
 609         }
 610         p_sav=p;
 611         p=p->next;
 612       }
 613
 614       /*
 615        * If we get here, we have a non-empty structure but not the right block
 616        * inside hence create a new structure.
 617        */
 618
 619       p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
 620
 621       //two entries because this library ignores indices starting in zerox
 622       p->constraintnum=matno;
 623       p->blocknum=blkno;
 624       p->numentries=1;
 625       p->next=NULL;
 626
 627       p->entries=calloc(p->numentries+1, sizeof(double));
 628       p->iindices=calloc(p->numentries+1, sizeof(int));
 629       p->jindices=calloc(p->numentries+1, sizeof(int));
 630
 631       p->entries[p->numentries]=ent;
 632       p->iindices[p->numentries]=indexi;
 633       p->jindices[p->numentries]=indexj;
 634
 635       p->blocksize=blocksize;
 636
 637       p_sav->next=p;
 638     }
 639   } else {
 640     if (ent != 0.0){
 641         int blksz=C->blocks[blkno].blocksize;
 642         if (C->blocks[blkno].blockcategory == DIAG){
 643             C->blocks[blkno].data.vec[indexi]=ent;
 644         }else{
 645           C->blocks[blkno].data.mat[ijtok(indexi,indexj,blksz)]=ent;
 646           C->blocks[blkno].data.mat[ijtok(indexj,indexi,blksz)]=ent;
 647         };
 648     };
 649
 650   }
 651 }
 652