lpkit.c

生成直角Steiner树的程序包

  if(lp->scaling_used)
    value /= lp->scale[lp->rows + column];
  if(value < lp->orig_lowbo[lp->rows + column])
    error("Upperbound must be >= lowerbound"); 
  lp->eta_valid = FALSE;
  lp->orig_upbo[lp->rows + column] = value;
}

void set_lowbo(lprec *lp, int column, REAL value)
{
  if(column > lp->columns || column < 1)
    error("Column out of range");
  if(lp->scaling_used)
    value /= lp->scale[lp->rows + column];
  if(value > lp->orig_upbo[lp->rows + column])
    error("Upperbound must be >= lowerbound"); 
  /*
    if(value < 0)
    error("Lower bound cannot be < 0");
    */
  lp->eta_valid = FALSE;
  lp->orig_lowbo[lp->rows + column] = value;
}

void set_bounds(lprec *lp, int column, REAL lower, REAL upper)
{
  if(column > lp->columns || column < 1)
    error("Column out of range");
  if(lp->scaling_used)
    {
      lower /= lp->scale[lp->rows + column];
      upper /= lp->scale[lp->rows + column];
    }
  if(lower > upper)
    error("Upperbound must be >= lowerbound");
  lp->eta_valid = FALSE;
  lp->orig_lowbo[lp->rows+column] = lower;
  lp->orig_upbo[lp->rows+column] = upper;
}

void set_int(lprec *lp, int column, short must_be_int)
{
  if(column > lp->columns || column < 1)
    error("Column out of range");
  lp->must_be_int[lp->rows + column] = must_be_int;
  if(must_be_int == TRUE)
    if(lp->columns_scaled)
      unscale_columns(lp);
}

void set_rh(lprec *lp, int row, REAL value)
{
  if(row > lp->rows || row < 0)
    error("Row out of Range");

  if ((row == 0) && (!lp->maximise))  /* setting of RHS of OF IS meaningful */
    value = -value;
  if(lp->scaling_used) {
    if(lp->ch_sign[row])
      lp->orig_rh[row] = -value * lp->scale[row];
    else
      lp->orig_rh[row] = value * lp->scale[row];
  }
  else
    if(lp->ch_sign[row])
      lp->orig_rh[row] = -value;
    else
      lp->orig_rh[row] = value;
  lp->eta_valid = FALSE;
} 

void set_rh_vec(lprec *lp, REAL *rh)
{
  int i;
  if(lp->scaling_used) {
    for(i = 1; i <= lp->rows; i++)
      if(lp->ch_sign[i])
        lp->orig_rh[i] = -rh[i]*lp->scale[i];
      else
        lp->orig_rh[i] = rh[i]*lp->scale[i];
  }
  else
    for(i = 1; i <= lp->rows; i++)
      if(lp->ch_sign[i])
        lp->orig_rh[i] = -rh[i];
      else
        lp->orig_rh[i] = rh[i];
  lp->eta_valid = FALSE;
}

void str_set_rh_vec(lprec *lp, char *rh_string)
{
  int  i;
  REAL *newrh;
  char *p, *newp;
  CALLOC(newrh, lp->rows + 1);
  p = rh_string;
 
  for(i = 1; i <= lp->rows; i++) {
    newrh[i] = (REAL) strtod(p, &newp);
    if(p == newp)
      error("Bad string in str_set_rh_vec");
    else
      p = newp; 
  }
  set_rh_vec(lp, newrh);
  free(newrh);
}


void set_maxim(lprec *lp)
{
  int i;
  if(lp->maximise == FALSE) {
    for(i = 0; i < lp->non_zeros; i++)
      if(lp->mat[i].row_nr == 0)
	lp->mat[i].value *= -1;
    lp->eta_valid = FALSE;
    lp->orig_rh[0] *= -1;
  } 
  lp->maximise = TRUE;
  lp->ch_sign[0] = TRUE;
}

void set_minim(lprec *lp)
{
  int i;
  if(lp->maximise == TRUE) {
    for(i = 0; i < lp->non_zeros; i++)
      if(lp->mat[i].row_nr == 0)
	lp->mat[i].value = -lp->mat[i].value;
    lp->eta_valid = FALSE;
    lp->orig_rh[0] *= -1;
  } 
  lp->maximise = FALSE;
  lp->ch_sign[0] = FALSE;
}

