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. */
10 #include "xbt/sysdep.h"
11 #include "xbt/mallocator.h"
12 #include "maxmin_private.h"
20 * SDP specific variables.
23 #include <declarations.h>
26 static void create_cross_link(struct constraintmatrix *myconstraints, int k);
28 static void addentry(struct constraintmatrix *constraints,
29 struct blockmatrix *, int matno, int blkno, int indexi, int indexj, double ent, int blocksize);
32 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
33 "Logging specific to SURF (sdp)");
34 XBT_LOG_NEW_SUBCATEGORY(surf_sdp_out, surf,
35 "Logging specific to SURF (sdp)");
37 ########################################################################
38 ######################## Simple Proportionnal fairness #################
39 ########################################################################
40 ####### Original problem ##########
53 # We assume in the following that n=2^K
55 ####### Standard SDP form #########
57 # A SDP can be written in the following standard form:
59 # (P) min c1*x1+c2*x2+...+cm*xm
60 # st F1*x1+F2*x2+...+Fm*xm-F0=X
63 # Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
64 # constraint X>0 means that X must be positive semidefinite positive.
66 ########## Equivalent SDP #########
70 # y(k,i) for k in [0,K] and i in [1,2^k]
75 # y(2,1) y(2,2) y(2,3) y(2,4)
76 # y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
77 # -------------------------------------------------------
78 # x_1 x_2 x_3 x_4 x_5 x_6 x_7 x_8
81 # Structure Constraints:
83 # forall k in [0,K-1], and i in [1,2^k]
84 # [ y(k, 2i-1) y(k-1, i) ]
85 # [ y(k-1, i ) y(k, 2i) ]
89 # Capacity Constraints:
91 # -(A.y(K,*))_l >= - b_l
93 # Positivity Constraints:
94 # forall k in [0,K], and i in [1,2^k]
97 # Latency Constraints:
98 # forall i in [1,2^k] and v in [0,m-1]
107 //typedef struct lmm_system {
109 // s_xbt_swag_t variable_set; /* a list of lmm_variable_t */
110 // s_xbt_swag_t constraint_set; /* a list of lmm_constraint_t */
111 // s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */
112 // s_xbt_swag_t saturated_variable_set; /* a list of lmm_variable_t */
113 // s_xbt_swag_t saturated_constraint_set; /* a list of lmm_constraint_t_t */
114 // xbt_mallocator_t variable_mallocator;
117 #define get_y(a,b) (pow(2,a)+b-1)
119 void sdp_solve(lmm_system_t sys)
123 * SDP mapping variables.
130 int nb_cnsts_capacity=0;
131 int nb_cnsts_struct=0;
132 int nb_cnsts_positivy=0;
133 int nb_cnsts_latency=0;
136 int total_block_size=0;
138 // FILE *sdpout = fopen("SDPA-printf.tmp","w");
140 double *tempdiag = NULL;
147 * CSDP library specific variables.
149 struct blockmatrix C;
150 struct blockmatrix X,Z;
154 struct constraintmatrix *constraints;
157 * Classic maxmin variables.
159 lmm_constraint_t cnst = NULL;
160 lmm_element_t elem = NULL;
161 xbt_swag_t cnst_list = NULL;
162 xbt_swag_t var_list = NULL;
163 xbt_swag_t elem_list = NULL;
164 lmm_variable_t var = NULL;
167 struct sparseblock *p;
174 if ( !(sys->modified))
178 * Initialize the var list variable with only the active variables.
179 * Associate an index in the swag variables.
181 var_list = &(sys->variable_set);
183 xbt_swag_foreach(var, var_list) {
184 if(var->weight != 0.0) i++;
190 DEBUG1("Variable set : %d", xbt_swag_size(var_list));
191 DEBUG1("Flows : %d", flows);
198 xbt_swag_foreach(var, var_list) {
205 cnst_list=&(sys->active_constraint_set);
206 DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
207 DEBUG1("Links : %d", links);
210 * Those fields are the top level description of the platform furnished in the xml file.
212 links = xbt_swag_size(&(sys->active_constraint_set));
215 * This number is found based on the tree structure explained on top.
217 tmp_k = (double) log((double)flows)/log(2.0);
218 K = (int) ceil(tmp_k);
219 //xbt_assert1(K!=0, "Invalide value of K (=%d) aborting.", K);
223 * The number of variables in the SDP program.
225 nb_var = get_y(K, pow(2,K));
226 DEBUG1("Number of variables in the SDP program : %d", nb_var);
230 * Find the size of each group of constraints.
232 nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
233 nb_cnsts_struct = (int)pow(2,K) - 1;
234 nb_cnsts_positivy = (int)pow(2,K);
235 nb_cnsts_latency = nb_var;
239 * The total number of constraints.
241 nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy + nb_cnsts_latency;
242 CDEBUG1(surf_sdp_out,"Number of constraints = %d", nb_cnsts);
243 DEBUG1("Number of constraints in the SDP program : %d", nb_cnsts);
246 * Keep track of which blocks have off diagonal entries.
248 isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
249 for (i=1; i<=nb_cnsts; i++)
252 C.nblocks = nb_cnsts;
253 C.blocks = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
254 constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
256 for(i = 1; i <= nb_var; i++){
257 constraints[i].blocks=NULL;
260 a = (double *)calloc(nb_var+1, sizeof(double));
269 * For each constraint block do.
271 for(i = 1; i <= nb_cnsts; i++){
274 * Structured blocks are size 2 and all others are size 1.
276 if(i <= nb_cnsts_struct){
277 total_block_size += block_size = 2;
280 total_block_size += block_size = 1;
281 CDEBUG0(surf_sdp_out,"1 ");
285 * All blocks are matrices.
287 C.blocks[block_num].blockcategory = MATRIX;
288 C.blocks[block_num].blocksize = block_size;
289 C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
294 CDEBUG0(surf_sdp_out," ");
298 * Creates de objective function array.
300 a = (double *)calloc((nb_var+1), sizeof(double));
302 for(i = 1; i <= nb_var; i++){
304 //CDEBUG0(surf_sdp_out,"-1 ");
307 //CDEBUG0(surf_sdp_out,"0 ");
314 * Structure contraint blocks.
318 for(k = 1; k <= K; k++){
319 for(i = 1; i <= pow(2,k-1); i++){
320 matno=get_y(k,2*i-1);
321 CDEBUG2(surf_sdp_out,"%d %d 1 1 1", matno , block_num);
322 addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
325 CDEBUG2(surf_sdp_out,"%d %d 2 2 1", matno , block_num);
326 addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
329 CDEBUG2(surf_sdp_out,"%d %d 1 2 1", matno , block_num);
330 addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);
333 CDEBUG2(surf_sdp_out,"%d %d 2 1 1", matno , block_num);
334 addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
336 isdiag[block_num] = 0;
343 * Capacity constraint block.
345 xbt_swag_foreach(cnst, cnst_list) {
347 CDEBUG2(surf_sdp_out,"0 %d 1 1 %f", block_num, - (cnst->bound));
348 addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);
350 elem_list = &(cnst->element_set);
351 xbt_swag_foreach(elem, elem_list){
352 if(elem->variable->weight <=0) break;
353 matno=get_y(K,elem->variable->index);
354 CDEBUG3(surf_sdp_out,"%d %d 1 1 %f", matno, block_num, - (elem->value));
355 addentry(constraints, &C, matno, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
363 * Positivy constraint blocks.
365 for(i = 1; i <= pow(2,K); i++){
367 CDEBUG2(surf_sdp_out,"%d %d 1 1 1", matno, block_num);
368 addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
372 * Latency constraint blocks.
374 xbt_swag_foreach(var, var_list) {
376 if(var->weight && var->bound > 0) {
377 matno = get_y(K,var->index);
378 CDEBUG3(surf_sdp_out,"%d %d 1 1 %f", matno, block_num, var->bound);
379 addentry(constraints, &C, matno, block_num, 1, 1, var->bound, C.blocks[block_num].blocksize);
384 * At this point, we'll stop to recognize whether any of the blocks
385 * are "hidden LP blocks" and correct the block type if needed.
387 for (i=1; i<=nb_cnsts; i++){
388 if ((C.blocks[i].blockcategory != DIAG) &&
389 (isdiag[i]==1) && (C.blocks[i].blocksize > 1)){
391 * We have a hidden diagonal block!
394 blocksz=C.blocks[i].blocksize;
395 tempdiag=(double *)calloc((blocksz+1), sizeof(double));
396 for (j=1; j<=blocksz; j++)
397 tempdiag[j]=C.blocks[i].data.mat[ijtok(j,j,blocksz)];
398 free(C.blocks[i].data.mat);
399 C.blocks[i].data.vec=tempdiag;
400 C.blocks[i].blockcategory=DIAG;
406 * Next, setup issparse and NULL out all nextbyblock pointers.
409 for (i=1; i<=k; i++) {
410 p=constraints[i].blocks;
413 * First, set issparse.
415 if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15)){
421 if (C.blocks[p->blocknum].blockcategory == DIAG)
425 * Setup the cross links.
435 * Create cross link reference.
437 create_cross_link(constraints, nb_var);
441 * Debuging print problem in SDPA format.
443 if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
444 DEBUG0("Printing SDPA...");
445 tmp=strdup("SURF-PROPORTIONNAL.sdpa");
446 write_prob(tmp,total_block_size,nb_var,C,a,constraints);
450 * Initialize parameters.
452 DEBUG0("Initializing solution...");
453 initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
460 DEBUG0("Calling the solver...");
462 stdout=fopen("/dev/null","w");
463 ret = easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
469 case 1: DEBUG0("SUCCESS The problem is primal infeasible");
472 case 2: DEBUG0("SUCCESS The problem is dual infeasible");
475 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.");
478 case 4: DEBUG0("Failure. Maximum number of iterations reached.");
481 case 5: DEBUG0("Failure. Stuck at edge of primal feasibility.");
486 if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
487 tmp=strdup("SURF-PROPORTIONNAL.sol");
488 write_sol(tmp,total_block_size, nb_var, X, y, Z);
492 * Write out the solution if necessary.
494 xbt_swag_foreach(cnst, cnst_list) {
496 elem_list = &(cnst->element_set);
497 xbt_swag_foreach(elem, elem_list) {
498 if(elem->variable->weight <=0) break;
500 i = (int)get_y(K, elem->variable->index);
501 elem->variable->value = y[i];
510 free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
518 if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
526 * Create the cross_link reference in order to have a byblock list.
528 void create_cross_link(struct constraintmatrix *myconstraints, int k){
532 struct sparseblock *p;
533 struct sparseblock *q;
535 struct sparseblock *prev;
541 for (i=1; i<=k; i++){
542 p=myconstraints[i].blocks;
544 if (p->nextbyblock == NULL){
548 * link in the remaining blocks.
550 for (j=i+1; j<=k; j++){
551 q=myconstraints[j].blocks;
554 if (q->blocknum == p->blocknum){
555 if (p->nextbyblock == NULL){
579 void addentry(struct constraintmatrix *constraints,
580 struct blockmatrix *C,
588 struct sparseblock *p;
589 struct sparseblock *p_sav;
591 p=constraints[matno].blocks;
596 * We haven't yet allocated any blocks.
598 p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
600 //two entries because this library ignores indices starting in zerox
601 p->constraintnum=matno;
606 p->entries=calloc(p->numentries+1, sizeof(double));
607 p->iindices=calloc(p->numentries+1, sizeof(int));
608 p->jindices=calloc(p->numentries+1, sizeof(int));
610 p->entries[p->numentries]=ent;
611 p->iindices[p->numentries]=indexi;
612 p->jindices[p->numentries]=indexj;
614 p->blocksize=blocksize;
616 constraints[matno].blocks=p;
619 * We have some existing blocks. See whether this block is already
623 if (p->blocknum == blkno){
625 * Found the right block.
627 p->constraintnum=matno;
629 p->numentries=p->numentries+1;
631 p->entries = realloc(p->entries, (p->numentries+1) * sizeof(double) );
632 p->iindices = realloc(p->iindices, (p->numentries+1) * sizeof(int) );
633 p->jindices = realloc(p->jindices, (p->numentries+1) * sizeof(int) );
635 p->entries[p->numentries]=ent;
636 p->iindices[p->numentries]=indexi;
637 p->jindices[p->numentries]=indexj;
646 * If we get here, we have a non-empty structure but not the right block
647 * inside hence create a new structure.
650 p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
652 //two entries because this library ignores indices starting in zerox
653 p->constraintnum=matno;
658 p->entries=calloc(p->numentries+1, sizeof(double));
659 p->iindices=calloc(p->numentries+1, sizeof(int));
660 p->jindices=calloc(p->numentries+1, sizeof(int));
662 p->entries[p->numentries]=ent;
663 p->iindices[p->numentries]=indexi;
664 p->jindices[p->numentries]=indexj;
666 p->blocksize=blocksize;
672 int blksz=C->blocks[blkno].blocksize;
673 if (C->blocks[blkno].blockcategory == DIAG){
674 C->blocks[blkno].data.vec[indexi]=ent;
676 C->blocks[blkno].data.mat[ijtok(indexi,indexj,blksz)]=ent;
677 C->blocks[blkno].data.mat[ijtok(indexj,indexi,blksz)]=ent;