 Algorithmique Numérique Distribuée Public GIT Repository
1 /*      \$Id\$     */
3 /* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved.     */
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. */
9 #include "xbt/sysdep.h"
10 #include "xbt/log.h"
11 #include "maxmin_private.h"
12 #include <stdlib.h>
14 #include <declarations.h>
16 #ifndef MATH
17 #include <math.h>
18 #endif
20 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
21                                 "Logging specific to SURF (sdp)");
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 */
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;
100 #define get_y(a,b) (pow(2,a)+b-1)
102 void sdp_solve(lmm_system_t sys)
103 {
104   lmm_variable_t var = NULL;
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;
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;
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   */
145   if ( !(sys->modified))
146     return;
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));
154   /*
155    * This  number is found based on the tree structure explained on top.
156    */
157   K = (int)log((double)flows)/log(2.0);
159   /*
160    * The number of variables in the SDP style.
161    */
162   nb_var = get_y(K, pow(2,K));
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);
171   /*
172    * The total number of constraints.
173    */
174   nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
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;
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));
188   for(i = 1; i <= nb_var; i++){
189     constraints[i].blocks=NULL;
190   }
192   a = (double *)calloc(nb_var+1, sizeof(double));
194   /*
195    * Block sizes.
196    */
197   block_num=1;
198   block_size=0;
200   /*
201    * For each constraint block do.
202    */
203   for(i = 1; i <= nb_cnsts; i++){
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     }
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));
223     block_num++;
224   }
226   fprintf(sdpout,"\n");
229   /*
230    * Creates de objective function array.
231    */
232   a = (double *)calloc((nb_var+1), sizeof(double));
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");
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);
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);
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);
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);
269       isdiag[block_num] = 0;
270       block_num++;
271     }
272   }
277 }