X-Git-Url: http://info.iut-bm.univ-fcomte.fr/pub/gitweb/simgrid.git/blobdiff_plain/065ecafc1ff3e06dc643403d2499e33ee9cb2ea7..dff9e15c44ab6340d27215957c56fa72fad246a2:/src/surf/sdp.c diff --git a/src/surf/sdp.c b/src/surf/sdp.c index 227d24fb4d..ed543c6487 100644 --- a/src/surf/sdp.c +++ b/src/surf/sdp.c @@ -6,24 +6,33 @@ * under the terms of the license (GNU LGPL) which comes with this package. */ -#include "xbt/sysdep.h" #include "xbt/log.h" +#include "xbt/sysdep.h" #include "maxmin_private.h" -#include - -#include +#include +#include #ifndef MATH #include #endif -XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf, - "Logging specific to SURF (sdp)"); +/* + * 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 ################# @@ -35,6 +44,10 @@ XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf, # (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 @@ -81,19 +94,23 @@ XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf, # 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 */ +// 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; @@ -101,177 +118,590 @@ XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_sdp, surf, void sdp_solve(lmm_system_t sys) { - lmm_variable_t var = NULL; /* * SDP mapping variables. */ - int K=0; + 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 block_num=0; - int block_size; - int *isdiag=NULL; - FILE *sdpout = fopen("SDPA-printf.tmp","w"); - int matno=0; - int i=0; - int j=0; - int k=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; + struct blockmatrix X, Z; double *y; double pobj, dobj; double *a; - struct constraintmatrix *constraints; - + 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; - double min_usage = -1; - */ + lmm_variable_t var = NULL; + + double tmp_k; + struct sparseblock *p; + char *tmp = NULL; + FILE *stdout_sav; + int ret; + - if ( !(sys->modified)) + + 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. */ - flows = xbt_swag_size(&(sys->variable_set)); links = xbt_swag_size(&(sys->active_constraint_set)); /* * This number is found based on the tree structure explained on top. */ - K = (int)log((double)flows)/log(2.0); - + 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 style. + * The number of variables in the SDP program. */ - nb_var = get_y(K, pow(2,K)); - + 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_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 = + 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; + 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; + 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)); + + a = (double *) calloc(nb_var + 1, sizeof(double)); /* * Block sizes. */ - block_num=1; - block_size=0; + block_num = 1; + block_size = 0; /* * For each constraint block do. */ - for(i = 1; i <= nb_cnsts; i++){ - + for (i = 1; i <= nb_cnsts; i++) { + /* * Structured blocks are size 2 and all others are size 1. */ - if(i <= nb_cnsts_struct){ - block_size = 2; - fprintf(sdpout,"2 "); - }else{ - block_size = 1; - fprintf(sdpout,"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)); - + C.blocks[block_num].data.mat = + (double *) calloc(block_size * block_size, sizeof(double)); + block_num++; } - fprintf(sdpout,"\n"); + 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){ - fprintf(sdpout,"-1 "); - a[i]=-1; - }else{ - fprintf(sdpout,"0 "); - a[i]=0; + 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; } } - fprintf(sdpout,"\n"); /* * 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); - fprintf(sdpout,"%d %d 1 1 1\n", matno , block_num); - addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize); - - matno=get_y(k,2*i); - fprintf(sdpout,"%d %d 2 2 1\n", matno , block_num); - addentry(constraints, &C, matno, block_num, 2, 2, 1.0, C.blocks[block_num].blocksize); - - matno=get_y(k-1,i); - fprintf(sdpout,"%d %d 1 2 1\n", matno , block_num); - addentry(constraints, &C, matno, block_num, 1, 2, 1.0, C.blocks[block_num].blocksize); - - matno=get_y(k-1,i); - fprintf(sdpout,"%d %d 2 1 1\n", matno , block_num); - addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize); - + 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; + }; + }; + + } }