Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Shut libcsdp mouth!
[simgrid.git] / src / surf / sdp.c
index 6d54e2e..35a49bf 100644 (file)
  * under the terms of the license (GNU LGPL) which comes with this package. */
 
 
-#include "xbt/sysdep.h"
 #include "xbt/log.h"
+#include "xbt/sysdep.h"
+#include "xbt/mallocator.h"
 #include "maxmin_private.h"
+
 #include <stdlib.h>
+#ifndef MATH
+#include <math.h>
+#endif
+
+/*
+ * SDP specific variables.
+ */
+#include <declarations.h>
+
+
+static void create_cross_link(struct constraintmatrix *myconstraints, int k);
+
+static void addentry(struct constraintmatrix *constraints, 
+             struct blockmatrix *, int matno, int blkno, int indexi, int indexj, double ent, int blocksize);
+
+
 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
                                "Logging specific to SURF (sdp)");
+XBT_LOG_NEW_SUBCATEGORY(surf_sdp_out, surf,
+                               "Logging specific to SURF (sdp)");
+/*
+########################################################################
+######################## Simple Proportionnal fairness #################
+########################################################################
+####### Original problem ##########
+#
+#  Max x_1*x_2*...*x_n
+#    (A.x)_1 <= b_1
+#    (A.x)_2 <= b_2
+#     ...
+#    (A.x)_m <= b_m
+#    x>=0
+#
+#  We assume in the following that n=2^K
+#
+####### Standard SDP form #########
+#
+# A SDP can be written in the following standard form:
+# 
+#   (P) min c1*x1+c2*x2+...+cm*xm
+#       st  F1*x1+F2*x2+...+Fm*xm-F0=X
+#                                 X >= 0
+#
+# Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
+# constraint X>0 means that X must be positive semidefinite positive.
+#
+########## Equivalent SDP #########
+#
+#  Variables:
+#
+#    y(k,i) for k in [0,K] and i in [1,2^k]
+#
+#  Structure :
+#                             y(0,1)
+#                y(1,1)                      y(1,2)
+#        y(2,1)        y(2,2)        y(2,3)        y(2,4)
+#    y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
+#    -------------------------------------------------------
+#     x_1     x_2    x_3    x_4    x_5    x_6    x_7    x_8
+#  
+#
+#  Structure Constraints:
+#
+#    forall k in [0,K-1], and i in [1,2^k]
+#      [ y(k,  2i-1)   y(k-1, i) ]
+#      [ y(k-1, i  )   y(k,  2i) ]
+#
+#    2^K-1
+#
+#  Capacity Constraints:
+#     forall l in [1,m]
+#        -(A.y(K,*))_l >= - b_l
+#
+#  Positivity Constraints:
+#    forall k in [0,K], and i in [1,2^k]
+#        y(k,i) >= 0
+#
+#  Minimize -y(0,1)
+*/
+
+//typedef struct lmm_system {
+//  int modified;
+//  s_xbt_swag_t variable_set; /* a list of lmm_variable_t */
+//  s_xbt_swag_t constraint_set;       /* a list of lmm_constraint_t */
+//  s_xbt_swag_t active_constraint_set;        /* a list of lmm_constraint_t */
+//  s_xbt_swag_t saturated_variable_set;       /* a list of lmm_variable_t */
+//  s_xbt_swag_t saturated_constraint_set;     /* a list of lmm_constraint_t_t */
+//  xbt_mallocator_t variable_mallocator;
+//} s_lmm_system_t;
+
+#define get_y(a,b) (pow(2,a)+b-1)
 
 void sdp_solve(lmm_system_t sys)
 {
-/*   lmm_variable_t var = NULL; */
-/*   lmm_constraint_t cnst = NULL; */
-/*   lmm_element_t elem = NULL; */
-/*   xbt_swag_t cnst_list = NULL; */
-/*   xbt_swag_t var_list = NULL; */
-/*   xbt_swag_t elem_list = NULL; */
-/*   double min_usage = -1; */
-
-  if (!(sys->modified))
+  lmm_variable_t var = NULL;
+
+  /*
+   * SDP mapping variables.
+   */
+  int K=0;
+  int nb_var = 0;
+  int nb_cnsts = 0;
+  int flows = 0;
+  int links = 0;
+  int nb_cnsts_capacity=0;
+  int nb_cnsts_struct=0;
+  int nb_cnsts_positivy=0;
+  int block_num=0;
+  int block_size; 
+  int total_block_size=0;
+  int *isdiag=NULL;
+  //  FILE *sdpout = fopen("SDPA-printf.tmp","w");
+  int blocksz = 0;
+  double *tempdiag = NULL;
+  int matno=0;
+  int i=0;
+  int j=0;
+  int k=0;
+  
+  /* 
+   * CSDP library specific variables.
+   */
+  struct blockmatrix C;
+  struct blockmatrix X,Z;
+  double *y;
+  double pobj, dobj;
+  double *a;
+  struct constraintmatrix *constraints;  
+  /*
+   * Classic maxmin variables.
+   */
+  lmm_constraint_t cnst = NULL;
+  lmm_element_t elem = NULL;
+  xbt_swag_t cnst_list = NULL;
+  xbt_swag_t var_list = NULL;
+  xbt_swag_t elem_list = NULL;
+
+  if ( !(sys->modified))
     return;
 
+  /* 
+   * Initialize the var list variable with only the active variables. 
+   * Associate an index in the swag variables.
+   */
+  var_list = &(sys->variable_set);
+  i=0;
+  xbt_swag_foreach(var, var_list) {
+    if(var->weight) i++;
+  }
+  flows=i;
+  DEBUG1("Variable set : %d", xbt_swag_size(var_list));
+  xbt_swag_foreach(var, var_list) {
+    var->value = 0.0;
+    if(var->weight) var->index = i--;
+  }
+
+  cnst_list=&(sys->active_constraint_set); 
+  DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
+
+  /*
+   * Those fields are the top level description of the platform furnished in the xml file.
+   */
+  links = xbt_swag_size(&(sys->active_constraint_set));
+
+  /* 
+   * This  number is found based on the tree structure explained on top.
+   */
+  double tmp_k;
+
+  tmp_k = (double) log((double)flows)/log(2.0);
+  K = (int) ceil(tmp_k);
+
+
+  /* 
+   * The number of variables in the SDP program. 
+   */
+  nb_var = get_y(K, pow(2,K));
+  DEBUG1("Number of variables in the SDP program : %d", nb_var);
+  CDEBUG1(surf_sdp_out,"%d", nb_var);
+  /*
+   * Find the size of each group of constraints.
+   */
+  nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
+  nb_cnsts_struct = (int)pow(2,K) - 1;
+  nb_cnsts_positivy = (int)pow(2,K);
+
+  /*
+   * The total number of constraints.
+   */
+  nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
+  //fprintf(sdpout,"%d\n", nb_cnsts);
+
+  /*
+   * Keep track of which blocks have off diagonal entries. 
+   */
+  isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
+  for (i=1; i<=nb_cnsts; i++)
+    isdiag[i]=1;
+
+  C.nblocks = nb_cnsts;
+  C.blocks  = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
+  constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
+
+  for(i = 1; i <= nb_var; i++){
+    constraints[i].blocks=NULL; 
+  }
+  
+  a = (double *)calloc(nb_var+1, sizeof(double)); 
+
+  /*
+   * Block sizes.
+   */
+  block_num=1;
+  block_size=0;
+
+  /*
+   * For each constraint block do.
+   */
+  for(i = 1; i <= nb_cnsts; i++){
+    
+    /*
+     * Structured blocks are size 2 and all others are size 1.
+     */
+    if(i <= nb_cnsts_struct){
+      total_block_size += block_size = 2;
+      //fprintf(sdpout,"2 ");
+    }else{
+      total_block_size += block_size = 1;
+      //fprintf(sdpout,"1 ");
+    }
+    
+    /*
+     * All blocks are matrices.
+     */
+    C.blocks[block_num].blockcategory = MATRIX;
+    C.blocks[block_num].blocksize = block_size;
+    C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
+    block_num++;
+  }
+
+  //fprintf(sdpout,"\n");
+
+
+  /*
+   * Creates de objective function array.
+   */
+  a = (double *)calloc((nb_var+1), sizeof(double));
+  
+  for(i = 1; i <= nb_var; i++){
+    if(get_y(0,1)==i){
+      //fprintf(sdpout,"-1 ");
+      a[i]=-1;
+    }else{
+      //fprintf(sdpout,"0 ");
+      a[i]=0;
+    }
+  }
+  //fprintf(sdpout,"\n");
+
+
+  /*
+   * Structure contraint blocks.
+   */
+  block_num = 1;
+  matno = 1;
+  for(k = 1; k <= K; k++){
+    for(i = 1; i <= pow(2,k-1); i++){
+      matno=get_y(k,2*i-1);
+      //fprintf(sdpout,"%d %d 1 1 1\n", matno  , block_num);
+      addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
+
+      matno=get_y(k,2*i);
+      //fprintf(sdpout,"%d %d 2 2 1\n", matno  , block_num);
+      addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
+
+      matno=get_y(k-1,i);
+      //fprintf(sdpout,"%d %d 1 2 1\n", matno  , block_num);
+      addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);      
+      
+      matno=get_y(k-1,i);
+      //fprintf(sdpout,"%d %d 2 1 1\n", matno  , block_num);
+      addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
+      
+      isdiag[block_num] = 0;
+      block_num++;
+    }
+  }
+
+  /*
+   * Capacity constraint block.
+   */
+  xbt_swag_foreach(cnst, cnst_list) {
+
+    //fprintf(sdpout,"0 %d 1 1 %f\n", block_num,  - (cnst->bound));    
+    addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);    
+
+
+    elem_list = &(cnst->element_set);
+    xbt_swag_foreach(elem, elem_list) {
+      if(elem->variable->weight <=0) break;
+      matno=get_y(K,elem->variable->index);
+      //fprintf(sdpout,"%d %d 1 1 %f\n", elem->variable->index, block_num, - (elem->value)); 
+      addentry(constraints, &C, matno, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
+    }
+    block_num++;
+  }
+
+  
+  /*
+   * Positivy constraint blocks.
+   */
+  for(i = 1; i <= pow(2,K); i++){
+    matno=get_y(K, i);
+    //fprintf(sdpout,"%d %d 1 1 1\n", matno, block_num);
+    addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
+    block_num++;
+  }
+
+
+  /*
+   * At this point, we'll stop to recognize whether any of the blocks
+   * are "hidden LP blocks"  and correct the block type if needed.
+   */
+  for (i=1; i<=nb_cnsts; i++){
+    if ((C.blocks[i].blockcategory != DIAG) && 
+       (isdiag[i]==1) && (C.blocks[i].blocksize > 1)){
+      /*
+       * We have a hidden diagonal block!
+       */
+      
+      //printf("Block %d is actually diagonal.\n",i);
+      
+      blocksz=C.blocks[i].blocksize;
+      tempdiag=(double *)calloc((blocksz+1), sizeof(double));
+      for (j=1; j<=blocksz; j++)
+       tempdiag[j]=C.blocks[i].data.mat[ijtok(j,j,blocksz)];
+      free(C.blocks[i].data.mat);
+      C.blocks[i].data.vec=tempdiag;
+      C.blocks[i].blockcategory=DIAG;
+    };
+  };
+  
+  
+  /*
+   * Next, setup issparse and NULL out all nextbyblock pointers.
+   */
+  struct sparseblock *p=NULL;
+  for (i=1; i<=k; i++) {
+    p=constraints[i].blocks;
+    while (p != NULL){
+         /*
+          * First, set issparse.
+          */
+         if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15)){
+           p->issparse=0;
+         }else{
+           p->issparse=1;
+         };
+         
+         if (C.blocks[p->blocknum].blockcategory == DIAG)
+           p->issparse=1;
+         
+         /*
+          * Setup the cross links.
+          */
+         
+         p->nextbyblock=NULL;
+         p=p->next;
+    };
+  };
+
+
+  /*
+   * Create cross link reference.
+   */
+  create_cross_link(constraints, nb_var);
+  
+  /*
+   * Debuging print problem in SDPA format.
+   */
+  if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
+    DEBUG0("Printing SDPA...\n");
+    char *tmp=strdup("SURF-PROPORTIONNAL.sdpa");
+    write_prob(tmp,total_block_size,nb_var,C,a,constraints);
+    free(tmp);
+  }
+
+  /*
+   * Initialize parameters.
+   */
+  DEBUG0("Initializing solution...\n");
+  initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);  
+  
+
+
+  /*
+   * Call the solver.
+   */
+  DEBUG0("Calling the solver...\n");
+  FILE *stdout_sav=stdout;
+  stdout=fopen("/dev/null","w");
+  int ret = easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
+  fclose(stdout);
+  stdout=stdout_sav;
+
+  switch(ret){
+  case 0:
+  case 1: DEBUG0("SUCCESS The problem is primal infeasible\n");
+          break;
+
+  case 2: DEBUG0("SUCCESS The problem is dual infeasible\n");
+          break;
+  case 3: DEBUG0("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");
+          break;
+
+  case 4: DEBUG0("Failure. Maximum number of iterations reached.");
+          break;
+
+  case 5: DEBUG0("Failure. Stuck at edge of primal feasibility.");
+          break;
+
+  }
+
+  if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
+    char *tmp=strdup("SURF-PROPORTIONNAL.sol");
+    write_sol(tmp,total_block_size, nb_var, X, y, Z);
+    free(tmp);
+  }
+
+  /*
+   * Write out the solution if necessary.
+   */
+  xbt_swag_foreach(cnst, cnst_list) {
+
+    elem_list = &(cnst->element_set);
+    xbt_swag_foreach(elem, elem_list) {
+      if(elem->variable->weight <=0) break;
+      
+      i = (int)get_y(K, elem->variable->index);
+      elem->variable->value = y[i]; 
+      
+    }
+  }
+
+
+  /*
+   * Free up memory.
+   */
+  free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
+
+  //  fclose(sdpout);
+  free(isdiag);
+  sys->modified = 0;
+
+  if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
+    lmm_print(sys);
+  }
+
 }
