+ /*
+ * 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) i++;
+ }
+ flows=i;
+ DEBUG1("Variable set : %d", xbt_swag_size(var_list));
+ 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));
+
+ /*
+ * 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.
+ */
+ double tmp_k;
+
+ tmp_k = (double) log((double)flows)/log(2.0);
+ K = (int) ceil(tmp_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);
+ CDEBUG1(surf_sdp_out,"%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);
+
+ /*
+ * The total number of constraints.
+ */
+ nb_cnsts = nb_cnsts_capacity + nb_cnsts_struct + nb_cnsts_positivy;
+ //fprintf(sdpout,"%d\n", 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;
+ //fprintf(sdpout,"2 ");
+ }else{
+ total_block_size += block_size = 1;
+ //fprintf(sdpout,"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++;
+ }
+
+ //fprintf(sdpout,"\n");
+
+
+ /*
+ * 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;
+ }
+ }
+ //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);
+
+ isdiag[block_num] = 0;
+ block_num++;
+ }
+ }
+
+
+ /*
+ * Capacity constraint block.
+ */
+ xbt_swag_foreach(cnst, cnst_list) {
+
+ //fprintf(sdpout,"0 %d 1 1 %f\n", 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);
+ //fprintf(sdpout,"%d %d 1 1 %f\n", elem->variable->index, 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);
+ //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);
+ block_num++;
+ }
+
+
+ /*
+ * 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!
+ */
+
+ //printf("Block %d is actually diagonal.\n",i);
+
+ 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.
+ */
+ struct sparseblock *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...\n");
+ char *tmp=strdup("SURF-PROPORTIONNAL.sdpa");
+ write_prob(tmp,total_block_size,nb_var,C,a,constraints);
+ free(tmp);
+ }
+
+ /*
+ * Initialize parameters.
+ */
+ DEBUG0("Initializing solution...\n");
+ initsoln(total_block_size, nb_var, C, a, constraints, &X, &y, &Z);
+
+
+
+ /*
+ * Call the solver.
+ */
+ DEBUG0("Calling the solver...\n");
+ FILE *stdout_sav=stdout;
+ stdout=fopen("/dev/null","w");
+ int 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\n");
+ break;
+
+ case 2: DEBUG0("SUCCESS The problem is dual infeasible\n");
+ 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.\n");
+ 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)) {
+ char *tmp=strdup("SURF-PROPORTIONNAL.sol");
+ write_sol(tmp,total_block_size, nb_var, X, y, Z);
+ free(tmp);
+ }
+
+ /*
+ * 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);
+
+ // fclose(sdpout);
+ free(isdiag);
+ sys->modified = 0;
+
+ if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
+ lmm_print(sys);
+ }
+
