lpkit.c

利用c语言编写

/* It is however faster ... */
static int find_mat_elm1(lprec *lp, int row, int column, int *insertpos)
{
  int low, high, mid, item;

  (*insertpos) = -1;

  if(row < 0 || row > lp->rows) {
    report(lp, IMPORTANT, "find_mat_elm: Row %d out of range", row);
    return(-1);
  }
  if(column < 1 || column > lp->columns) {
    report(lp, IMPORTANT, "find_mat_elm: Column %d out of range", column);
    return(-1);
  }

  low = lp->col_end[column-1];
  (*insertpos) = low;
  high = lp->col_end[column] - 1;
  if(low > high)
    return(-2);

 /* Do binary search logic */
  mid = (low+high) / 2;
  item = lp->mat[mid].row_nr;
  while(high - low > LINEARSEARCH) {
    if(item < row) {
      low = mid + 1;
      mid = (low+high) / 2;
      item = lp->mat[mid].row_nr;
    }
    else if(item > row) {
      high = mid - 1;
      mid = (low+high) / 2;
      item = lp->mat[mid].row_nr;
    }
    else {
      low = mid;
      high = mid;
    }
  }

 /* Do linear scan search logic */
  if((high > low) && (high - low <= LINEARSEARCH)) {
    item = lp->mat[low].row_nr;
    while((low < high) && (item < row)) {
      low++;
      item = lp->mat[low].row_nr;
    }
    if(item == row)
      high = low;
  }

  (*insertpos) = low;
  if((low == high) && (row == item))
    return(low);
  else {
    if((low < lp->col_end[column]) && (lp->mat[low].row_nr < row))
      (*insertpos)++;
    return(-2);
  }
}

#if defined _DEBUG || defined DEBUG
/* derrived from the original code because find_mat_elm from KJELL not always returns correct results ... */
static int find_mat_elm2(lprec *lp, int row, int column, int *insertpos)
{
  int elmnr, col_end;
  matrec *matel;

  if(row < 0 || row > lp->rows) {
    report(lp, SEVERE, "find_mat_elm: Row %d out of range", row);
    return(-1);
  }
  if(column < 1 || column > lp->columns) {
    report(lp, SEVERE, "find_mat_elm: Column %d out of range", column);
    return(-1);
  }

  elmnr = lp->col_end[column - 1];
  matel = lp->mat + elmnr;
  col_end = lp->col_end[column];
  while(elmnr < col_end && matel->row_nr < /* != */ row) {
    matel++;
    elmnr++;
  }

  *insertpos = elmnr;

  if((elmnr == col_end) || (matel->row_nr > row))
    elmnr = -2;

  return(elmnr);
}

static int find_mat_elm(lprec *lp, int row, int column, int *insertpos)
{
  int elmnr1, elmnr2, insertpos1, insertpos2;

  /* check spare matrix order */

  {
    int column,elmnr,col_end,row;
    matrec *matel;

    for (column=1;column<=lp->columns;column++) {
	elmnr = lp->col_end[column - 1];
	matel = lp->mat + elmnr;
	col_end = lp->col_end[column];
	row=-1;
	while(elmnr < col_end) {
	    if ((row!=-1) && (matel->row_nr<=row))
	      fprintf(stderr,"Error\n");
	    row=matel->row_nr;
	    matel++;
	    elmnr++;
	}
    }
  }


  elmnr1 = find_mat_elm1(lp, row, column, &insertpos1);
  elmnr2 = find_mat_elm2(lp, row, column, &insertpos2);
  if ((elmnr1 != elmnr2) || (insertpos1 != insertpos2)) {
    fprintf(stderr,"Error: find_mat_elm1 returns different result from find_mat_elm2\n");
  }

  *insertpos = insertpos1;
  return(elmnr1);
}
#else

#define find_mat_elm find_mat_elm1

#endif

