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 "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,
30 int indexi, int indexj, double ent, int blocksize);
33 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf,
34 "Logging specific to SURF (sdp)");
35 XBT_LOG_NEW_SUBCATEGORY(surf_sdp_out, surf, "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)
191 DEBUG1("Variable set : %d", xbt_swag_size(var_list));
192 DEBUG1("Flows : %d", flows);
199 xbt_swag_foreach(var, var_list) {
206 cnst_list = &(sys->active_constraint_set);
207 DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list));
208 DEBUG1("Links : %d", links);
211 * Those fields are the top level description of the platform furnished in the xml file.
213 links = xbt_swag_size(&(sys->active_constraint_set));
216 * This number is found based on the tree structure explained on top.
218 tmp_k = (double) log((double) flows) / log(2.0);
219 K = (int) ceil(tmp_k);
220 //xbt_assert1(K!=0, "Invalide value of K (=%d) aborting.", K);
224 * The number of variables in the SDP program.
226 nb_var = get_y(K, pow(2, K));
227 DEBUG1("Number of variables in the SDP program : %d", nb_var);
231 * Find the size of each group of constraints.
233 nb_cnsts_capacity = links + ((int) pow(2, K)) - flows;
234 nb_cnsts_struct = (int) pow(2, K) - 1;
235 nb_cnsts_positivy = (int) pow(2, K);
236 nb_cnsts_latency = nb_var;
240 * The total number of constraints.
243 nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy +
245 CDEBUG1(surf_sdp_out, "Number of constraints = %d", nb_cnsts);
246 DEBUG1("Number of constraints in the SDP program : %d", nb_cnsts);
249 * Keep track of which blocks have off diagonal entries.
251 isdiag = (int *) calloc((nb_cnsts + 1), sizeof(int));
252 for (i = 1; i <= nb_cnsts; i++)
255 C.nblocks = nb_cnsts;
257 (struct blockrec *) calloc((C.nblocks + 1), sizeof(struct blockrec));
259 (struct constraintmatrix *) calloc((nb_var + 1),
260 sizeof(struct constraintmatrix));
262 for (i = 1; i <= nb_var; i++) {
263 constraints[i].blocks = NULL;
266 a = (double *) calloc(nb_var + 1, sizeof(double));
275 * For each constraint block do.
277 for (i = 1; i <= nb_cnsts; i++) {
280 * Structured blocks are size 2 and all others are size 1.
282 if (i <= nb_cnsts_struct) {
283 total_block_size += block_size = 2;
286 total_block_size += block_size = 1;
287 CDEBUG0(surf_sdp_out, "1 ");
291 * All blocks are matrices.
293 C.blocks[block_num].blockcategory = MATRIX;
294 C.blocks[block_num].blocksize = block_size;
295 C.blocks[block_num].data.mat =
296 (double *) calloc(block_size * block_size, sizeof(double));
301 CDEBUG0(surf_sdp_out, " ");
305 * Creates de objective function array.
307 a = (double *) calloc((nb_var + 1), sizeof(double));
309 for (i = 1; i <= nb_var; i++) {
310 if (get_y(0, 1) == i) {
311 //CDEBUG0(surf_sdp_out,"-1 ");
314 //CDEBUG0(surf_sdp_out,"0 ");
321 * Structure contraint blocks.
325 for (k = 1; k <= K; k++) {
326 for (i = 1; i <= pow(2, k - 1); i++) {
327 matno = get_y(k, 2 * i - 1);
328 CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num);
329 addentry(constraints, &C, matno, block_num, 1, 1, 1.0,
330 C.blocks[block_num].blocksize);
332 matno = get_y(k, 2 * i);
333 CDEBUG2(surf_sdp_out, "%d %d 2 2 1", matno, block_num);
334 addentry(constraints, &C, matno, block_num, 2, 2, 1.0,
335 C.blocks[block_num].blocksize);
337 matno = get_y(k - 1, i);
338 CDEBUG2(surf_sdp_out, "%d %d 1 2 1", matno, block_num);
339 addentry(constraints, &C, matno, block_num, 1, 2, 1.0,
340 C.blocks[block_num].blocksize);
342 matno = get_y(k - 1, i);
343 CDEBUG2(surf_sdp_out, "%d %d 2 1 1", matno, block_num);
344 addentry(constraints, &C, matno, block_num, 2, 1, 1.0,
345 C.blocks[block_num].blocksize);
347 isdiag[block_num] = 0;
354 * Capacity constraint block.
356 xbt_swag_foreach(cnst, cnst_list) {
358 CDEBUG2(surf_sdp_out, "0 %d 1 1 %f", block_num, -(cnst->bound));
359 addentry(constraints, &C, 0, block_num, 1, 1, -(cnst->bound),
360 C.blocks[block_num].blocksize);
362 elem_list = &(cnst->element_set);
363 xbt_swag_foreach(elem, elem_list) {
364 if (elem->variable->weight <= 0)
366 matno = get_y(K, elem->variable->index);
367 CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num, -(elem->value));
368 addentry(constraints, &C, matno, block_num, 1, 1, -(elem->value),
369 C.blocks[block_num].blocksize);
377 * Positivy constraint blocks.
379 for (i = 1; i <= pow(2, K); i++) {
381 CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num);
382 addentry(constraints, &C, matno, block_num, 1, 1, 1.0,
383 C.blocks[block_num].blocksize);
387 * Latency constraint blocks.
389 xbt_swag_foreach(var, var_list) {
391 if (var->weight && var->bound > 0) {
392 matno = get_y(K, var->index);
393 CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num, var->bound);
394 addentry(constraints, &C, matno, block_num, 1, 1, var->bound,
395 C.blocks[block_num].blocksize);
400 * At this point, we'll stop to recognize whether any of the blocks
401 * are "hidden LP blocks" and correct the block type if needed.
403 for (i = 1; i <= nb_cnsts; i++) {
404 if ((C.blocks[i].blockcategory != DIAG) &&
405 (isdiag[i] == 1) && (C.blocks[i].blocksize > 1)) {
407 * We have a hidden diagonal block!
410 blocksz = C.blocks[i].blocksize;
411 tempdiag = (double *) calloc((blocksz + 1), sizeof(double));
412 for (j = 1; j <= blocksz; j++)
413 tempdiag[j] = C.blocks[i].data.mat[ijtok(j, j, blocksz)];
414 free(C.blocks[i].data.mat);
415 C.blocks[i].data.vec = tempdiag;
416 C.blocks[i].blockcategory = DIAG;
422 * Next, setup issparse and NULL out all nextbyblock pointers.
425 for (i = 1; i <= k; i++) {
426 p = constraints[i].blocks;
429 * First, set issparse.
431 if (((p->numentries) > 0.25 * (p->blocksize))
432 && ((p->numentries) > 15)) {
438 if (C.blocks[p->blocknum].blockcategory == DIAG)
442 * Setup the cross links.
445 p->nextbyblock = NULL;
452 * Create cross link reference.
454 create_cross_link(constraints, nb_var);
458 * Debuging print problem in SDPA format.
460 if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
461 DEBUG0("Printing SDPA...");
462 tmp = strdup("SURF-PROPORTIONNAL.sdpa");
463 write_prob(tmp, total_block_size, nb_var, C, a, constraints);
467 * Initialize parameters.
469 DEBUG0("Initializing solution...");
470 initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
477 DEBUG0("Calling the solver...");
479 stdout = fopen("/dev/null", "w");
481 easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y,
489 DEBUG0("SUCCESS The problem is primal infeasible");
493 DEBUG0("SUCCESS The problem is dual infeasible");
498 ("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 DEBUG0("Failure. Maximum number of iterations reached.");
506 DEBUG0("Failure. Stuck at edge of primal feasibility.");
511 if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
512 tmp = strdup("SURF-PROPORTIONNAL.sol");
513 write_sol(tmp, total_block_size, nb_var, X, y, Z);
517 * Write out the solution if necessary.
519 xbt_swag_foreach(cnst, cnst_list) {
521 elem_list = &(cnst->element_set);
522 xbt_swag_foreach(elem, elem_list) {
523 if (elem->variable->weight <= 0)
526 i = (int) get_y(K, elem->variable->index);
527 elem->variable->value = y[i];
536 free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z);
544 if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
552 * Create the cross_link reference in order to have a byblock list.
554 void create_cross_link(struct constraintmatrix *myconstraints, int k)
559 struct sparseblock *p;
560 struct sparseblock *q;
562 struct sparseblock *prev;
568 for (i = 1; i <= k; i++) {
569 p = myconstraints[i].blocks;
571 if (p->nextbyblock == NULL) {
575 * link in the remaining blocks.
577 for (j = i + 1; j <= k; j++) {
578 q = myconstraints[j].blocks;
581 if (q->blocknum == p->blocknum) {
582 if (p->nextbyblock == NULL) {
584 q->nextbyblock = NULL;
587 prev->nextbyblock = q;
588 q->nextbyblock = NULL;
605 void addentry(struct constraintmatrix *constraints,
606 struct blockmatrix *C,
608 int blkno, int indexi, int indexj, double ent, int blocksize)
610 struct sparseblock *p;
611 struct sparseblock *p_sav;
613 p = constraints[matno].blocks;
618 * We haven't yet allocated any blocks.
620 p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock));
622 //two entries because this library ignores indices starting in zerox
623 p->constraintnum = matno;
628 p->entries = calloc(p->numentries + 1, sizeof(double));
629 p->iindices = calloc(p->numentries + 1, sizeof(int));
630 p->jindices = calloc(p->numentries + 1, sizeof(int));
632 p->entries[p->numentries] = ent;
633 p->iindices[p->numentries] = indexi;
634 p->jindices[p->numentries] = indexj;
636 p->blocksize = blocksize;
638 constraints[matno].blocks = p;
641 * We have some existing blocks. See whether this block is already
645 if (p->blocknum == blkno) {
647 * Found the right block.
649 p->constraintnum = matno;
651 p->numentries = p->numentries + 1;
654 realloc(p->entries, (p->numentries + 1) * sizeof(double));
656 realloc(p->iindices, (p->numentries + 1) * sizeof(int));
658 realloc(p->jindices, (p->numentries + 1) * sizeof(int));
660 p->entries[p->numentries] = ent;
661 p->iindices[p->numentries] = indexi;
662 p->jindices[p->numentries] = indexj;
671 * If we get here, we have a non-empty structure but not the right block
672 * inside hence create a new structure.
675 p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock));
677 //two entries because this library ignores indices starting in zerox
678 p->constraintnum = matno;
683 p->entries = calloc(p->numentries + 1, sizeof(double));
684 p->iindices = calloc(p->numentries + 1, sizeof(int));
685 p->jindices = calloc(p->numentries + 1, sizeof(int));
687 p->entries[p->numentries] = ent;
688 p->iindices[p->numentries] = indexi;
689 p->jindices[p->numentries] = indexj;
691 p->blocksize = blocksize;
697 int blksz = C->blocks[blkno].blocksize;
698 if (C->blocks[blkno].blockcategory == DIAG) {
699 C->blocks[blkno].data.vec[indexi] = ent;
701 C->blocks[blkno].data.mat[ijtok(indexi, indexj, blksz)] = ent;
702 C->blocks[blkno].data.mat[ijtok(indexj, indexi, blksz)] = ent;