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.
22 #include <declarations.h>
25 static void create_cross_link(struct constraintmatrix *myconstraints, int k);
27 static void addentry(struct constraintmatrix *constraints,
28 struct blockmatrix *, int matno, int blkno, int indexi, int indexj, double ent, int blocksize);
31 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
32 "Logging specific to SURF (sdp)");
34 ########################################################################
35 ######################## Simple Proportionnal fairness #################
36 ########################################################################
37 ####### Original problem ##########
46 # We assume in the following that n=2^K
48 ####### Standard SDP form #########
50 # A SDP can be written in the following standard form:
52 # (P) min c1*x1+c2*x2+...+cm*xm
53 # st F1*x1+F2*x2+...+Fm*xm-F0=X
56 # Where F1, F2, ..., Fm are symetric matrixes of size n by n. The
57 # constraint X>0 means that X must be positive semidefinite positive.
59 ########## Equivalent SDP #########
63 # y(k,i) for k in [0,K] and i in [1,2^k]
68 # y(2,1) y(2,2) y(2,3) y(2,4)
69 # y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8)
70 # -------------------------------------------------------
71 # x_1 x_2 x_3 x_4 x_5 x_6 x_7 x_8
74 # Structure Constraints:
76 # forall k in [0,K-1], and i in [1,2^k]
77 # [ y(k, 2i-1) y(k-1, i) ]
78 # [ y(k-1, i ) y(k, 2i) ]
82 # Capacity Constraints:
84 # -(A.y(K,*))_l >= - b_l
86 # Positivity Constraints:
87 # forall k in [0,K], and i in [1,2^k]
93 //typedef struct lmm_system {
95 // s_xbt_swag_t variable_set; /* a list of lmm_variable_t */
96 // s_xbt_swag_t constraint_set; /* a list of lmm_constraint_t */
97 // s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */
98 // s_xbt_swag_t saturated_variable_set; /* a list of lmm_variable_t */
99 // s_xbt_swag_t saturated_constraint_set; /* a list of lmm_constraint_t_t */
100 // xbt_mallocator_t variable_mallocator;
103 #define get_y(a,b) (pow(2,a)+b-1)
105 void sdp_solve(lmm_system_t sys)
107 lmm_variable_t var = NULL;
110 * SDP mapping variables.
117 int nb_cnsts_capacity=0;
118 int nb_cnsts_struct=0;
119 int nb_cnsts_positivy=0;
122 int total_block_size=0;
124 FILE *sdpout = fopen("SDPA-printf.tmp","w");
126 double *tempdiag = NULL;
133 * CSDP library specific variables.
135 struct blockmatrix C;
136 struct blockmatrix X,Z;
140 struct constraintmatrix *constraints;
143 * Classic maxmin variables.
145 lmm_constraint_t cnst = NULL;
146 lmm_element_t elem = NULL;
147 xbt_swag_t cnst_list = NULL;
148 xbt_swag_t var_list = NULL;
149 xbt_swag_t elem_list = NULL;
151 if ( !(sys->modified))
155 * Initialize the var list variable with only the active variables.
156 * Associate an index in the swag variables.
159 var_list = &(sys->variable_set);
160 DEBUG1("Variable set : %d", xbt_swag_size(var_list));
162 xbt_swag_foreach(var, var_list) {
164 if(var->weight) var->index = i++;
167 cnst_list=&(sys->saturated_constraint_set);
168 DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
171 * Those fields are the top level description of the platform furnished in the xml file.
174 links = xbt_swag_size(&(sys->active_constraint_set));
177 * This number is found based on the tree structure explained on top.
181 tmp_k = (double) log((double)flows)/log(2.0);
182 K = (int) ceil(tmp_k);
186 * The number of variables in the SDP program.
188 nb_var = get_y(K, pow(2,K));
189 fprintf(sdpout,"%d\n", nb_var);
192 * Find the size of each group of constraints.
194 nb_cnsts_capacity = links + ((int)pow(2,K)) - flows;
195 nb_cnsts_struct = (int)pow(2,K) - 1;
196 nb_cnsts_positivy = (int)pow(2,K);
199 * The total number of constraints.
201 nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
202 fprintf(sdpout,"%d\n", nb_cnsts);
205 * Keep track of which blocks have off diagonal entries.
207 isdiag=(int *)calloc((nb_cnsts+1), sizeof(int));
208 for (i=1; i<=nb_cnsts; i++)
211 C.nblocks = nb_cnsts;
212 C.blocks = (struct blockrec *) calloc((C.nblocks+1),sizeof(struct blockrec));
213 constraints = (struct constraintmatrix *)calloc((nb_var+1),sizeof(struct constraintmatrix));
215 for(i = 1; i <= nb_var; i++){
216 constraints[i].blocks=NULL;
219 a = (double *)calloc(nb_var+1, sizeof(double));
228 * For each constraint block do.
230 for(i = 1; i <= nb_cnsts; i++){
233 * Structured blocks are size 2 and all others are size 1.
235 if(i <= nb_cnsts_struct){
236 total_block_size += block_size = 2;
237 fprintf(sdpout,"2 ");
239 total_block_size += block_size = 1;
240 fprintf(sdpout,"1 ");
244 * All blocks are matrices.
246 C.blocks[block_num].blockcategory = MATRIX;
247 C.blocks[block_num].blocksize = block_size;
248 C.blocks[block_num].data.mat = (double *)calloc(block_size * block_size, sizeof(double));
253 fprintf(sdpout,"\n");
257 * Creates de objective function array.
259 a = (double *)calloc((nb_var+1), sizeof(double));
261 for(i = 1; i <= nb_var; i++){
263 fprintf(sdpout,"-1 ");
266 fprintf(sdpout,"0 ");
270 fprintf(sdpout,"\n");
274 * Structure contraint blocks.
278 for(k = 1; k <= K; k++){
279 for(i = 1; i <= pow(2,k-1); i++){
280 matno=get_y(k,2*i-1);
281 fprintf(sdpout,"%d %d 1 1 1\n", matno , block_num);
282 addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
285 fprintf(sdpout,"%d %d 2 2 1\n", matno , block_num);
286 addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize);
289 fprintf(sdpout,"%d %d 1 2 1\n", matno , block_num);
290 addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize);
293 fprintf(sdpout,"%d %d 2 1 1\n", matno , block_num);
294 addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
296 isdiag[block_num] = 0;
303 * Capacity constraint block.
305 xbt_swag_foreach(cnst, cnst_list) {
307 fprintf(sdpout,"0 %d 1 1 %d\n", block_num, (int) - (cnst->bound));
308 addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);
311 elem_list = &(cnst->element_set);
312 xbt_swag_foreach(elem, elem_list) {
313 if(elem->variable->weight <=0) break;
314 fprintf(sdpout,"%d %d 1 1 %d\n", elem->variable->index, block_num, (int) - (elem->variable->value));
315 addentry(constraints, &C, elem->variable->index, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
323 * Positivy constraint blocks.
325 for(i = 1; i <= pow(2,K); i++){
327 fprintf(sdpout,"%d %d 1 1 1\n", matno, block_num);
328 addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
334 * At this point, we'll stop to recognize whether any of the blocks
335 * are "hidden LP blocks" and correct the block type if needed.
337 for (i=1; i<=nb_cnsts; i++){
338 if ((C.blocks[i].blockcategory != DIAG) &&
339 (isdiag[i]==1) && (C.blocks[i].blocksize > 1)){
341 * We have a hidden diagonal block!
344 printf("Block %d is actually diagonal.\n",i);
346 blocksz=C.blocks[i].blocksize;
347 tempdiag=(double *)calloc((blocksz+1), sizeof(double));
348 for (j=1; j<=blocksz; j++)
349 tempdiag[j]=C.blocks[i].data.mat[ijtok(j,j,blocksz)];
350 free(C.blocks[i].data.mat);
351 C.blocks[i].data.vec=tempdiag;
352 C.blocks[i].blockcategory=DIAG;
358 * Next, setup issparse and NULL out all nextbyblock pointers.
360 struct sparseblock *p=NULL;
361 for (i=1; i<=k; i++) {
362 p=constraints[i].blocks;
365 * First, set issparse.
367 if (((p->numentries) > 0.25*(p->blocksize)) && ((p->numentries) > 15)){
373 if (C.blocks[p->blocknum].blockcategory == DIAG)
377 * Setup the cross links.
387 * Create cross link reference.
389 create_cross_link(constraints, nb_var);
393 * Debuging print problem in SDPA format.
395 printf("Printing SDPA...\n");
396 if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
397 char *tmp=strdup("SDPA.tmp");
398 write_prob(tmp, nb_cnsts, nb_var, C, a, constraints);
399 //int write_prob(char *fname, int n, int k, struct blockmatrix C, double *a, struct constraintmatrix *constraints);
404 * Initialize parameters.
406 printf("Initializing solution...\n");
407 initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
414 printf("Calling the solver...\n");
415 int ret = easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
419 case 1: printf("SUCCESS The problem is primal infeasible\n");
422 case 2: printf("SUCCESS The problem is dual infeasible\n");
425 case 3: printf("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");
428 case 4: printf("Failure. Maximum number of iterations reached.");
431 case 5: printf("Failure. Stuck at edge of primal feasibility.");
438 * Write out the solution if necessary.
440 xbt_swag_foreach(cnst, cnst_list) {
442 elem_list = &(cnst->element_set);
443 xbt_swag_foreach(elem, elem_list) {
444 if(elem->variable->weight <=0) break;
446 i = (int)get_y(K, elem->variable->index);
447 elem->variable->value = y[i];
456 free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
462 if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
470 * Create the cross_link reference in order to have a byblock list.
472 void create_cross_link(struct constraintmatrix *myconstraints, int k){
476 struct sparseblock *p;
477 struct sparseblock *q;
479 struct sparseblock *prev;
485 for (i=1; i<=k; i++){
486 p=myconstraints[i].blocks;
488 if (p->nextbyblock == NULL){
492 * link in the remaining blocks.
494 for (j=i+1; j<=k; j++){
495 q=myconstraints[j].blocks;
498 if (q->blocknum == p->blocknum){
499 if (p->nextbyblock == NULL){
523 void addentry(struct constraintmatrix *constraints,
524 struct blockmatrix *C,
532 struct sparseblock *p;
533 struct sparseblock *p_sav;
535 p=constraints[matno].blocks;
540 * We haven't yet allocated any blocks.
542 p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
544 //two entries because this library ignores indices starting in zerox
545 p->constraintnum=matno;
550 p->entries=calloc(p->numentries+1, sizeof(double));
551 p->iindices=calloc(p->numentries+1, sizeof(int));
552 p->jindices=calloc(p->numentries+1, sizeof(int));
554 p->entries[p->numentries]=ent;
555 p->iindices[p->numentries]=indexi;
556 p->jindices[p->numentries]=indexj;
558 p->blocksize=blocksize;
560 constraints[matno].blocks=p;
563 * We have some existing blocks. See whether this block is already
567 if (p->blocknum == blkno){
569 * Found the right block.
571 p->constraintnum=matno;
573 p->numentries=p->numentries+1;
575 p->entries = realloc(p->entries, (p->numentries+1) * sizeof(double) );
576 p->iindices = realloc(p->iindices, (p->numentries+1) * sizeof(int) );
577 p->jindices = realloc(p->jindices, (p->numentries+1) * sizeof(int) );
579 p->entries[p->numentries]=ent;
580 p->iindices[p->numentries]=indexi;
581 p->jindices[p->numentries]=indexj;
590 * If we get here, we have a non-empty structure but not the right block
591 * inside hence create a new structure.
594 p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
596 //two entries because this library ignores indices starting in zerox
597 p->constraintnum=matno;
602 p->entries=calloc(p->numentries+1, sizeof(double));
603 p->iindices=calloc(p->numentries+1, sizeof(int));
604 p->jindices=calloc(p->numentries+1, sizeof(int));
606 p->entries[p->numentries]=ent;
607 p->iindices[p->numentries]=indexi;
608 p->jindices[p->numentries]=indexj;
610 p->blocksize=blocksize;
616 int blksz=C->blocks[blkno].blocksize;
617 if (C->blocks[blkno].blockcategory == DIAG){
618 C->blocks[blkno].data.vec[indexi]=ent;
620 C->blocks[blkno].data.mat[ijtok(indexi,indexj,blksz)]=ent;
621 C->blocks[blkno].data.mat[ijtok(indexj,indexi,blksz)]=ent;