int set_matrix(lprec *lp, int Row, int Column, REAL Value, MYBOOL doscale)
{
  int elmnr, lastelm, i;

  /* This function is 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 */

  /* find out if we already have such an entry, or return insertion point */
  i = find_mat_elm(lp, Row, Column, &elmnr);
  if(i == -1)
    return(FALSE);

  if (lp->basis[Column] == TRUE && Row > 0)
    lp->basis_valid = FALSE;
  lp->eta_valid = FALSE;

  if((elmnr != lp->col_end[Column]) && (lp->mat[elmnr].row_nr == Row)) {
    /* there is an existing entry */
    if(my_abs(Value) > lp->epsel) { /* we replace it by something non-zero */
      if(doscale && lp->scaling_used) {
	if(lp->ch_sign[Row])
	  lp->mat[elmnr].value = -Value * lp->scale[Row] * lp->scale[lp->rows + Column];
	else
	  lp->mat[elmnr].value = Value * lp->scale[Row] * lp->scale[lp->rows + 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 up tail end of the matrix */
      lp->non_zeros--;           /* Don't cross array border - Moved by KE */
      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]--;
    }
  }
  else if(my_abs(Value) > lp->epsel) {
    /* no existing entry. make new one only if not nearly zero */
    /* check if more space is needed for matrix */
    if (!inc_mat_space(lp, 1))
      return(FALSE);

    /* Shift down tail end of the matrix by one */
    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(doscale && lp->scaling_used) {
      if(lp->ch_sign[Row])
	lp->mat[elmnr].value = -Value * lp->scale[Row] * lp->scale[lp->rows + Column];
      else
	lp->mat[elmnr].value = Value * lp->scale[Row] * lp->scale[lp->rows + 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++;
  }
  return(TRUE);
}

int set_mat(lprec *lp, int Row, int Column, REAL Value)
{
  return(set_matrix(lp, Row, Column, Value, lp->scaling_used));
}

int set_obj_fn(lprec *lp, REAL *row)
{
  int i;

  for(i = 1; i <= lp->columns; i++)
    if (!set_mat(lp, 0, i, row[i]))
      return(FALSE);
  return(TRUE);
}

int str_set_obj_fn(lprec *lp, char *row)
{
  int  i, ok = TRUE;
  REAL *arow;
  char *p, *newp;

  if (CALLOC(arow, lp->columns + 1) == NULL) {
    lp->spx_status = OUT_OF_MEMORY;
    ok = FALSE;
  }
  else {
    p = row;
    for(i = 1; i <= lp->columns; i++) {
      arow[i] = (REAL) strtod(p, &newp);
      if(p == newp) {
	report(lp, IMPORTANT, "str_set_obj_fn: Bad string");
	ok = FALSE;
	lp->spx_status = IGNORED;
	break;
      }
      else
	p = newp;
    }
    if (ok)
      if (!set_obj_fn(lp, arow))
        ok = FALSE;
    free(arow);
  }
  return(ok);
}

int add_constraint(lprec *lp, REAL *row, short constr_type, REAL rh)
{
  int    i, j, stcol, ok = TRUE;
  int    elmnr, orignr, newnr;
  MYBOOL   *addto;

  if(!(constr_type == LE || constr_type == GE || constr_type == EQ)) {
    report(lp, SEVERE, "add_constraint: Invalid %d constraint type", constr_type);
    return(FALSE);
  }

  if (MALLOC(addto, lp->columns + 1) == NULL) {
    lp->spx_status = OUT_OF_MEMORY;
    return(FALSE);
  }

  newnr = 0;
  for(i = 1; i <= lp->columns; i++)
    if(row[i] != 0) {
      addto[i] = TRUE;
      newnr++;
    }
    else
      addto[i] = FALSE;

  orignr = lp->non_zeros;
  lp->non_zeros += newnr;

  if (!inc_mat_space(lp, 0)) {
    free(addto);
    return(FALSE);
  }
  lp->rows++;
  lp->sum++;
  if (!inc_row_space(lp)) {
    free(addto);
    return(FALSE);
  }

  if(lp->scale != NULL) {
    /* shift column scale */
    for(i = lp->sum; i > lp->rows; i--)
      lp->scale[i] = lp->scale[i - 1];
    /* insert default new row scalar */
    lp->scale[lp->rows] = 1;
  }

  if(lp->scaling_used && lp->columns_scaled)
    for(i = 1; i <= lp->columns; i++)
      row[i] *= lp->scale[lp->rows + lp->nr_lagrange + i];

  if(constr_type == GE)
    lp->ch_sign[lp->rows] = TRUE;
  else
    lp->ch_sign[lp->rows] = FALSE;

  elmnr = lp->non_zeros - 1;
  orignr--;
  for(j = lp->columns; j > 0; j--) {
    stcol = lp->col_end[j] - 1;
    lp->col_end[j] = elmnr + 1;
   /* Add a new non-zero entry */
    if(addto[j]) {
      lp->mat[elmnr].row_nr = lp->rows;
      if(lp->ch_sign[lp->rows])
  	lp->mat[elmnr].value = -row[j];
      else
  	lp->mat[elmnr].value = row[j];
      elmnr--;
    }

   /* Check if we are finished */
    if(elmnr <= orignr)
      break;

   /* Shift previous column entries down */
    for(i = stcol; i >= lp->col_end[j-1]; i--) {
      lp->mat[elmnr] = lp->mat[orignr];
      orignr--;
      elmnr--;
    }
  }

  free(addto);

  for(i = lp->sum; i > lp->rows; i--) {
    lp->orig_upbo[i]   = lp->orig_upbo[i - 1];
    lp->orig_lowbo[i]  = lp->orig_lowbo[i - 1];
    lp->basis[i]       = lp->basis[i - 1];
    lp->lower[i]       = lp->lower[i - 1];
  }

  /* changed from i <= lp->rows to i < lp->rows, MB */
  for(i = 1 ; i < lp->rows; i++)
    if(lp->bas[i] >= lp->rows)
      lp->bas[i]++;

  if(constr_type == LE || constr_type == GE) {
    lp->orig_upbo[lp->rows] = lp->infinite;
  }
  else if(constr_type == EQ) {
    lp->orig_upbo[lp->rows] = 0;
  }
  else {
    report(lp, SEVERE, "add_constraint: Wrong constraint type");
    return(FALSE);
  }

  lp->orig_lowbo[lp->rows] = 0;

  if(constr_type == GE && rh != 0)
    lp->orig_rh[lp->rows] = -rh;
  else
    lp->orig_rh[lp->rows] = rh;

  lp->row_end_valid = FALSE;

  lp->bas[lp->rows] = lp->rows;
  lp->basis[lp->rows] = TRUE;
  lp->lower[lp->rows] = TRUE;
  lp->eta_valid = FALSE;
  return(ok);
}

int str_add_constraint(lprec *lp,
			char *row_string,
			short constr_type,
			REAL rh)
{
  int  i, ok = TRUE;
  REAL *aRow;
  char *p, *newp;

  if (CALLOC(aRow, lp->columns + 1) == NULL) {
    ok = FALSE;
    lp->spx_status = OUT_OF_MEMORY;
  }
  else {
    p = row_string;

    for(i = 1; i <= lp->columns; i++) {
      aRow[i] = (REAL) strtod(p, &newp);
      if(p == newp) {
	report(lp, IMPORTANT, "str_add_constraint: Bad string");
	ok = FALSE;
	lp->spx_status = IGNORED;
	break;
      }
      else
	p = newp;
    }
    if(ok)
      if (!add_constraint(lp, aRow, constr_type, rh))
        ok = FALSE;
    free(aRow);
  }
  return(ok);
}

int del_constraint(lprec *lp, int del_row)
{
  int i, j;
#if 0
  int k
#endif
  int elmnr;
  int startcol;

  if(del_row<1 || del_row>lp->rows) {
    report(lp, IMPORTANT, "del_constraint: Attempt to delete constraint %d that does not exist", del_row);
    return(FALSE);
  }

#if 0
  k = lp->var_to_orig[del_row];
#endif

  elmnr = 0;
  startcol = 0;
  for(i = 1; i <= lp->columns; i++) {
    for(j = startcol; j < lp->col_end[i]; j++) {
      if(lp->mat[j].row_nr != del_row) {
	lp->mat[elmnr] = lp->mat[j];
	if(lp->mat[elmnr].row_nr > del_row)
	  lp->mat[elmnr].row_nr--;
	elmnr++;
      }
      else
	lp->non_zeros--;
    }
    startcol = lp->col_end[i];
    lp->col_end[i] = elmnr;
  }
  for(i = del_row; i < lp->rows; i++) {
    lp->orig_rh[i] = lp->orig_rh[i + 1];
    lp->ch_sign[i] = lp->ch_sign[i + 1];
    lp->bas[i] = lp->bas[i + 1];
    if(lp->names_used)
      lp->row_name[i] = lp->row_name[i + 1];
  }
  for(i = 1; i < lp->rows; i++)
    if(lp->bas[i] > del_row)
      lp->bas[i]--;

  for(i = del_row; i < lp->sum; i++) {
    lp->lower[i] = lp->lower[i + 1];
    lp->basis[i] = lp->basis[i + 1];
    lp->orig_upbo[i] = lp->orig_upbo[i + 1];
    lp->orig_lowbo[i] = lp->orig_lowbo[i + 1];
    if(lp->scaling_used)
      lp->scale[i] = lp->scale[i + 1];
#if 0
    lp->var_to_orig[i] = lp->var_to_orig[i + 1];
#endif
  }

#if 0
  lp->orig_to_var[k] = -del_row;
  for(i = k + 1; i <= lp->orig_rows; i++)
    if(lp->orig_to_var[i] > del_row)
      lp->orig_to_var[i]--;
#endif

  lp->rows--;
  lp->sum--;

  lp->row_end_valid = FALSE;
  lp->eta_valid     = FALSE;
  lp->basis_valid   = FALSE;
  return(TRUE);
}

int add_lag_con(lprec *lp, REAL *row, short con_type, REAL rhs)
{
  int i, ok = TRUE;
  REAL sign;

  if(con_type == LE || con_type == EQ)
    sign = 1;
  else if(con_type == GE)
    sign = -1;
  else {
    report(lp, IMPORTANT, "add_lag_con: con_type %d not implemented", con_type);
    return(FALSE);
  }

  lp->nr_lagrange++;
  if(lp->nr_lagrange == 1) {
    lp->lag_row = NULL;
    lp->lag_rhs = NULL;
    lp->lambda = NULL;
    lp->lag_con_type = NULL;
    if ((CALLOC(lp->lag_row, lp->nr_lagrange) == NULL) ||
        (CALLOC(lp->lag_rhs, lp->nr_lagrange) == NULL) ||
        (CALLOC(lp->lambda, lp->nr_lagrange) == NULL) ||
        (CALLOC(lp->lag_con_type, lp->nr_lagrange) == NULL)
       ) {
      FREE(lp->lag_con_type);
      FREE(lp->lambda);
      FREE(lp->lag_rhs);
      FREE(lp->lag_row);
      lp->spx_status = OUT_OF_MEMORY;
      ok = FALSE;
    }
  }
  else {
    if ((REALLOC(lp->lag_row, lp->nr_lagrange) == NULL) ||
        (REALLOC(lp->lag_rhs, lp->nr_lagrange) == NULL) ||