void set_constr_type(lprec *lp, int row, short con_type)
{
  int i;
  if(row > lp->rows || row < 1)
    error("Row out of Range");
  if(con_type == REL_EQ) {
    lp->orig_upbo[row] = 0;
    lp->basis_valid = FALSE;
    if(lp->ch_sign[row]) {
      for(i = 0; i < lp->non_zeros; i++)
	if(lp->mat[i].row_nr == row)
	  lp->mat[i].value *= -1;
      lp->eta_valid = FALSE;
      lp->ch_sign[row] = FALSE;
      if(lp->orig_rh[row] != 0)
	lp->orig_rh[row] *= -1;
    }
  }
  else if(con_type == REL_LE) {
    lp->orig_upbo[row] = lp->infinite;
    lp->basis_valid = FALSE;
    if(lp->ch_sign[row]) {
      for(i = 0; i < lp->non_zeros; i++)
	if(lp->mat[i].row_nr == row)
	  lp->mat[i].value *= -1;
      lp->eta_valid = FALSE;
      lp->ch_sign[row] = FALSE;
      if(lp->orig_rh[row] != 0)
	lp->orig_rh[row] *= -1;
    }
  }
  else if(con_type == REL_GE) {
    lp->orig_upbo[row] = lp->infinite;
    lp->basis_valid = FALSE;
    if(!lp->ch_sign[row]) {
      for(i = 0; i < lp->non_zeros; i++)
	if(lp->mat[i].row_nr == row)
	  lp->mat[i].value *= -1;
      lp->eta_valid = FALSE;
      lp->ch_sign[row] = TRUE;
      if(lp->orig_rh[row] != 0)
	lp->orig_rh[row] *= -1;
    }
  } 
  else
    error("Constraint type not (yet) implemented");
}

REAL mat_elm(lprec *lp, int row, int column)
{
  REAL value;
  int elmnr;
  if(row < 0 || row > lp->rows)
    error("Row out of range in mat_elm");
  if(column < 1 || column > lp->columns)
    error("Column out of range in mat_elm");
  value = 0;
  elmnr = lp->col_end[column-1];
  while(lp->mat[elmnr].row_nr != row && elmnr < lp->col_end[column])
    elmnr++;
  if(elmnr != lp->col_end[column]) {
    value = lp->mat[elmnr].value;
    if(lp->ch_sign[row])
      value = -value;
    if(lp->scaling_used)
      value /= lp->scale[row] * lp->scale[lp->rows + column];
  }
  return(value);
}


void get_row(lprec *lp, int row_nr, REAL *row)
{
  int i, j;

  if(row_nr <0 || row_nr > lp->rows)
    error("Row nr. out of range in get_row");
  for(i = 1; i <= lp->columns; i++) {
    row[i] = 0;
    for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
      if(lp->mat[j].row_nr == row_nr)
	row[i] = lp->mat[j].value;
    if(lp->scaling_used)
      row[i] /= lp->scale[lp->rows + i] * lp->scale[row_nr];
  }
  if(lp->ch_sign[row_nr])
    for(i = 0; i <= lp->columns; i++)
      if(row[i] != 0)
        row[i] = -row[i];
}

void get_column(lprec *lp, int col_nr, REAL *column)
{
  int i;

  if(col_nr < 1 || col_nr > lp->columns)
    error("Col. nr. out of range in get_column");
  for(i = 0; i <= lp->rows; i++)
    column[i] = 0;
  for(i = lp->col_end[col_nr - 1]; i < lp->col_end[col_nr]; i++)
    column[lp->mat[i].row_nr] = lp->mat[i].value;
  for(i = 0; i <= lp->rows; i++)
    if(column[i] != 0) {
      if(lp->ch_sign[i])
	column[i] *= -1;
      if(lp->scaling_used)
	column[i] /= (lp->scale[i] * lp->scale[lp->rows + col_nr]);
    }
}

void get_reduced_costs(lprec *lp, REAL *rc)
{
  int varnr, i, j, row;
  REAL f;

  if(!lp->basis_valid)
    error("Not a valid basis in get_reduced_costs");

  if(!lp->eta_valid)
    invert(lp);  

  for(i = 1; i <= lp->sum; i++)
    rc[i] = 0;
  rc[0] = 1;

  btran(lp, rc);

  for(i = 1; i <= lp->columns; i++) {
    varnr = lp->rows + i;
    if(!lp->basis[varnr])
      if(lp->upbo[varnr] > 0) {
	if(lp->scaling_used) {
	  f = 0;
	  for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
	    row = lp->mat[j].row_nr;
	    f += rc[row] * lp->mat[j].value
		  * (lp->scale[row] * lp->scale[varnr]);
	  }
	  rc[varnr] = f;
	}
	else {
	  f = 0;
	  for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	    f += rc[lp->mat[j].row_nr] * lp->mat[j].value;
	  rc[varnr] = f;
	}
      }
  }
  for(i = 1; i <= lp->sum; i++)
    my_round(rc[i], lp->epsd);
}

