lpkit.c

工具箱 (整数线性规划工具箱) Matlab中使用

#include "lpkit.h"
#include "lpglob.h"
#include <stdlib.h>
#include <stdarg.h>
#include <stdio.h>
#include <string.h>

#define HASHSIZE 10007

/* Globals */
int     Rows;
int     Columns;
int     Sum;
int     Non_zeros;
int     Level;

REAL	Trej;

short   Maximise;
REAL    Extrad;

int     Warn_count; /* used in CHECK version of rounding macro */

void error(char *format, ...)
{
  va_list ap;

  va_start(ap, format);
  vfprintf(stderr, format, ap);
  fputc('\n', stderr);
  va_end(ap);

  exit(EXIT_FAILURE);
}

lprec *make_lp(int rows, int columns)
{
  lprec *newlp;
  int i, sum;  

  if(rows < 0 || columns < 0)
    error("rows < 0 or columns < 0");

  sum = rows + columns;

  CALLOC(newlp, 1);

  strcpy(newlp->lp_name, "unnamed");

  newlp->verbose = FALSE;
  newlp->print_duals = FALSE;
  newlp->print_sol = FALSE;
  newlp->debug = FALSE;
  newlp->print_at_invert = FALSE;
  newlp->trace = FALSE;

  newlp->rows = rows;
  newlp->columns = columns;
  newlp->sum = sum;
  newlp->rows_alloc = rows;
  newlp->columns_alloc = columns;
  newlp->sum_alloc = sum;
  newlp->names_used = FALSE;

  newlp->obj_bound = DEF_INFINITE;
  newlp->infinite = DEF_INFINITE;
  newlp->epsilon = DEF_EPSILON;
  newlp->epsb = DEF_EPSB;
  newlp->epsd = DEF_EPSD;
  newlp->epsel = DEF_EPSEL;
  newlp->non_zeros = 0;
  newlp->mat_alloc = 1;
  CALLOC(newlp->mat, newlp->mat_alloc);
  CALLOC(newlp->col_no, newlp->mat_alloc + 1);
  CALLOC(newlp->col_end, columns + 1);
  CALLOC(newlp->row_end, rows + 1);
  newlp->row_end_valid = FALSE;
  CALLOC(newlp->orig_rh, rows + 1);
  CALLOC(newlp->rh, rows + 1);
  CALLOC(newlp->rhs, rows + 1);

  CALLOC(newlp->must_be_int, sum + 1);
  for(i = 0; i <= sum; i++)
    newlp->must_be_int[i] = FALSE;

  CALLOC(newlp->orig_upbo, sum + 1);
  for(i = 0; i <= sum; i++)
    newlp->orig_upbo[i] = newlp->infinite;

  CALLOC(newlp->upbo, sum + 1);
  CALLOC(newlp->orig_lowbo, sum + 1);
  CALLOC(newlp->lowbo, sum + 1);

  newlp->basis_valid = TRUE;
  CALLOC(newlp->bas, rows+1);
  CALLOC(newlp->basis, sum + 1);
  CALLOC(newlp->lower, sum + 1);

  for(i = 0; i <= rows; i++) {
    newlp->bas[i] = i;
    newlp->basis[i] = TRUE;
  }

  for(i = rows + 1; i <= sum; i++)
    newlp->basis[i] = FALSE;


  for(i = 0 ; i <= sum; i++)
    newlp->lower[i] = TRUE;
 
  newlp->eta_valid = TRUE;
  newlp->eta_size = 0;
  newlp->eta_alloc = INITIAL_MAT_SIZE;
  newlp->max_num_inv = DEFNUMINV;

  newlp->nr_lagrange = 0;

  CALLOC(newlp->eta_value, newlp->eta_alloc);
  CALLOC(newlp->eta_row_nr, newlp->eta_alloc);

  /* +1 reported by Christian Rank */
  CALLOC(newlp->eta_col_end, newlp->rows_alloc + newlp->max_num_inv + 1);

  newlp->bb_rule = FIRST_NI;
  newlp->break_at_int = FALSE;
  newlp->break_value = 0;

  newlp->iter = 0;
  newlp->total_iter = 0;

  CALLOC(newlp->solution, sum + 1);
  CALLOC(newlp->best_solution, sum + 1);
  CALLOC(newlp->duals, rows + 1);

  newlp->maximise = FALSE;
  newlp->floor_first = TRUE;

  newlp->scaling_used = FALSE;
  newlp->columns_scaled = FALSE;

  CALLOC(newlp->ch_sign, rows + 1);

  for(i = 0; i <= rows; i++)
    newlp->ch_sign[i] = FALSE;

  newlp->valid = FALSE; 

  /* create two hash tables for names */
  newlp->rowname_hashtab = create_hash_table(HASHSIZE);
  newlp->colname_hashtab = create_hash_table(HASHSIZE);

  return(newlp);
}

