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