X-Git-Url: http://bilbo.iut-bm.univ-fcomte.fr/pub/gitweb/simgrid.git/blobdiff_plain/ddd401baf17993eeac40c56d6082066777535c78..18bbef07996d903b23d1d769ca16d412cee1d5b2:/src/surf/sdp.c diff --git a/src/surf/sdp.c b/src/surf/sdp.c deleted file mode 100644 index ed543c6487..0000000000 --- a/src/surf/sdp.c +++ /dev/null @@ -1,707 +0,0 @@ -/* $Id$ */ - -/* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved. */ - -/* This program is free software; you can redistribute it and/or modify it - * under the terms of the license (GNU LGPL) which comes with this package. */ - - -#include "xbt/log.h" -#include "xbt/sysdep.h" -#include "maxmin_private.h" - -#include -#include -#ifndef MATH -#include -#endif - -/* - * SDP specific variables. - */ -#define NOSHORTS -#include - - -static void create_cross_link(struct constraintmatrix *myconstraints, int k); - -static void addentry(struct constraintmatrix *constraints, - struct blockmatrix *, int matno, int blkno, - int indexi, int indexj, double ent, int blocksize); - - -XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf, - "Logging specific to SURF (sdp)"); -XBT_LOG_NEW_SUBCATEGORY(surf_sdp_out, surf, "Logging specific to SURF (sdp)"); -/* -######################################################################## -######################## Simple Proportionnal fairness ################# -######################################################################## -####### Original problem ########## -# -# Max x_1*x_2*...*x_n -# (A.x)_1 <= b_1 -# (A.x)_2 <= b_2 -# ... -# (A.x)_m <= b_m -# x_1 <= c_1 -# x_2 <= c_2 -# ... -# x_m <= c_m -# x>=0 -# -# We assume in the following that n=2^K -# -####### Standard SDP form ######### -# -# A SDP can be written in the following standard form: -# -# (P) min c1*x1+c2*x2+...+cm*xm -# st F1*x1+F2*x2+...+Fm*xm-F0=X -# X >= 0 -# -# Where F1, F2, ..., Fm are symetric matrixes of size n by n. The -# constraint X>0 means that X must be positive semidefinite positive. -# -########## Equivalent SDP ######### -# -# Variables: -# -# y(k,i) for k in [0,K] and i in [1,2^k] -# -# Structure : -# y(0,1) -# y(1,1) y(1,2) -# y(2,1) y(2,2) y(2,3) y(2,4) -# y(3,1) y(3,2) y(3,3) y(3,4) y(3,5) y(3,6) y(3,7) y(3,8) -# ------------------------------------------------------- -# x_1 x_2 x_3 x_4 x_5 x_6 x_7 x_8 -# -# -# Structure Constraints: -# -# forall k in [0,K-1], and i in [1,2^k] -# [ y(k, 2i-1) y(k-1, i) ] -# [ y(k-1, i ) y(k, 2i) ] -# -# 2^K-1 -# -# Capacity Constraints: -# forall l in [1,m] -# -(A.y(K,*))_l >= - b_l -# -# Positivity Constraints: -# forall k in [0,K], and i in [1,2^k] -# y(k,i) >= 0 -# -# Latency Constraints: -# forall i in [1,2^k] and v in [0,m-1] -# if(i <= m-1){ -# y(k-1, i) <= bound -# }else{ -# y(k-1, i) <= 1 -# } -# Minimize -y(0,1) -*/ - -//typedef struct lmm_system { -// int modified; -// s_xbt_swag_t variable_set; /* a list of lmm_variable_t */ -// s_xbt_swag_t constraint_set; /* a list of lmm_constraint_t */ -// s_xbt_swag_t active_constraint_set; /* a list of lmm_constraint_t */ -// s_xbt_swag_t saturated_variable_set; /* a list of lmm_variable_t */ -// s_xbt_swag_t saturated_constraint_set; /* a list of lmm_constraint_t_t */ -// xbt_mallocator_t variable_mallocator; -//} s_lmm_system_t; - -#define get_y(a,b) (pow(2,a)+b-1) - -void sdp_solve(lmm_system_t sys) -{ - - /* - * SDP mapping variables. - */ - int K = 0; - int nb_var = 0; - int nb_cnsts = 0; - int flows = 0; - int links = 0; - int nb_cnsts_capacity = 0; - int nb_cnsts_struct = 0; - int nb_cnsts_positivy = 0; - int nb_cnsts_latency = 0; - int block_num = 0; - int block_size; - int total_block_size = 0; - int *isdiag = NULL; - // FILE *sdpout = fopen("SDPA-printf.tmp","w"); - int blocksz = 0; - double *tempdiag = NULL; - int matno = 0; - int i = 0; - int j = 0; - int k = 0; - - /* - * CSDP library specific variables. - */ - struct blockmatrix C; - struct blockmatrix X, Z; - double *y; - double pobj, dobj; - double *a; - struct constraintmatrix *constraints; - - /* - * Classic maxmin variables. - */ - lmm_constraint_t cnst = NULL; - lmm_element_t elem = NULL; - xbt_swag_t cnst_list = NULL; - xbt_swag_t var_list = NULL; - xbt_swag_t elem_list = NULL; - lmm_variable_t var = NULL; - - double tmp_k; - struct sparseblock *p; - char *tmp = NULL; - FILE *stdout_sav; - int ret; - - - - if (!(sys->modified)) - return; - - /* - * Initialize the var list variable with only the active variables. - * Associate an index in the swag variables. - */ - var_list = &(sys->variable_set); - i = 0; - xbt_swag_foreach(var, var_list) { - if (var->weight != 0.0) - i++; - } - - - - flows = i; - DEBUG1("Variable set : %d", xbt_swag_size(var_list)); - DEBUG1("Flows : %d", flows); - - if (flows == 0) { - return; - } - - - xbt_swag_foreach(var, var_list) { - var->value = 0.0; - if (var->weight) { - var->index = i--; - } - } - - cnst_list = &(sys->active_constraint_set); - DEBUG1("Active constraints : %d", xbt_swag_size(cnst_list)); - DEBUG1("Links : %d", links); - - /* - * Those fields are the top level description of the platform furnished in the xml file. - */ - links = xbt_swag_size(&(sys->active_constraint_set)); - - /* - * This number is found based on the tree structure explained on top. - */ - tmp_k = (double) log((double) flows) / log(2.0); - K = (int) ceil(tmp_k); - //xbt_assert1(K!=0, "Invalide value of K (=%d) aborting.", K); - - - /* - * The number of variables in the SDP program. - */ - nb_var = get_y(K, pow(2, K)); - DEBUG1("Number of variables in the SDP program : %d", nb_var); - - - /* - * Find the size of each group of constraints. - */ - nb_cnsts_capacity = links + ((int) pow(2, K)) - flows; - nb_cnsts_struct = (int) pow(2, K) - 1; - nb_cnsts_positivy = (int) pow(2, K); - nb_cnsts_latency = nb_var; - - - /* - * The total number of constraints. - */ - nb_cnsts = - nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy + - nb_cnsts_latency; - CDEBUG1(surf_sdp_out, "Number of constraints = %d", nb_cnsts); - DEBUG1("Number of constraints in the SDP program : %d", nb_cnsts); - - /* - * Keep track of which blocks have off diagonal entries. - */ - isdiag = (int *) calloc((nb_cnsts + 1), sizeof(int)); - for (i = 1; i <= nb_cnsts; i++) - isdiag[i] = 1; - - C.nblocks = nb_cnsts; - C.blocks = - (struct blockrec *) calloc((C.nblocks + 1), sizeof(struct blockrec)); - constraints = - (struct constraintmatrix *) calloc((nb_var + 1), - sizeof(struct constraintmatrix)); - - for (i = 1; i <= nb_var; i++) { - constraints[i].blocks = NULL; - } - - a = (double *) calloc(nb_var + 1, sizeof(double)); - - /* - * Block sizes. - */ - block_num = 1; - block_size = 0; - - /* - * For each constraint block do. - */ - for (i = 1; i <= nb_cnsts; i++) { - - /* - * Structured blocks are size 2 and all others are size 1. - */ - if (i <= nb_cnsts_struct) { - total_block_size += block_size = 2; - DEBUG0("2 "); - } else { - total_block_size += block_size = 1; - CDEBUG0(surf_sdp_out, "1 "); - } - - /* - * All blocks are matrices. - */ - C.blocks[block_num].blockcategory = MATRIX; - C.blocks[block_num].blocksize = block_size; - C.blocks[block_num].data.mat = - (double *) calloc(block_size * block_size, sizeof(double)); - - block_num++; - } - - CDEBUG0(surf_sdp_out, " "); - - - /* - * Creates de objective function array. - */ - a = (double *) calloc((nb_var + 1), sizeof(double)); - - for (i = 1; i <= nb_var; i++) { - if (get_y(0, 1) == i) { - //CDEBUG0(surf_sdp_out,"-1 "); - a[i] = -1; - } else { - //CDEBUG0(surf_sdp_out,"0 "); - a[i] = 0; - } - } - - - /* - * Structure contraint blocks. - */ - block_num = 1; - matno = 1; - for (k = 1; k <= K; k++) { - for (i = 1; i <= pow(2, k - 1); i++) { - matno = get_y(k, 2 * i - 1); - CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num); - addentry(constraints, &C, matno, block_num, 1, 1, 1.0, - C.blocks[block_num].blocksize); - - matno = get_y(k, 2 * i); - CDEBUG2(surf_sdp_out, "%d %d 2 2 1", matno, block_num); - addentry(constraints, &C, matno, block_num, 2, 2, 1.0, - C.blocks[block_num].blocksize); - - matno = get_y(k - 1, i); - CDEBUG2(surf_sdp_out, "%d %d 1 2 1", matno, block_num); - addentry(constraints, &C, matno, block_num, 1, 2, 1.0, - C.blocks[block_num].blocksize); - - matno = get_y(k - 1, i); - CDEBUG2(surf_sdp_out, "%d %d 2 1 1", matno, block_num); - addentry(constraints, &C, matno, block_num, 2, 1, 1.0, - C.blocks[block_num].blocksize); - - isdiag[block_num] = 0; - block_num++; - } - } - - - /* - * Capacity constraint block. - */ - xbt_swag_foreach(cnst, cnst_list) { - - CDEBUG2(surf_sdp_out, "0 %d 1 1 %f", block_num, -(cnst->bound)); - addentry(constraints, &C, 0, block_num, 1, 1, -(cnst->bound), - C.blocks[block_num].blocksize); - - elem_list = &(cnst->element_set); - xbt_swag_foreach(elem, elem_list) { - if (elem->variable->weight <= 0) - break; - matno = get_y(K, elem->variable->index); - CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num, -(elem->value)); - addentry(constraints, &C, matno, block_num, 1, 1, -(elem->value), - C.blocks[block_num].blocksize); - - } - block_num++; - } - - - /* - * Positivy constraint blocks. - */ - for (i = 1; i <= pow(2, K); i++) { - matno = get_y(K, i); - CDEBUG2(surf_sdp_out, "%d %d 1 1 1", matno, block_num); - addentry(constraints, &C, matno, block_num, 1, 1, 1.0, - C.blocks[block_num].blocksize); - block_num++; - } - /* - * Latency constraint blocks. - */ - xbt_swag_foreach(var, var_list) { - var->value = 0.0; - if (var->weight && var->bound > 0) { - matno = get_y(K, var->index); - CDEBUG3(surf_sdp_out, "%d %d 1 1 %f", matno, block_num, var->bound); - addentry(constraints, &C, matno, block_num, 1, 1, var->bound, - C.blocks[block_num].blocksize); - } - } - - /* - * At this point, we'll stop to recognize whether any of the blocks - * are "hidden LP blocks" and correct the block type if needed. - */ - for (i = 1; i <= nb_cnsts; i++) { - if ((C.blocks[i].blockcategory != DIAG) && - (isdiag[i] == 1) && (C.blocks[i].blocksize > 1)) { - /* - * We have a hidden diagonal block! - */ - - blocksz = C.blocks[i].blocksize; - tempdiag = (double *) calloc((blocksz + 1), sizeof(double)); - for (j = 1; j <= blocksz; j++) - tempdiag[j] = C.blocks[i].data.mat[ijtok(j, j, blocksz)]; - free(C.blocks[i].data.mat); - C.blocks[i].data.vec = tempdiag; - C.blocks[i].blockcategory = DIAG; - }; - }; - - - /* - * Next, setup issparse and NULL out all nextbyblock pointers. - */ - p = NULL; - for (i = 1; i <= k; i++) { - p = constraints[i].blocks; - while (p != NULL) { - /* - * First, set issparse. - */ - if (((p->numentries) > 0.25 * (p->blocksize)) - && ((p->numentries) > 15)) { - p->issparse = 0; - } else { - p->issparse = 1; - }; - - if (C.blocks[p->blocknum].blockcategory == DIAG) - p->issparse = 1; - - /* - * Setup the cross links. - */ - - p->nextbyblock = NULL; - p = p->next; - }; - }; - - - /* - * Create cross link reference. - */ - create_cross_link(constraints, nb_var); - - - /* - * Debuging print problem in SDPA format. - */ - if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) { - DEBUG0("Printing SDPA..."); - tmp = strdup("SURF-PROPORTIONNAL.sdpa"); - write_prob(tmp, total_block_size, nb_var, C, a, constraints); - } - - /* - * Initialize parameters. - */ - DEBUG0("Initializing solution..."); - initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z); - - - - /* - * Call the solver. - */ - DEBUG0("Calling the solver..."); - stdout_sav = stdout; - stdout = fopen("/dev/null", "w"); - ret = - easy_sdp(total_block_size, nb_var, C, a, constraints, 0.0, &X, &y, - &Z, &pobj, &dobj); - fclose(stdout); - stdout = stdout_sav; - - switch (ret) { - case 0: - case 1: - DEBUG0("SUCCESS The problem is primal infeasible"); - break; - - case 2: - DEBUG0("SUCCESS The problem is dual infeasible"); - break; - - 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."); - break; - - case 4: - DEBUG0("Failure. Maximum number of iterations reached."); - break; - - case 5: - DEBUG0("Failure. Stuck at edge of primal feasibility."); - break; - - } - - if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) { - tmp = strdup("SURF-PROPORTIONNAL.sol"); - write_sol(tmp, total_block_size, nb_var, X, y, Z); - } - - /* - * Write out the solution if necessary. - */ - xbt_swag_foreach(cnst, cnst_list) { - - elem_list = &(cnst->element_set); - xbt_swag_foreach(elem, elem_list) { - if (elem->variable->weight <= 0) - break; - - i = (int) get_y(K, elem->variable->index); - elem->variable->value = y[i]; - - } - } - - - /* - * Free up memory. - */ - free_prob(total_block_size, nb_var, C, a, constraints, X, y, Z); - - free(isdiag); - free(tempdiag); - free(tmp); - - sys->modified = 0; - - if (XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) { - lmm_print(sys); - } - -} - - -/* - * Create the cross_link reference in order to have a byblock list. - */ -void create_cross_link(struct constraintmatrix *myconstraints, int k) -{ - - int i, j; - int blk; - struct sparseblock *p; - struct sparseblock *q; - - struct sparseblock *prev; - - /* - * Now, cross link. - */ - prev = NULL; - for (i = 1; i <= k; i++) { - p = myconstraints[i].blocks; - while (p != NULL) { - if (p->nextbyblock == NULL) { - blk = p->blocknum; - - /* - * link in the remaining blocks. - */ - for (j = i + 1; j <= k; j++) { - q = myconstraints[j].blocks; - - while (q != NULL) { - if (q->blocknum == p->blocknum) { - if (p->nextbyblock == NULL) { - p->nextbyblock = q; - q->nextbyblock = NULL; - prev = q; - } else { - prev->nextbyblock = q; - q->nextbyblock = NULL; - prev = q; - } - break; - } - q = q->next; - } - } - } - p = p->next; - } - } -} - - - - -void addentry(struct constraintmatrix *constraints, - struct blockmatrix *C, - int matno, - int blkno, int indexi, int indexj, double ent, int blocksize) -{ - struct sparseblock *p; - struct sparseblock *p_sav; - - p = constraints[matno].blocks; - - if (matno != 0.0) { - if (p == NULL) { - /* - * We haven't yet allocated any blocks. - */ - p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock)); - - //two entries because this library ignores indices starting in zerox - p->constraintnum = matno; - p->blocknum = blkno; - p->numentries = 1; - p->next = NULL; - - p->entries = calloc(p->numentries + 1, sizeof(double)); - p->iindices = calloc(p->numentries + 1, sizeof(int)); - p->jindices = calloc(p->numentries + 1, sizeof(int)); - - p->entries[p->numentries] = ent; - p->iindices[p->numentries] = indexi; - p->jindices[p->numentries] = indexj; - - p->blocksize = blocksize; - - constraints[matno].blocks = p; - } else { - /* - * We have some existing blocks. See whether this block is already - * in the chain. - */ - while (p != NULL) { - if (p->blocknum == blkno) { - /* - * Found the right block. - */ - p->constraintnum = matno; - p->blocknum = blkno; - p->numentries = p->numentries + 1; - - p->entries = - realloc(p->entries, (p->numentries + 1) * sizeof(double)); - p->iindices = - realloc(p->iindices, (p->numentries + 1) * sizeof(int)); - p->jindices = - realloc(p->jindices, (p->numentries + 1) * sizeof(int)); - - p->entries[p->numentries] = ent; - p->iindices[p->numentries] = indexi; - p->jindices[p->numentries] = indexj; - - return; - } - p_sav = p; - p = p->next; - } - - /* - * If we get here, we have a non-empty structure but not the right block - * inside hence create a new structure. - */ - - p = (struct sparseblock *) calloc(1, sizeof(struct sparseblock)); - - //two entries because this library ignores indices starting in zerox - p->constraintnum = matno; - p->blocknum = blkno; - p->numentries = 1; - p->next = NULL; - - p->entries = calloc(p->numentries + 1, sizeof(double)); - p->iindices = calloc(p->numentries + 1, sizeof(int)); - p->jindices = calloc(p->numentries + 1, sizeof(int)); - - p->entries[p->numentries] = ent; - p->iindices[p->numentries] = indexi; - p->jindices[p->numentries] = indexj; - - p->blocksize = blocksize; - - p_sav->next = p; - } - } else { - if (ent != 0.0) { - int blksz = C->blocks[blkno].blocksize; - if (C->blocks[blkno].blockcategory == DIAG) { - C->blocks[blkno].data.vec[indexi] = ent; - } else { - C->blocks[blkno].data.mat[ijtok(indexi, indexj, blksz)] = ent; - C->blocks[blkno].data.mat[ijtok(indexj, indexi, blksz)] = ent; - }; - }; - - } -}