Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
reindent.
[simgrid.git] / 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 #define NOSHORTS
23 #include <declarations.h>
24
25
26 static void create_cross_link(struct constraintmatrix *myconstraints,
27                               int k);
28
29 static void addentry(struct constraintmatrix *constraints,
30                      struct blockmatrix *, int matno, int blkno,
31                      int indexi, int indexj, double ent, int blocksize);
32
33
34 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
35                                 "Logging specific to SURF (sdp)");
36 XBT_LOG_NEW_SUBCATEGORY(surf_sdp_out, surf,
37                         "Logging specific to SURF (sdp)");
38 /*
39 ########################################################################
40 ######################## Simple Proportionnal fairness #################
41 ########################################################################
42 ####### Original problem ##########
43 #
44 #  Max x_1*x_2*...*x_n
45 #    (A.x)_1 <= b_1
46 #    (A.x)_2 <= b_2
47 #     ...
48 #    (A.x)_m <= b_m
49 #        x_1 <= c_1
50 #        x_2 <= c_2
51 #         ...
52 #        x_m <= c_m
53 #    x>=0
54 #
55 #  We assume in the following that n=2^K
56 #
57 ####### Standard SDP form #########
58 #
59 # A SDP can be written in the following standard form:
60
61 #   (P) min c1*x1+c2*x2+...+cm*xm
62 #       st  F1*x1+F2*x2+...+Fm*xm-F0=X
63 #                                 X >= 0
64 #
65 # Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
66 # constraint X>0 means that X must be positive semidefinite positive.
67 #
68 ########## Equivalent SDP #########
69 #
70 #  Variables:
71 #
72 #    y(k,i) for k in [0,K] and i in [1,2^k]
73 #
74 #  Structure :
75 #                             y(0,1)
76 #                y(1,1)                      y(1,2)
77 #        y(2,1)        y(2,2)        y(2,3)        y(2,4)
78 #    y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
79 #    -------------------------------------------------------
80 #     x_1     x_2    x_3    x_4    x_5    x_6    x_7    x_8
81 #  
82 #
83 #  Structure Constraints:
84 #
85 #    forall k in [0,K-1], and i in [1,2^k]
86 #      [ y(k,  2i-1)   y(k-1, i) ]
87 #      [ y(k-1, i  )   y(k,  2i) ]
88 #
89 #    2^K-1
90 #
91 #  Capacity Constraints:
92 #     forall l in [1,m]
93 #        -(A.y(K,*))_l >= - b_l
94 #
95 #  Positivity Constraints:
96 #    forall k in [0,K], and i in [1,2^k]
97 #        y(k,i) >= 0
98 #
99 #  Latency Constraints:
100 #    forall i in [1,2^k] and v in [0,m-1]
101 #        if(i <= m-1){
102 #            y(k-1, i) <= bound
103 #        }else{
104 #            y(k-1, i) <= 1
105 #        }
106 #  Minimize -y(0,1)
107 */
108
109 //typedef struct lmm_system {
110 //  int modified;
111 //  s_xbt_swag_t variable_set;  /* a list of lmm_variable_t */
112 //  s_xbt_swag_t constraint_set;        /* a list of lmm_constraint_t */
113 //  s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */
114 //  s_xbt_swag_t saturated_variable_set;        /* a list of lmm_variable_t */
115 //  s_xbt_swag_t saturated_constraint_set;      /* a list of lmm_constraint_t_t */
116 //  xbt_mallocator_t variable_mallocator;
117 //} s_lmm_system_t;
118
119 #define get_y(a,b) (pow(2,a)+b-1)
120
121 void sdp_solve(lmm_system_t sys)
122 {
123
124   /*
125    * SDP mapping variables.
126    */
127   int K = 0;
128   int nb_var = 0;
129   int nb_cnsts = 0;
130   int flows = 0;
131   int links = 0;
132   int nb_cnsts_capacity = 0;
133   int nb_cnsts_struct = 0;
134   int nb_cnsts_positivy = 0;
135   int nb_cnsts_latency = 0;
136   int block_num = 0;
137   int block_size;
138   int total_block_size = 0;
139   int *isdiag = NULL;
140   //  FILE *sdpout = fopen("SDPA-printf.tmp","w");
141   int blocksz = 0;
142   double *tempdiag = NULL;
143   int matno = 0;
144   int i = 0;
145   int j = 0;
146   int k = 0;
147
148   /* 
149    * CSDP library specific variables.
150    */
151   struct blockmatrix C;
152   struct blockmatrix X, Z;
153   double *y;
154   double pobj, dobj;
155   double *a;
156   struct constraintmatrix *constraints;
157
158   /*
159    * Classic maxmin variables.
160    */
161   lmm_constraint_t cnst = NULL;
162   lmm_element_t elem = NULL;
163   xbt_swag_t cnst_list = NULL;
164   xbt_swag_t var_list = NULL;
165   xbt_swag_t elem_list = NULL;
166   lmm_variable_t var = NULL;
167
168   double tmp_k;
169   struct sparseblock *p;
170   char *tmp = NULL;
171   FILE *stdout_sav;
172   int ret;
173
174
175
176   if (!(sys->modified))
177     return;
178
179   /* 
180    * Initialize the var list variable with only the active variables. 
181    * Associate an index in the swag variables.
182    */
183   var_list = &(sys->variable_set);
184   i = 0;
185   xbt_swag_foreach(var, var_list) {
186     if (var->weight != 0.0)
187       i++;
188   }
189
190
191
192   flows = i;
193   DEBUG1("Variable set : %d", xbt_swag_size(var_list));
194   DEBUG1("Flows : %d", flows);
195
196   if (flows == 0) {
197     return;
198   }
199
200
201   xbt_swag_foreach(var, var_list) {
202     var->value = 0.0;
203     if (var->weight) {
204       var->index = i--;
205     }
206   }
207
208   cnst_list = &(sys->active_constraint_set);
209   DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
210   DEBUG1("Links : %d", links);
211
212   /*
213    * Those fields are the top level description of the platform furnished in the xml file.
214    */
215   links = xbt_swag_size(&(sys->active_constraint_set));
216
217   /* 
218    * This  number is found based on the tree structure explained on top.
219    */
220   tmp_k = (double) log((double) flows) / log(2.0);
221   K = (int) ceil(tmp_k);
222   //xbt_assert1(K!=0, "Invalide value of K (=%d) aborting.", K);
223
224
225   /* 
226    * The number of variables in the SDP program. 
227    */
228   nb_var = get_y(K, pow(2, K));
229   DEBUG1("Number of variables in the SDP program : %d", nb_var);
230
231
232   /*
233    * Find the size of each group of constraints.
234    */
235   nb_cnsts_capacity = links + ((int) pow(2, K)) - flows;
236   nb_cnsts_struct = (int) pow(2, K) - 1;
237   nb_cnsts_positivy = (int) pow(2, K);
238   nb_cnsts_latency = nb_var;
239
240
241   /*
242    * The total number of constraints.
243    */
244   nb_cnsts =
245       nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy +
246       nb_cnsts_latency;
247   CDEBUG1(surf_sdp_out, "Number of constraints = %d", nb_cnsts);
248   DEBUG1("Number of constraints in the SDP program : %d", nb_cnsts);
249
250   /*
251    * Keep track of which blocks have off diagonal entries. 
252    */
253   isdiag = (int *) calloc((nb_cnsts + 1), sizeof(int));
254   for (i = 1; i <= nb_cnsts; i++)
255     isdiag[i] = 1;
256
257   C.nblocks = nb_cnsts;
258   C.blocks =
259       (struct blockrec *) calloc((C.nblocks + 1), sizeof(struct blockrec));
260   constraints =
261       (struct constraintmatrix *) calloc((nb_var + 1),
262                                          sizeof(struct constraintmatrix));
263
264   for (i = 1; i <= nb_var; i++) {
265     constraints[i].blocks = NULL;
266   }
267
268   a = (double *) calloc(nb_var + 1, sizeof(double));
269
270   /*
271    * Block sizes.
272    */
273   block_num = 1;
274   block_size = 0;
275
276   /*
277    * For each constraint block do.
278    */
279   for (i = 1; i <= nb_cnsts; i++) {
280
281     /*
282      * Structured blocks are size 2 and all others are size 1.
283      */
284     if (i <= nb_cnsts_struct) {
285       total_block_size += block_size = 2;
286       DEBUG0("2 ");
287     } else {
288       total_block_size += block_size = 1;
289       CDEBUG0(surf_sdp_out, "1 ");
290     }
291
292     /*
293      * All blocks are matrices.
294      */
295     C.blocks[block_num].blockcategory = MATRIX;
296     C.blocks[block_num].blocksize = block_size;
297     C.blocks[block_num].data.mat =
298         (double *) calloc(block_size * block_size, sizeof(double));
299
300     block_num++;
301   }
302
303   CDEBUG0(surf_sdp_out, " ");
304
305
306   /*
307    * Creates de objective function array.
308    */
309   a = (double *) calloc((nb_var + 1), sizeof(double));
310
311   for (i = 1; i <= nb_var; i++) {
312     if (get_y(0, 1) == i) {
313       //CDEBUG0(surf_sdp_out,"-1 ");
314       a[i] = -1;
315     } else {
316       //CDEBUG0(surf_sdp_out,"0 ");
317       a[i] = 0;
318     }
319   }
320
321
322   /*
323    * Structure contraint blocks.
324    */
325   block_num = 1;
326   matno = 1;
327   for (k = 1; k <= K; k++) {
328     for (i = 1; i <= pow(2, k - 1); i++) {
329       matno = get_y(k, 2 * i - 1);
330       CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num);
331       addentry(constraints, &C, matno, block_num, 1, 1, 1.0,
332                C.blocks[block_num].blocksize);
333
334       matno = get_y(k, 2 * i);
335       CDEBUG2(surf_sdp_out, "%d %d 2 2 1", matno, block_num);
336       addentry(constraints, &C, matno, block_num, 2, 2, 1.0,
337                C.blocks[block_num].blocksize);
338
339       matno = get_y(k - 1, i);
340       CDEBUG2(surf_sdp_out, "%d %d 1 2 1", matno, block_num);
341       addentry(constraints, &C, matno, block_num, 1, 2, 1.0,
342                C.blocks[block_num].blocksize);
343
344       matno = get_y(k - 1, i);
345       CDEBUG2(surf_sdp_out, "%d %d 2 1 1", matno, block_num);
346       addentry(constraints, &C, matno, block_num, 2, 1, 1.0,
347                C.blocks[block_num].blocksize);
348
349       isdiag[block_num] = 0;
350       block_num++;
351     }
352   }
353
354
355   /*
356    * Capacity constraint block.
357    */
358   xbt_swag_foreach(cnst, cnst_list) {
359
360     CDEBUG2(surf_sdp_out, "0 %d 1 1 %f", block_num, -(cnst->bound));
361     addentry(constraints, &C, 0, block_num, 1, 1, -(cnst->bound),
362              C.blocks[block_num].blocksize);
363
364     elem_list = &(cnst->element_set);
365     xbt_swag_foreach(elem, elem_list) {
366       if (elem->variable->weight <= 0)
367         break;
368       matno = get_y(K, elem->variable->index);
369       CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num,
370               -(elem->value));
371       addentry(constraints, &C, matno, block_num, 1, 1, -(elem->value),
372                C.blocks[block_num].blocksize);
373
374     }
375     block_num++;
376   }
377
378
379   /*
380    * Positivy constraint blocks.
381    */
382   for (i = 1; i <= pow(2, K); i++) {
383     matno = get_y(K, i);
384     CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num);
385     addentry(constraints, &C, matno, block_num, 1, 1, 1.0,
386              C.blocks[block_num].blocksize);
387     block_num++;
388   }
389   /*
390    * Latency constraint blocks.
391    */
392   xbt_swag_foreach(var, var_list) {
393     var->value = 0.0;
394     if (var->weight && var->bound > 0) {
395       matno = get_y(K, var->index);
396       CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num, var->bound);
397       addentry(constraints, &C, matno, block_num, 1, 1, var->bound,
398                C.blocks[block_num].blocksize);
399     }
400   }
401
402   /*
403    * At this point, we'll stop to recognize whether any of the blocks
404    * are "hidden LP blocks"  and correct the block type if needed.
405    */
406   for (i = 1; i <= nb_cnsts; i++) {
407     if ((C.blocks[i].blockcategory != DIAG) &&
408         (isdiag[i] == 1) && (C.blocks[i].blocksize > 1)) {
409       /*
410        * We have a hidden diagonal block!
411        */
412
413       blocksz = C.blocks[i].blocksize;
414       tempdiag = (double *) calloc((blocksz + 1), sizeof(double));
415       for (j = 1; j <= blocksz; j++)
416         tempdiag[j] = C.blocks[i].data.mat[ijtok(j, j, blocksz)];
417       free(C.blocks[i].data.mat);
418       C.blocks[i].data.vec = tempdiag;
419       C.blocks[i].blockcategory = DIAG;
420     };
421   };
422
423
424   /*
425    * Next, setup issparse and NULL out all nextbyblock pointers.
426    */
427   p = NULL;
428   for (i = 1; i <= k; i++) {
429     p = constraints[i].blocks;
430     while (p != NULL) {
431       /*
432        * First, set issparse.
433        */
434       if (((p->numentries) > 0.25 * (p->blocksize))
435           && ((p->numentries) > 15)) {
436         p->issparse = 0;
437       } else {
438         p->issparse = 1;
439       };
440
441       if (C.blocks[p->blocknum].blockcategory == DIAG)
442         p->issparse = 1;
443
444       /*
445        * Setup the cross links.
446        */
447
448       p->nextbyblock = NULL;
449       p = p->next;
450     };
451   };
452
453
454   /*
455    * Create cross link reference.
456    */
457   create_cross_link(constraints, nb_var);
458
459
460   /*
461    * Debuging print problem in SDPA format.
462    */
463   if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
464     DEBUG0("Printing SDPA...");
465     tmp = strdup("SURF-PROPORTIONNAL.sdpa");
466     write_prob(tmp, total_block_size, nb_var, C, a, constraints);
467   }
468
469   /*
470    * Initialize parameters.
471    */
472   DEBUG0("Initializing solution...");
473   initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
474
475
476
477   /*
478    * Call the solver.
479    */
480   DEBUG0("Calling the solver...");
481   stdout_sav = stdout;
482   stdout = fopen("/dev/null", "w");
483   ret =
484       easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y,
485                &Z, &pobj, &dobj);
486   fclose(stdout);
487   stdout = stdout_sav;
488
489   switch (ret) {
490   case 0:
491   case 1:
492     DEBUG0("SUCCESS The problem is primal infeasible");
493     break;
494
495   case 2:
496     DEBUG0("SUCCESS The problem is dual infeasible");
497     break;
498
499   case 3:
500     DEBUG0
501         ("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.");
502     break;
503
504   case 4:
505     DEBUG0("Failure. Maximum number of iterations reached.");
506     break;
507
508   case 5:
509     DEBUG0("Failure. Stuck at edge of primal feasibility.");
510     break;
511
512   }
513
514   if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
515     tmp = strdup("SURF-PROPORTIONNAL.sol");
516     write_sol(tmp, total_block_size, nb_var, X, y, Z);
517   }
518
519   /*
520    * Write out the solution if necessary.
521    */
522   xbt_swag_foreach(cnst, cnst_list) {
523
524     elem_list = &(cnst->element_set);
525     xbt_swag_foreach(elem, elem_list) {
526       if (elem->variable->weight <= 0)
527         break;
528
529       i = (int) get_y(K, elem->variable->index);
530       elem->variable->value = y[i];
531
532     }
533   }
534
535
536   /*
537    * Free up memory.
538    */
539   free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
540
541   free(isdiag);
542   free(tempdiag);
543   free(tmp);
544
545   sys->modified = 0;
546
547   if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
548     lmm_print(sys);
549   }
550
551 }
552
553
554 /*
555  * Create the cross_link reference in order to have a byblock list.
556  */
557 void create_cross_link(struct constraintmatrix *myconstraints, int k)
558 {
559
560   int i, j;
561   int blk;
562   struct sparseblock *p;
563   struct sparseblock *q;
564
565   struct sparseblock *prev;
566
567   /*
568    * Now, cross link.
569    */
570   prev = NULL;
571   for (i = 1; i <= k; i++) {
572     p = myconstraints[i].blocks;
573     while (p != NULL) {
574       if (p->nextbyblock == NULL) {
575         blk = p->blocknum;
576
577         /*
578          * link in the remaining blocks.
579          */
580         for (j = i + 1; j <= k; j++) {
581           q = myconstraints[j].blocks;
582
583           while (q != NULL) {
584             if (q->blocknum == p->blocknum) {
585               if (p->nextbyblock == NULL) {
586                 p->nextbyblock = q;
587                 q->nextbyblock = NULL;
588                 prev = q;
589               } else {
590                 prev->nextbyblock = q;
591                 q->nextbyblock = NULL;
592                 prev = q;
593               }
594               break;
595             }
596             q = q->next;
597           }
598         }
599       }
600       p = p->next;
601     }
602   }
603 }
604
605
606
607
608 void addentry(struct constraintmatrix *constraints,
609               struct blockmatrix *C,
610               int matno,
611               int blkno, int indexi, int indexj, double ent, int blocksize)
612 {
613   struct sparseblock *p;
614   struct sparseblock *p_sav;
615
616   p = constraints[matno].blocks;
617
618   if (matno != 0.0) {
619     if (p == NULL) {
620       /*
621        * We haven't yet allocated any blocks.
622        */
623       p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock));
624
625       //two entries because this library ignores indices starting in zerox
626       p->constraintnum = matno;
627       p->blocknum = blkno;
628       p->numentries = 1;
629       p->next = NULL;
630
631       p->entries = calloc(p->numentries + 1, sizeof(double));
632       p->iindices = calloc(p->numentries + 1, sizeof(int));
633       p->jindices = calloc(p->numentries + 1, sizeof(int));
634
635       p->entries[p->numentries] = ent;
636       p->iindices[p->numentries] = indexi;
637       p->jindices[p->numentries] = indexj;
638
639       p->blocksize = blocksize;
640
641       constraints[matno].blocks = p;
642     } else {
643       /*
644        * We have some existing blocks.  See whether this block is already
645        * in the chain.
646        */
647       while (p != NULL) {
648         if (p->blocknum == blkno) {
649           /*
650            * Found the right block.
651            */
652           p->constraintnum = matno;
653           p->blocknum = blkno;
654           p->numentries = p->numentries + 1;
655
656           p->entries =
657               realloc(p->entries, (p->numentries + 1) * sizeof(double));
658           p->iindices =
659               realloc(p->iindices, (p->numentries + 1) * sizeof(int));
660           p->jindices =
661               realloc(p->jindices, (p->numentries + 1) * sizeof(int));
662
663           p->entries[p->numentries] = ent;
664           p->iindices[p->numentries] = indexi;
665           p->jindices[p->numentries] = indexj;
666
667           return;
668         }
669         p_sav = p;
670         p = p->next;
671       }
672
673       /*
674        * If we get here, we have a non-empty structure but not the right block
675        * inside hence create a new structure.
676        */
677
678       p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock));
679
680       //two entries because this library ignores indices starting in zerox
681       p->constraintnum = matno;
682       p->blocknum = blkno;
683       p->numentries = 1;
684       p->next = NULL;
685
686       p->entries = calloc(p->numentries + 1, sizeof(double));
687       p->iindices = calloc(p->numentries + 1, sizeof(int));
688       p->jindices = calloc(p->numentries + 1, sizeof(int));
689
690       p->entries[p->numentries] = ent;
691       p->iindices[p->numentries] = indexi;
692       p->jindices[p->numentries] = indexj;
693
694       p->blocksize = blocksize;
695
696       p_sav->next = p;
697     }
698   } else {
699     if (ent != 0.0) {
700       int blksz = C->blocks[blkno].blocksize;
701       if (C->blocks[blkno].blockcategory == DIAG) {
702         C->blocks[blkno].data.vec[indexi] = ent;
703       } else {
704         C->blocks[blkno].data.mat[ijtok(indexi, indexj, blksz)] = ent;
705         C->blocks[blkno].data.mat[ijtok(indexj, indexi, blksz)] = ent;
706       };
707     };
708
709   }
710 }