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