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