#include "xbt/log.h"
#include "xbt/sysdep.h"
-#include "xbt/mallocator.h"
#include "maxmin_private.h"
#include <stdlib.h>
/*
* SDP specific variables.
*/
+#define NOSHORTS
#include <declarations.h>
-static void create_cross_link(struct constraintmatrix *myconstraints, int k);
+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);
+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)");
+ "Logging specific to SURF (sdp)");
/*
########################################################################
######################## Simple Proportionnal fairness #################
# (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
# 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;
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 total_block_size=0;
- int *isdiag=NULL;
+ 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;
-
+ 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.
*/
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;
+ char *tmp = NULL;
FILE *stdout_sav;
int ret;
- if ( !(sys->modified))
+
+ if (!(sys->modified))
return;
/*
* Associate an index in the swag variables.
*/
var_list = &(sys->variable_set);
- i=0;
+ i = 0;
xbt_swag_foreach(var, var_list) {
- if(var->weight) i++;
+ if (var->weight != 0.0)
+ i++;
}
- flows=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--;
+ if (var->weight) {
+ var->index = i--;
+ }
}
- cnst_list=&(sys->active_constraint_set);
+ 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.
/*
* This number is found based on the tree structure explained on top.
*/
- tmp_k = (double) 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 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;
- CDEBUG1(surf_sdp_out,"Number of constraints = %d", nb_cnsts);
+ 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){
+ if (i <= nb_cnsts_struct) {
total_block_size += block_size = 2;
DEBUG0("2 ");
- }else{
+ } else {
total_block_size += block_size = 1;
- CDEBUG0(surf_sdp_out,"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++;
}
- CDEBUG0(surf_sdp_out," ");
+ 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;
+ 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;
}
}
- CDEBUG0(surf_sdp_out,"\n");
/*
*/
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\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);
- CDEBUG2(surf_sdp_out,"%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);
- CDEBUG2(surf_sdp_out,"%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);
- CDEBUG2(surf_sdp_out,"%d %d 2 1 1\n", matno , block_num);
- addentry(constraints, &C, matno, block_num, 2, 1, 1.0, C.blocks[block_num].blocksize);
-
+ 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\n", block_num, - (cnst->bound));
- addentry(constraints, &C, 0, block_num, 1, 1, - (cnst->bound) , C.blocks[block_num].blocksize);
+ 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\n", elem->variable->index, block_num, - (elem->value));
- addentry(constraints, &C, matno, block_num, 1, 1, - (elem->value), C.blocks[block_num].blocksize);
-
+ 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\n", matno, block_num);
- addentry(constraints, &C, matno, block_num, 1, 1, 1.0, C.blocks[block_num].blocksize);
+ 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)){
+ 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)];
+
+ 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;
+ 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;
+ 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)) {
+ 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);
- free(tmp);
+ 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);
-
+ 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);
+ 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;
+ stdout = stdout_sav;
- switch(ret){
+ switch (ret) {
case 0:
- case 1: DEBUG0("SUCCESS The problem is primal infeasible");
- break;
+ case 1:
+ DEBUG0("SUCCESS The problem is primal infeasible");
+ break;
+
+ case 2:
+ DEBUG0("SUCCESS The problem is dual 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 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 4:
+ DEBUG0("Failure. Maximum number of iterations reached.");
+ break;
- case 5: DEBUG0("Failure. Stuck at edge of primal feasibility.");
- 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);
- free(tmp);
+ 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);
}
/*
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];
-
+ 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);
- free(a);
- free(y);
-
- // fclose(sdpout);
- free(isdiag);
- sys->modified = 0;
-
- if(XBT_LOG_ISENABLED(surf_sdp, xbt_log_priority_debug)) {
- lmm_print(sys);
- }
-
+
+ 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){
+void create_cross_link(struct constraintmatrix *myconstraints, int k)
+{
int i, j;
int blk;
/*
* 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;
-
+ 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;
+ 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;
+ q = q->next;
}
}
}
- p=p->next;
+ p = p->next;
}
}
}
-
-void addentry(struct constraintmatrix *constraints,
- struct blockmatrix *C,
+
+void addentry(struct constraintmatrix *constraints,
+ struct blockmatrix *C,
int matno,
- int blkno,
- int indexi,
- int indexj,
- double ent,
- int blocksize)
+ int blkno, int indexi, int indexj, double ent, int blocksize)
{
struct sparseblock *p;
struct sparseblock *p_sav;
- p=constraints[matno].blocks;
-
+ p = constraints[matno].blocks;
+
if (matno != 0.0) {
- if (p == NULL){
+ if (p == NULL) {
/*
* We haven't yet allocated any blocks.
*/
- p=(struct sparseblock *)calloc(1, sizeof(struct sparseblock));
-
+ 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;
+ 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){
+ 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;
-
+ 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;
+ 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));
-
+
+ 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;
+ 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;
- };
+ 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;
+ };
};
-
+
}
}
-