void get_slack_vars(lprec *lp, REAL *slacks)
{
  int    i, basi;

  if(!lp->basis_valid)
    error("Not a valid basis in get_slack_vars");

  /* zero all results of rows */
  memset(slacks, '\0', (lp->rows + 1) * sizeof(REAL));

  /* Lower bounds are always zero for slack variables.  Their upper	*/
  /* bounds are always infinity (except for equality constraints).	*/
  /* However they NEVER appear out of the basis except at their lower	*/
  /* bounds.								*/

  if(lp->scaling_used) {
    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi <= lp->rows)
	slacks[basi] = lp->rhs[i] * lp->scale[basi];
    }
  }
  else {
    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi <= lp->rows)
	slacks[basi] = lp->rhs[i];
    }
  }
}

short is_feasible(lprec *lp, REAL *values)
{
  int i, elmnr;
  REAL *this_rhs;
  REAL dist;

  if(lp->scaling_used) {
    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(   values[i - lp->rows] < lp->orig_lowbo[i] * lp->scale[i]
	 || values[i - lp->rows] > lp->orig_upbo[i]  * lp->scale[i])
	return(FALSE);
  }
  else {
    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(   values[i - lp->rows] < lp->orig_lowbo[i]
	 || values[i - lp->rows] > lp->orig_upbo[i])
	return(FALSE);
  }
  CALLOC(this_rhs, lp->rows + 1);
    for(i = 1; i <= lp->columns; i++)
      for(elmnr = lp->col_end[i - 1]; elmnr < lp->col_end[i]; elmnr++)
	this_rhs[lp->mat[elmnr].row_nr] += lp->mat[elmnr].value * values[i]; 
  for(i = 1; i <= lp->rows; i++) {
    dist = lp->orig_rh[i] - this_rhs[i];
    my_round(dist, 0.001) /* ugly constant, MB */
      if((lp->orig_upbo[i] == 0 && dist != 0) || dist < 0) {
	free(this_rhs);
	return(FALSE);
      }     
  } 
  free(this_rhs);
  return(TRUE);
}

/* fixed by Enrico Faggiolo */
short column_in_lp(lprec *lp, REAL *testcolumn)
{
  int i, j;
  int nz, ident;
  REAL value;

  for(nz = 0, i = 0; i <= lp->rows; i++)
    if(my_abs(testcolumn[i]) > lp->epsel) nz++;

  if(lp->scaling_used)
    for(i = 1; i <= lp->columns; i++) {
      ident = nz;
      for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
	value = lp->mat[j].value;
	if(lp->ch_sign[lp->mat[j].row_nr])
	  value = -value;
	value /= lp->scale[lp->rows + i];
	value /= lp->scale[lp->mat[j].row_nr];
	value -= testcolumn[lp->mat[j].row_nr];
	if(my_abs(value) > lp->epsel)
	  break;
	ident--;
	if(ident == 0)
	  return(TRUE);
      }
    }
  else
    for(i = 1; i <= lp->columns; i++) {
      ident = nz;
      for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
	value = lp->mat[j].value;
	if(lp->ch_sign[lp->mat[j].row_nr])
	  value = -value;
	value -= testcolumn[lp->mat[j].row_nr];
	if(my_abs(value) > lp->epsel)
	  break;
	ident--;
	if(ident == 0)
	  return(TRUE);
      }
    }
  return(FALSE);
}