void delete_lp(lprec *lp)
{
  int i; 

  if(lp->names_used) {
    free(lp->row_name);
    free(lp->col_name);
  }

  free(lp->mat);
  free(lp->col_no);
  free(lp->col_end);
  free(lp->row_end);
  free(lp->orig_rh);
  free(lp->rh);
  free(lp->rhs);
  free(lp->must_be_int);
  free(lp->orig_upbo);
  free(lp->orig_lowbo);
  free(lp->upbo);
  free(lp->lowbo);
  free(lp->bas);
  free(lp->basis);
  free(lp->lower);
  free(lp->eta_value);
  free(lp->eta_row_nr);
  free(lp->eta_col_end);
  free(lp->solution);
  free(lp->best_solution);
  free(lp->duals);
  free(lp->ch_sign);
  if(lp->scaling_used)
    free(lp->scale);
  if(lp->nr_lagrange > 0) {
    free(lp->lag_rhs);
    free(lp->lambda);
    free(lp->lag_con_type);
    for(i = 0; i < lp->nr_lagrange; i++)
      free(lp->lag_row[i]);
    free(lp->lag_row);
  }

  free_hash_table(lp->rowname_hashtab);
  free_hash_table(lp->colname_hashtab);

  free(lp);

}  

lprec *copy_lp(lprec *lp)
{
  lprec *newlp;
  int i, rowsplus, colsplus, sumplus;
  
  rowsplus = lp->rows_alloc + 1;
  colsplus = lp->columns_alloc + 1;
  sumplus  = lp->sum_alloc + 1;

  MALLOCCPY(newlp, lp, 1); /* copy all non pointers */

  if(newlp->names_used) {
    MALLOCCPY(newlp->col_name, lp->col_name, colsplus);
    MALLOCCPY(newlp->row_name, lp->row_name, rowsplus);
  }

  newlp->rowname_hashtab = copy_hash_table(lp->rowname_hashtab);
  newlp->colname_hashtab = copy_hash_table(lp->colname_hashtab);

  MALLOCCPY(newlp->mat, lp->mat, newlp->mat_alloc);
  MALLOCCPY(newlp->col_end, lp->col_end, colsplus);
  MALLOCCPY(newlp->col_no, lp->col_no, newlp->mat_alloc + 1);
  MALLOCCPY(newlp->row_end, lp->row_end, rowsplus);
  MALLOCCPY(newlp->orig_rh, lp->orig_rh, rowsplus);
  MALLOCCPY(newlp->rh, lp->rh, rowsplus);
  MALLOCCPY(newlp->rhs, lp->rhs, rowsplus);
  MALLOCCPY(newlp->must_be_int, lp->must_be_int, sumplus);
  MALLOCCPY(newlp->orig_upbo, lp->orig_upbo, sumplus);
  MALLOCCPY(newlp->orig_lowbo, lp->orig_lowbo, sumplus);
  MALLOCCPY(newlp->upbo, lp->upbo, sumplus);
  MALLOCCPY(newlp->lowbo, lp->lowbo, sumplus);
  MALLOCCPY(newlp->bas, lp->bas, rowsplus);
  MALLOCCPY(newlp->basis, lp->basis, sumplus);
  MALLOCCPY(newlp->lower, lp->lower, sumplus);
  MALLOCCPY(newlp->eta_value, lp->eta_value, lp->eta_alloc);
  MALLOCCPY(newlp->eta_row_nr, lp->eta_row_nr, lp->eta_alloc);
  MALLOCCPY(newlp->eta_col_end, lp->eta_col_end,
	    lp->rows_alloc + lp->max_num_inv + 1);
  MALLOCCPY(newlp->solution, lp->solution, sumplus);
  MALLOCCPY(newlp->best_solution, lp->best_solution, sumplus);
  MALLOCCPY(newlp->duals, lp->duals, rowsplus);
  MALLOCCPY(newlp->ch_sign, lp->ch_sign, rowsplus);

  if(newlp->scaling_used)
    MALLOCCPY(newlp->scale, lp->scale, sumplus);

  if(newlp->nr_lagrange > 0) {
    MALLOCCPY(newlp->lag_rhs, lp->lag_rhs, newlp->nr_lagrange);
    MALLOCCPY(newlp->lambda, lp->lambda, newlp->nr_lagrange);
    MALLOCCPY(newlp->lag_con_type, lp->lag_con_type, newlp->nr_lagrange);
    MALLOC(newlp->lag_row, newlp->nr_lagrange);
    for(i = 0; i < newlp->nr_lagrange; i++)
      MALLOCCPY(newlp->lag_row[i], lp->lag_row[i], colsplus);
  }
  return(newlp);
}

