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