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