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