Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
01e22449c2070d0e538d4d89c033f9a1e4c29635
[simgrid.git] / src / surf / lagrange.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  * Modelling the proportional fairness using the Lagrange Optimization 
10  * Approach. For a detailed description see:
11  * "ssh://username@scm.gforge.inria.fr/svn/memo/people/pvelho/lagrange/ppf.ps".
12  */
13 #include "xbt/log.h"
14 #include "xbt/sysdep.h"
15 #include "xbt/mallocator.h"
16 #include "maxmin_private.h"
17
18 #include <stdlib.h>
19 #ifndef MATH
20 #include <math.h>
21 #endif
22
23 #define LAMBDA_STEP 0.01
24
25
26 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
27                                 "Logging specific to SURF (lagrange)");
28
29 XBT_LOG_NEW_SUBCATEGORY(surf_writelambda, surf,
30                                 "Generates the lambda.in file. WARNING: the size of this file might be a few GBs.");
31
32 void lagrange_solve(lmm_system_t sys);
33
34 void lagrange_solve(lmm_system_t sys)
35 {
36   /*
37    * Lagrange Variables.
38    */
39   int max_iterations= 1000000;
40   double epsilon_min_error = 0.00001;
41   double overall_error = 1;
42   double sigma_step = LAMBDA_STEP;
43   //double capacity_error=0, bound_error=0;
44   int watch_out = 0;
45
46   /*
47    * Variables to manipulate the data structure proposed to model the maxmin
48    * fairness. See docummentation for more details.
49    */
50   xbt_swag_t elem_list = NULL;
51   //lmm_element_t elem = NULL;
52   lmm_element_t elem1 = NULL;
53
54
55   xbt_swag_t cnst_list = NULL;
56   //lmm_constraint_t cnst = NULL;
57   lmm_constraint_t cnst1 = NULL;
58   //lmm_constraint_t cnst2 = NULL;
59
60
61   xbt_swag_t var_list = NULL;
62   lmm_variable_t var1 = NULL;
63   lmm_variable_t var2 = NULL;
64
65   /*
66    * Auxiliar variables.
67    */
68   int iteration=0;
69   double mu_partial=0;
70   double lambda_partial=0;
71   double tmp=0;
72   int i,j;
73   FILE *gnuplot_file=NULL;
74   //char print_buf[1024];
75   //char *trace_buf=xbt_malloc0(sizeof(char));
76   //double sum;
77
78
79   DEBUG0("Iterative method configuration snapshot =====>");
80   DEBUG1("#### Maximum number of iterations : %d", max_iterations);
81   DEBUG1("#### Minimum error tolerated      : %e", epsilon_min_error);
82   DEBUG1("#### Step                         : %e", sigma_step);
83
84
85   if ( !(sys->modified))
86     return;
87
88   /* 
89    * Initialize the var list variable with only the active variables. 
90    * Associate an index in the swag variables. Initialize mu.
91    */
92   var_list = &(sys->variable_set);
93   i=0;
94   xbt_swag_foreach(var1, var_list) {
95     if((var1->bound < 0.0) || (var1->weight <= 0.0)){
96       DEBUG1("#### NOTE var1(%d) is a boundless variable", i);
97       var1->mu = -1.0;
98     } else{ 
99       var1->mu = 1.0;
100       var1->new_mu = 2.0;
101     }
102     DEBUG2("#### var1(%d)->mu:  %e", i, var1->mu);
103     DEBUG2("#### var1(%d)->weight: %e", i, var1->weight);
104     i++;
105   }
106
107   /* 
108    * Initialize lambda.
109    */
110   cnst_list=&(sys->active_constraint_set); 
111   xbt_swag_foreach(cnst1, cnst_list) {
112     cnst1->lambda = 1.0;
113     cnst1->new_lambda = 2.0;
114     DEBUG2("#### cnst1(%p)->lambda:  %e", cnst1, cnst1->lambda);
115   }
116
117   if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
118     gnuplot_file = fopen("lambda.in", "w");
119     fprintf(gnuplot_file, "# iteration    lambda1  lambda2 lambda3 ... lambdaP\n");
120   }
121
122   
123   /*
124    * While doesn't reach a minimun error or a number maximum of iterations.
125    */
126   while(overall_error > epsilon_min_error && iteration < max_iterations){
127     iteration++;
128
129     /*                        d Dual
130      * Compute the value of ----------- (\lambda^k, \mu^k) this portion
131      *                       d \mu_i^k
132      * of code depends on function f(x).
133      */
134     var_list = &(sys->variable_set);
135     xbt_swag_foreach(var1, var_list) {
136       mu_partial = 0;
137       if((var1->bound >= 0) && (var1->weight > 0) ){
138         //for each link with capacity cnsts[i] that uses flow of variable var1 do
139         for(i=0; i<var1->cnsts_number; i++)
140           mu_partial += (var1->cnsts[i].constraint)->lambda + var1->mu;
141        
142         mu_partial = -1.0 / mu_partial + var1->bound;
143         var1->new_mu = var1->mu - sigma_step * mu_partial; 
144
145         if(var1->new_mu < 0){
146           var1->new_mu = 0;
147         }
148       }
149     }
150
151
152     /*                         d Dual
153      * Compute the value of ------------- (\lambda^k, \mu^k) this portion
154      *                      d \lambda_i^k
155      * of code depends on function f(x).
156      */
157     j=0;
158     if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
159       fprintf(gnuplot_file, "\n%d",iteration);
160     }
161     xbt_swag_foreach(cnst1, cnst_list) {
162       j++;
163
164       lambda_partial = 0;
165       
166       elem_list = &(cnst1->element_set);
167       watch_out=0;
168       xbt_swag_foreach(elem1, elem_list) {
169    
170         var2 = elem1->variable;
171         
172         if(var2->weight<=0) continue;
173
174         tmp = 0;
175
176         for(i=0; i<var2->cnsts_number; i++){
177           tmp += (var2->cnsts[i].constraint)->lambda;
178         }
179         if(var2->bound > 0)
180           tmp += var2->mu;
181
182         
183         if(tmp==0) break;
184
185         if (tmp==cnst1->lambda)
186           watch_out=1;
187         lambda_partial += (-1.0 / tmp);
188       }
189
190       if(tmp == 0) 
191         cnst1->new_lambda = LAMBDA_STEP;
192       else {
193         lambda_partial += cnst1->bound;
194         if(watch_out && (lambda_partial>0)) {
195           /* INFO6("Watch Out (%d) %p! lambda_partial: %e; lambda : %e ; (%e %e) \n",iteration, cnst1,  */
196           /*            lambda_partial, cnst1->lambda, cnst1->lambda / 2, */
197           /*            cnst1->lambda - sigma_step * lambda_partial); */
198
199           if(cnst1->lambda < 0) WARN2("Value of cnst1->lambda(%p) = %e < 0", cnst1, cnst1->lambda);
200           if((cnst1->lambda - sigma_step * lambda_partial) < 0) WARN1("Value of lambda_new = %e < 0", (cnst1->lambda - sigma_step * lambda_partial));
201
202           if(cnst1->lambda - sigma_step * lambda_partial < cnst1->lambda / 2)
203             cnst1->new_lambda = cnst1->lambda / 2;
204           else
205             cnst1->new_lambda = cnst1->lambda - sigma_step * lambda_partial;
206         } else
207           cnst1->new_lambda = cnst1->lambda - sigma_step * lambda_partial;
208         if(cnst1->new_lambda < 0){
209           cnst1->new_lambda = 0;
210         }
211       }
212
213       if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
214         fprintf(gnuplot_file, "  %e", cnst1->lambda);
215       }
216
217     }
218
219
220     /*
221      * Now computes the values of each variable (\rho) based on
222      * the values of \lambda and \mu.
223      */
224     overall_error=0;
225     xbt_swag_foreach(var1, var_list) {
226       if(var1->weight <=0) 
227         var1->value = 0.0;
228       else {
229         tmp = 0;
230         for(i=0; i<var1->cnsts_number; i++){
231           tmp += (var1->cnsts[i].constraint)->lambda;
232           if(var1->bound > 0) 
233             tmp+=var1->mu;
234         }
235         
236         //computes de overall_error
237         if(overall_error < fabs(var1->value - 1.0/tmp)){
238           overall_error = fabs(var1->value - 1.0/tmp);
239         }
240
241         var1->value = 1.0 / tmp;
242       }
243       
244     }
245
246
247     /* Updating lambda's and mu's */  
248     xbt_swag_foreach(var1, var_list)
249       if(!((var1->bound > 0.0) || (var1->weight <= 0.0)))
250         var1->mu = var1->new_mu;
251     
252     
253     xbt_swag_foreach(cnst1, cnst_list)
254       cnst1->lambda = cnst1->new_lambda;
255   }
256
257
258
259
260   //verify the KKT property
261   xbt_swag_foreach(cnst1, cnst_list){
262     tmp = 0;
263     elem_list = &(cnst1->element_set);
264     xbt_swag_foreach(elem1, elem_list) {
265       var1 = elem1->variable;
266       if(var1->weight<=0) continue;
267       tmp += var1->value;
268     }
269
270     tmp = tmp - cnst1->bound;
271  
272
273     if(tmp != 0 ||  cnst1->lambda != 0){
274       WARN4("The link %s(%p) doesn't match the KKT property, value expected (=0) got (lambda=%e) (sum_rho=%e)", (char *)cnst1->id, cnst1, cnst1->lambda, tmp);
275     }
276     
277   }
278
279     
280   xbt_swag_foreach(var1, var_list){
281     if(var1->bound <= 0 || var1->weight <= 0) continue;
282     tmp = 0;
283     tmp = (var1->value - var1->bound);
284
285     
286     if(tmp != 0 ||  var1->mu != 0){
287       WARN4("The flow %s(%p) doesn't match the KKT property, value expected (=0) got (lambda=%e) (sum_rho=%e)", (char *)var1->id, var1, var1->mu, tmp);
288     }
289
290   }
291
292
293
294
295   if(overall_error <= epsilon_min_error){
296     DEBUG1("The method converge in %d iterations.", iteration);
297   }else{
298     WARN1("Method reach %d iterations, which is the maxmimun number of iterations allowed.", iteration);
299   }
300
301
302   if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
303     fclose(gnuplot_file);
304   }
305
306
307
308
309
310 /*   /\* */
311 /*    * Now computes the values of each variable (\rho) based on */
312 /*    * the values of \lambda and \mu. */
313 /*    *\/ */
314 /*   var_list = &(sys->variable_set); */
315 /*   xbt_swag_foreach(var1, var_list) { */
316 /*     tmp = 0; */
317 /*     for(i=0; i<var1->cnsts_number; i++){ */
318 /*       elem1 = &(var1->cnsts[i]); */
319 /*       tmp += (elem1->constraint)->lambda + var1->mu; */
320 /*     } */
321 /*     var1->weight = 1 / tmp; */
322
323 /*     DEBUG2("var1->weight (id=%s) : %e", (char *)var1->id, var1->weight); */
324 /*   } */
325
326
327
328
329 }