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/sysdep.h"
  10 #include "xbt/log.h"
  11 #include "maxmin_private.h"
  12 #include <stdlib.h>
  13
  14 #include <declarations.h>
  15
  16 #ifndef MATH
  17 #include <math.h>
  18 #endif
  19
  20 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
  21                                 "Logging specific to SURF (sdp)");
  22
  23
  24
  25
  26
  27 /*
  28 ########################################################################
  29 ######################## Simple Proportionnal fairness #################
  30 ########################################################################
  31 ####### Original problem ##########
  32 #
  33 #  Max x_1*x_2*...*x_n
  34 #    (A.x)_1 <= b_1
  35 #    (A.x)_2 <= b_2
  36 #     ...
  37 #    (A.x)_m <= b_m
  38 #    x>=0
  39 #
  40 #  We assume in the following that n=2^K
  41 #
  42 ####### Standard SDP form #########
  43 #
  44 # A SDP can be written in the following standard form:
  45 #
  46 #   (P) min c1*x1+c2*x2+...+cm*xm
  47 #       st  F1*x1+F2*x2+...+Fm*xm-F0=X
  48 #                                 X >= 0
  49 #
  50 # Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
  51 # constraint X>0 means that X must be positive semidefinite positive.
  52 #
  53 ########## Equivalent SDP #########
  54 #
  55 #  Variables:
  56 #
  57 #    y(k,i) for k in [0,K] and i in [1,2^k]
  58 #
  59 #  Structure :
  60 #                             y(0,1)
  61 #                y(1,1)                      y(1,2)
  62 #        y(2,1)        y(2,2)        y(2,3)        y(2,4)
  63 #    y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
  64 #    -------------------------------------------------------
  65 #     x_1     x_2    x_3    x_4    x_5    x_6    x_7    x_8
  66 #
  67 #
  68 #  Structure Constraints:
  69 #
  70 #    forall k in [0,K-1], and i in [1,2^k]
  71 #      [ y(k,  2i-1)   y(k-1, i) ]
  72 #      [ y(k-1, i  )   y(k,  2i) ]
  73 #
  74 #    2^K-1
  75 #
  76 #  Capacity Constraints:
  77 #     forall l in [1,m]
  78 #        -(A.y(K,*))_l >= - b_l
  79 #
  80 #  Positivity Constraints:
  81 #    forall k in [0,K], and i in [1,2^k]
  82 #        y(k,i) >= 0
  83 #
  84 #  Minimize -y(0,1)
  85 */
  86
  87
  88
  89
  90 //typedef struct lmm_system {
  91 //  int modified;
  92 //  s_xbt_swag_t variable_set;  /* a list of lmm_variable_t */
  93 //  s_xbt_swag_t constraint_set;        /* a list of lmm_constraint_t */
  94 //  s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */
  95 //  s_xbt_swag_t saturated_variable_set;        /* a list of lmm_variable_t */
  96 //  s_xbt_swag_t saturated_constraint_set;      /* a list of lmm_constraint_t_t */
  97 //  xbt_mallocator_t variable_mallocator;
  98 //} s_lmm_system_t;
  99
 100 #define get_y(a,b) (pow(2,a)+b-1)
 101
 102 void sdp_solve(lmm_system_t sys)
 103 {
 104   lmm_variable_t var = NULL;
 105
 106   /*
 107    * SDP mapping variables.
 108    */
 109   int K=0;
 110   int nb_var = 0;
 111   int nb_cnsts = 0;
 112   int flows = 0;
 113   int links = 0;
 114   int nb_cnsts_capacity=0;
 115   int nb_cnsts_struct=0;
 116   int nb_cnsts_positivy=0;
 117   int block_num=0;
 118   int block_size;
 119   int *isdiag=NULL;
 120   FILE *sdpout = fopen("SDPA-printf.tmp","w");
 121   int matno=0;
 122   int i=0;
 123   int j=0;
 124   int k=0;
 125
 126   /*
 127    * CSDP library specific variables.
 128    */
 129   struct blockmatrix C;
 130   struct blockmatrix X,Z;
 131   double *y;
 132   double pobj, dobj;
 133   double *a;
 134   struct constraintmatrix *constraints;
 135
 136   /*
 137   lmm_constraint_t cnst = NULL;
 138   lmm_element_t elem = NULL;
 139   xbt_swag_t cnst_list = NULL;
 140   xbt_swag_t var_list = NULL;
 141   xbt_swag_t elem_list = NULL;
 142   double min_usage = -1;
 143   */
 144
 145   if ( !(sys->modified))
 146     return;
 147
 148   /*
 149    * Those fields are the top level description of the platform furnished in the xml file.
 150    */
 151   flows = xbt_swag_size(&(sys->variable_set));
 152   links = xbt_swag_size(&(sys->active_constraint_set));
 153
 154   /*
 155    * This  number is found based on the tree structure explained on top.
 156    */
 157   K = (int)log((double)flows)/log(2.0);
 158
 159   /*
 160    * The number of variables in the SDP style.
 161    */
 162   nb_var = get_y(K, pow(2,K));
 163
 164   /*
 165    * Find the size of each group of constraints.
 166    */
 167   nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
 168   nb_cnsts_struct = (int)pow(2,K) - 1;
 169   nb_cnsts_positivy = (int)pow(2,K);
 170
 171   /*
 172    * The total number of constraints.
 173    */
 174   nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
 175
 176
 177   /*
 178    * Keep track of which blocks have off diagonal entries.
 179    */
 180   isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
 181   for (i=1; i<=nb_cnsts; i++)
 182     isdiag[i]=1;
 183
 184   C.nblocks = nb_cnsts;
 185   C.blocks  = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
 186   constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
 187
 188   for(i = 1; i <= nb_var; i++){
 189     constraints[i].blocks=NULL;
 190   }
 191
 192   a = (double *)calloc(nb_var+1, sizeof(double));
 193
 194   /*
 195    * Block sizes.
 196    */
 197   block_num=1;
 198   block_size=0;
 199
 200   /*
 201    * For each constraint block do.
 202    */
 203   for(i = 1; i <= nb_cnsts; i++){
 204
 205     /*
 206      * Structured blocks are size 2 and all others are size 1.
 207      */
 208     if(i <= nb_cnsts_struct){
 209       block_size = 2;
 210       fprintf(sdpout,"2 ");
 211     }else{
 212       block_size = 1;
 213       fprintf(sdpout,"1 ");
 214     }
 215
 216     /*
 217      * All blocks are matrices.
 218      */
 219     C.blocks[block_num].blockcategory = MATRIX;
 220     C.blocks[block_num].blocksize = block_size;
 221     C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
 222
 223     block_num++;
 224   }
 225
 226   fprintf(sdpout,"\n");
 227
 228
 229   /*
 230    * Creates de objective function array.
 231    */
 232   a = (double *)calloc((nb_var+1), sizeof(double));
 233
 234   for(i = 1; i <= nb_var; i++){
 235     if(get_y(0,1)==i){
 236       fprintf(sdpout,"-1 ");
 237       a[i]=-1;
 238     }else{
 239       fprintf(sdpout,"0 ");
 240       a[i]=0;
 241     }
 242   }
 243   fprintf(sdpout,"\n");
 244
 245
 246   /*
 247    * Structure contraint blocks.
 248    */
 249   block_num = 1;
 250   matno = 1;
 251   for(k = 1; k <= K; k++){
 252     for(i = 1; i <= pow(2,k-1); i++){
 253       matno=get_y(k,2*i-1);
 254       fprintf(sdpout,"%d %d 1 1 1\n", matno  , block_num);
 255       addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
 256
 257       matno=get_y(k,2*i);
 258       fprintf(sdpout,"%d %d 2 2 1\n", matno  , block_num);
 259       addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
 260
 261       matno=get_y(k-1,i);
 262       fprintf(sdpout,"%d %d 1 2 1\n", matno  , block_num);
 263       addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);
 264
 265       matno=get_y(k-1,i);
 266       fprintf(sdpout,"%d %d 2 1 1\n", matno  , block_num);
 267       addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
 268
 269       isdiag[block_num] = 0;
 270       block_num++;
 271     }
 272   }
 273
 274
 275
 276
 277 }