lpkit.c

matlab整数规划工具箱

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 == 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 == 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 == 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;
  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) {
	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);
}   

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);
  if (lp->columns_scaled) {
    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] /
	  lp->scale[lp->rows + i];
  }
  else {
    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] == GE)
	printf(">= ");
      else if(lp->lag_con_type[i] == LE)
	printf("<= ");
      else if(lp->lag_con_type[i] == 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))
    return((REAL)1);

  /* scale = 1 / pow(10, (log10(min) + log10(max)) / 2); */
  /* Jon van Reet noticed: can be simplified to: */
  scale = 1 / sqrt(min * max);

  return(scale);
}

void unscale_columns(lprec *lp)
{
  int i, j;

  /* unscale mat */
  for(j = 1; j <= lp->columns; j++)
    for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
      lp->mat[i].value /= lp->scale[lp->rows + j];

  /* unscale bounds as well */
  for(i = lp->rows + 1; i <= lp->sum; i++) { /* was < */ /* changed by PN */
    if(lp->orig_lowbo[i] != 0)
      lp->orig_lowbo[i] *= lp->scale[i];
    if(lp->orig_upbo[i] != lp->infinite)
      lp->orig_upbo[i] *= lp->scale[i];
  }
    
  for(i = lp->rows + 1; i<= lp->sum; i++)
    lp->scale[i] = 1;
  lp->columns_scaled = FALSE;
  lp->eta_valid = FALSE;
}

void unscale(lprec *lp)
{
  int i, j;
  
  if(lp->scaling_used) {
    /* unscale mat */
    for(j = 1; j <= lp->columns; j++)
      for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	lp->mat[i].value /= lp->scale[lp->rows + j];

    /* unscale bounds */
    for(i = lp->rows + 1; i <= lp->sum; i++) { /* was < */ /* changed by PN */
      if(lp->orig_lowbo[i] != 0)
	lp->orig_lowbo[i] *= lp->scale[i];
      if(lp->orig_upbo[i] != lp->infinite)
	lp->orig_upbo[i] *= lp->scale[i];
    }
    
    /* unscale the matrix */
    for(j = 1; j <= lp->columns; j++)
      for(i = lp->col_end[j-1]; i < lp->col_end[j]; i++)
	lp->mat[i].value /= lp->scale[lp->mat[i].row_nr];

    /* unscale the rhs! */
    for(i = 0; i <= lp->rows; i++)
      lp->orig_rh[i] /= lp->scale[i];

    /* and don't forget to unscale the upper and lower bounds ... */
    for(i = 0; i <= lp->rows; i++) {
      if(lp->orig_lowbo[i] != 0)
	lp->orig_lowbo[i] /= lp->scale[i];
      if(lp->orig_upbo[i] != lp->infinite)
	lp->orig_upbo[i] /= lp->scale[i];
    }

    free(lp->scale);
    lp->scaling_used = FALSE;
    lp->eta_valid = FALSE;
  }
}


void auto_scale(lprec *lp)
{
  int i, j, row_nr;
  REAL *row_max, *row_min, *scalechange, absval;
  REAL col_max, col_min;
  
  if(!lp->scaling_used) {
    MALLOC(lp->scale, lp->sum_alloc + 1);
    for(i = 0; i <= lp->sum; i++)
      lp->scale[i] = 1;
  }
  
  MALLOC(row_max, lp->rows + 1);
  MALLOC(row_min, lp->rows + 1);