void inc_mat_space(lprec *lp, int maxextra)
{
   if(lp->non_zeros + maxextra >= lp->mat_alloc) {
     /* let's allocate at least INITIAL_MAT_SIZE  entries */
     if(lp->mat_alloc < INITIAL_MAT_SIZE)
       lp->mat_alloc = INITIAL_MAT_SIZE;
     
     /* increase the size by 50% each time it becomes too small */
     while(lp->non_zeros + maxextra >= lp->mat_alloc)
       lp->mat_alloc *= 1.5;

     REALLOC(lp->mat, lp->mat_alloc);
     REALLOC(lp->col_no, lp->mat_alloc + 1);
   }
}
 
void inc_row_space(lprec *lp)
{
  if(lp->rows > lp->rows_alloc) {
    lp->rows_alloc = lp->rows+10;
    lp->sum_alloc  = lp->rows_alloc + lp->columns_alloc;
    REALLOC(lp->orig_rh, lp->rows_alloc + 1);
    REALLOC(lp->rh, lp->rows_alloc + 1);
    REALLOC(lp->rhs, lp->rows_alloc + 1);
    REALLOC(lp->orig_upbo, lp->sum_alloc + 1);
    REALLOC(lp->upbo, lp->sum_alloc + 1);
    REALLOC(lp->orig_lowbo, lp->sum_alloc + 1);
    REALLOC(lp->lowbo, lp->sum_alloc + 1);
    REALLOC(lp->solution, lp->sum_alloc + 1);
    REALLOC(lp->best_solution, lp->sum_alloc + 1);
    REALLOC(lp->row_end, lp->rows_alloc + 1);
    REALLOC(lp->basis, lp->sum_alloc + 1);
    REALLOC(lp->lower, lp->sum_alloc + 1);
    REALLOC(lp->must_be_int, lp->sum_alloc + 1);
    REALLOC(lp->bas, lp->rows_alloc + 1);
    REALLOC(lp->duals, lp->rows_alloc + 1);
    REALLOC(lp->ch_sign, lp->rows_alloc + 1);
    REALLOC(lp->eta_col_end, lp->rows_alloc + lp->max_num_inv + 1);
    if(lp->names_used)
      REALLOC(lp->row_name, lp->rows_alloc + 1);
    if(lp->scaling_used)
      REALLOC(lp->scale, lp->sum_alloc + 1);
  }
}

void inc_col_space(lprec *lp)
{
  if(lp->columns >= lp->columns_alloc) {
    lp->columns_alloc = lp->columns + 10;
    lp->sum_alloc = lp->rows_alloc + lp->columns_alloc;
    REALLOC(lp->must_be_int, lp->sum_alloc + 1);
    REALLOC(lp->orig_upbo, lp->sum_alloc + 1);
    REALLOC(lp->upbo, lp->sum_alloc + 1);
    REALLOC(lp->orig_lowbo, lp->sum_alloc + 1);
    REALLOC(lp->lowbo, lp->sum_alloc + 1);
    REALLOC(lp->solution, lp->sum_alloc + 1);
    REALLOC(lp->best_solution, lp->sum_alloc + 1);
    REALLOC(lp->basis, lp->sum_alloc + 1);
    REALLOC(lp->lower, lp->sum_alloc + 1);
    if(lp->names_used)
      REALLOC(lp->col_name, lp->columns_alloc + 1);
    if(lp->scaling_used)
      REALLOC(lp->scale, lp->sum_alloc + 1);
    REALLOC(lp->col_end, lp->columns_alloc + 1);
  }
}