+
+
+/*
+ * Create the cross_link reference in order to have a byblock list.
+ */
+void create_cross_link(struct constraintmatrix *myconstraints, int k){
+
+  int i, j;
+  int blk;
+  struct sparseblock *p;
+  struct sparseblock *q;
+
+  struct sparseblock *prev;
+
+  /*
+   * Now, cross link.
+   */
+  prev=NULL;
+  for (i=1; i<=k; i++){
+    p=myconstraints[i].blocks;
+    while (p != NULL){
+      if (p->nextbyblock == NULL){
+       blk=p->blocknum;
+       
+       /*
+        * link in the remaining blocks.
+        */
+       for (j=i+1; j<=k; j++){
+         q=myconstraints[j].blocks;
+                 
+         while (q != NULL){
+           if (q->blocknum == p->blocknum){
+             if (p->nextbyblock == NULL){
+               p->nextbyblock=q;
+               q->nextbyblock=NULL;
+               prev=q;
+             }
+             else{
+               prev->nextbyblock=q;
+               q->nextbyblock=NULL;
+               prev=q;
+             }
+             break;
+           }
+           q=q->next;
+         }
+       }
+      }
+      p=p->next;
+    }
+  }
+}
+
+
+
+void addentry(struct constraintmatrix *constraints, 
+             struct blockmatrix *C, 
+             int matno,
+             int blkno,
+             int indexi,
+             int indexj,
+             double ent, 
+             int blocksize)
+{
+  struct sparseblock *p;
+  struct sparseblock *p_sav;
+
+  p=constraints[matno].blocks;
+  
+  if (matno != 0.0) {
+    if (p == NULL){
+      /*
+       * We haven't yet allocated any blocks.
+       */
+      p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
+      
+      //two entries because this library ignores indices starting in zerox
+      p->constraintnum=matno;
+      p->blocknum=blkno;
+      p->numentries=1;
+      p->next=NULL;
+
+      p->entries=calloc(p->numentries+1, sizeof(double));
+      p->iindices=calloc(p->numentries+1, sizeof(int));
+      p->jindices=calloc(p->numentries+1, sizeof(int));
+   
+      p->entries[p->numentries]=ent;
+      p->iindices[p->numentries]=indexi;
+      p->jindices[p->numentries]=indexj;
+
+      p->blocksize=blocksize;
+
+      constraints[matno].blocks=p;
+    } else {
+      /*
+       * We have some existing blocks.  See whether this block is already
+       * in the chain.
+       */
+      while (p != NULL){
+       if (p->blocknum == blkno){
+         /*
+          * Found the right block.
+          */
+         p->constraintnum=matno;
+         p->blocknum=blkno;
+         p->numentries=p->numentries+1;
+
+         p->entries = realloc(p->entries,  (p->numentries+1) * sizeof(double) );
+         p->iindices = realloc(p->iindices, (p->numentries+1) * sizeof(int) );
+         p->jindices = realloc(p->jindices, (p->numentries+1) * sizeof(int) );
+   
+         p->entries[p->numentries]=ent;
+         p->iindices[p->numentries]=indexi;
+         p->jindices[p->numentries]=indexj;
+       
+         return;
+       }
+       p_sav=p;
+       p=p->next;
+      }
+
+      /*
+       * If we get here, we have a non-empty structure but not the right block
+       * inside hence create a new structure.
+       */
+      
+      p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
+        
+      //two entries because this library ignores indices starting in zerox
+      p->constraintnum=matno;
+      p->blocknum=blkno;
+      p->numentries=1;
+      p->next=NULL;
+
+      p->entries=calloc(p->numentries+1, sizeof(double));
+      p->iindices=calloc(p->numentries+1, sizeof(int));
+      p->jindices=calloc(p->numentries+1, sizeof(int));
+   
+      p->entries[p->numentries]=ent;
+      p->iindices[p->numentries]=indexi;
+      p->jindices[p->numentries]=indexj;
+
+      p->blocksize=blocksize;
+    
+      p_sav->next=p;
+    }
+  } else {
+    if (ent != 0.0){
+       int blksz=C->blocks[blkno].blocksize;
+       if (C->blocks[blkno].blockcategory == DIAG){
+           C->blocks[blkno].data.vec[indexi]=ent;
+       }else{
+         C->blocks[blkno].data.mat[ijtok(indexi,indexj,blksz)]=ent;
+         C->blocks[blkno].data.mat[ijtok(indexj,indexi,blksz)]=ent;
+       };
+    };
+    
+  }
+}
+