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