void print_lp(lprec *lp)
{
  int i, j;
  REAL *fatmat;
  CALLOC(fatmat, (lp->rows + 1) * lp->columns);
  for(i = 1; i <= lp->columns; i++)
    for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
      fatmat[(i - 1) * (lp->rows + 1) + lp->mat[j].row_nr] = lp->mat[j].value;

  printf("problem name: %s\n", lp->lp_name);
  printf("          ");
  for(j = 1; j <= lp->columns; j++)
    if(lp->names_used)
      printf("%8s ", lp->col_name[j]);
    else
      printf("Var[%3d] ", j);
  if(lp->maximise)
    {
      printf("\nMaximise  ");
      for(j = 0; j < lp->columns; j++)
	printf("%8g ", (double) -fatmat[j * (lp->rows + 1)]);
    }
  else
    {
      printf("\nMinimize  ");
      for(j = 0; j < lp->columns; j++)
	printf("%8g ", (double) fatmat[j * (lp->rows + 1)]);
    }
  printf("\n");
  for(i = 1; i <= lp->rows; i++) {
    if(lp->names_used)
      printf("%-9s ", lp->row_name[i]);
    else
      printf("Row[%3d]  ", i);
    for(j = 0; j < lp->columns; j++)
      if(lp->ch_sign[i] && fatmat[j*(lp->rows + 1) + i] != 0)
	printf("%8g ", (double)-fatmat[j * (lp->rows+1) + i]);
      else
	printf("%8g ", (double)fatmat[j * (lp->rows + 1) + i]);
    if(lp->orig_upbo[i] != 0) {
      if(lp->ch_sign[i])
	printf(">= ");
      else
	printf("<= ");
    }
    else
      printf(" = ");
    if(lp->ch_sign[i])
      printf("%8g", (double)-lp->orig_rh[i]);
    else
      printf("%8g", (double)lp->orig_rh[i]);
    if(lp->orig_lowbo[i] != 0) {
      printf("  %s = %8g", (lp->ch_sign[i]) ? "lowbo" : "upbo",
	     (double)lp->orig_lowbo[i]);
    }
    if((lp->orig_upbo[i] != lp->infinite) && (lp->orig_upbo[i] != 0.0)) {
      printf("  %s = %8g", (lp->ch_sign[i]) ? "upbo" : "lowbo",
	     (double)lp->orig_upbo[i]);
    }
    printf("\n");
  }
  printf("Type      ");
  for(i = 1; i <= lp->columns; i++)
    if(lp->must_be_int[lp->rows + i] == TRUE)
      printf("     Int ");
    else
      printf("    Real ");
  printf("\nupbo      ");
  for(i = 1; i <= lp->columns; i++)
    if(lp->orig_upbo[lp->rows + i] == lp->infinite)
      printf(" Infinite");
    else
      printf("%8g ", (double)lp->orig_upbo[lp->rows + i]);
  printf("\nlowbo     ");
  for(i = 1; i <= lp->columns; i++)
    printf("%8g ", (double)lp->orig_lowbo[lp->rows + i]);
  printf("\n");
  for(i = 0; i < lp->nr_lagrange; i++) {
    printf("lag[%d]  ", i);
    for(j = 1; j <= lp->columns; j++)
      printf("%8g ", (double)lp->lag_row[i][j]);
    if(lp->orig_upbo[i] == lp->infinite) {
      if(lp->lag_con_type[i] == REL_GE)
	printf(">= ");
      else if(lp->lag_con_type[i] == REL_LE)
	printf("<= ");
      else if(lp->lag_con_type[i] == REL_EQ)
	printf(" = ");
    }
    printf("%8g\n", (double)lp->lag_rhs[i]);
  }

  free(fatmat);
}  

void set_row_name(lprec *lp, int row, nstring new_name)
{
  int i;
  hashelem *hp;

  if(!lp->names_used) {
    CALLOC(lp->row_name, lp->rows_alloc + 1);
    CALLOC(lp->col_name, lp->columns_alloc + 1);
    lp->names_used = TRUE;
    for(i = 0; i <= lp->rows; i++)
      sprintf(lp->row_name[i], "r_%d", i);
    for(i = 1; i <= lp->columns; i++)
      sprintf(lp->col_name[i], "var_%d", i);
  }
  strcpy(lp->row_name[row], new_name);
  hp = puthash(lp->row_name[row], lp->rowname_hashtab);
  hp->index = row;
}

void set_col_name(lprec *lp, int column, nstring new_name)
{
  int i;
  hashelem *hp;
 
  if(!lp->names_used) {
    CALLOC(lp->row_name, lp->rows_alloc + 1);
    CALLOC(lp->col_name, lp->columns_alloc + 1);
    lp->names_used = TRUE;
    for(i = 0; i <= lp->rows; i++)
      sprintf(lp->row_name[i], "r_%d", i);
    for(i = 1; i <= lp->columns; i++)
      sprintf(lp->col_name[i], "var_%d", i);
  }
  strcpy(lp->col_name[column], new_name);
  hp = puthash(lp->col_name[column], lp->colname_hashtab);
  hp->index = column;
}

static REAL minmax_to_scale(REAL min, REAL max)
{
  REAL scale;

  /* should do something sensible when min or max is 0, MB */
  if((min == 0) || (max == 0))