void set_mat(lprec *lp, int Row, int Column, REAL Value)
{
  int elmnr, lastelm, i;

  /* This function is very inefficient if used to add new matrix entries in
     other places than at the end of the matrix. OK for replacing existing
     non-zero values */

  if(Row > lp->rows || Row < 0)
    error("Row out of range");
  if(Column > lp->columns || Column < 1)
    error("Column out of range");

  /* scaling is performed twice? MB */
  if(lp->scaling_used)
    Value *= lp->scale[Row] * lp->scale[lp->rows + Column];
  
  if (lp->basis[Column] == TRUE && Row > 0)
    lp->basis_valid = FALSE;
  lp->eta_valid = FALSE;

  /* find out if we already have such an entry */
  elmnr = lp->col_end[Column - 1];
  while((elmnr < lp->col_end[Column]) && (lp->mat[elmnr].row_nr != Row))
    elmnr++;
  
  if((elmnr != lp->col_end[Column]) && (lp->mat[elmnr].row_nr == Row)) {
    /* there is an existing entry */
    if(my_abs(Value) > lp->epsilon) { /* we replace it by something non-zero */
      if (lp->scaling_used) {
	if(lp->ch_sign[Row])
	  lp->mat[elmnr].value = -Value * lp->scale[Row] * lp->scale[Column];
	else
	  lp->mat[elmnr].value = Value * lp->scale[Row] * lp->scale[Column];
      }
      else { /* no scaling */
	if(lp->ch_sign[Row])
	  lp->mat[elmnr].value = -Value;
	else
	  lp->mat[elmnr].value = Value;
      }
    }
    else { /* setting existing non-zero entry to zero. Remove the entry */
      /* this might remove an entire column, or leave just a bound. No
	 nice solution for that yet */
      
      /* Shift the matrix */
      lastelm = lp->non_zeros; 
      for(i = elmnr; i < lastelm ; i++)
	lp->mat[i] = lp->mat[i + 1];
      for(i = Column; i <= lp->columns; i++)
	lp->col_end[i]--;
      
      lp->non_zeros--;
    }
  }
  else if(my_abs(Value) > lp->epsilon) {
    /* no existing entry. make new one only if not nearly zero */
    /* check if more space is needed for matrix */
    inc_mat_space(lp, 1);
    
    /* Shift the matrix */
    lastelm = lp->non_zeros; 
    for(i = lastelm; i > elmnr ; i--)
      lp->mat[i] = lp->mat[i - 1];
    for(i = Column; i <= lp->columns; i++)
      lp->col_end[i]++;
    
    /* Set new element */
    lp->mat[elmnr].row_nr = Row;
    
    if (lp->scaling_used) {
      if(lp->ch_sign[Row])
	lp->mat[elmnr].value = -Value * lp->scale[Row] * lp->scale[Column];
      else
	lp->mat[elmnr].value = Value * lp->scale[Row] * lp->scale[Column];
    }
    else /* no scaling */
      {
	if(lp->ch_sign[Row])
	  lp->mat[elmnr].value = -Value;
	else
	  lp->mat[elmnr].value = Value;
      }
    
    lp->row_end_valid = FALSE;
    
    lp->non_zeros++;
  }      
}

void set_obj_fn(lprec *lp, REAL *row)
{
  int i;
  for(i = 1; i <= lp->columns; i++)
    set_mat(lp, 0, i, row[i]);
}

void str_set_obj_fn(lprec *lp, char *row)
{
  int  i;
  REAL *arow;
  char *p, *newp;
  CALLOC(arow, lp->columns + 1);
  p = row;
  for(i = 1; i <= lp->columns; i++) {
    arow[i] = (REAL) strtod(p, &newp);
    if(p == newp)
      error("Bad string in str_set_obj_fn");
    else
      p = newp; 
  }
  set_obj_fn(lp, arow);
  free(arow);
}


void add_constraint(lprec *lp, REAL *row, short constr_type, REAL rh)
{
  matrec *newmat;
  int  i, j;
  int  elmnr;
  int  stcol;
  int  *addtoo;

  MALLOC(addtoo, lp->columns + 1);

    for(i = 1; i <= lp->columns; i++)
      if(row[i] != 0) {