lpkit.c

生成直角Steiner树的程序包

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

  if(constr_type == REL_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;
}

void str_add_constraint(lprec *lp,
			char *row_string,
			short constr_type,
			REAL rh)
{
  int  i;
  REAL *aRow;
  char *p, *newp;
  CALLOC(aRow, lp->columns + 1);
  p = row_string;
 
  for(i = 1; i <= lp->columns; i++) {
    aRow[i] = (REAL) strtod(p, &newp);
    if(p == newp)
      error("Bad string in str_add_constr");
    else
      p = newp; 
  }
  add_constraint(lp, aRow, constr_type, rh);
  free(aRow);
}

void add_rows(lprec *lp,
	      int ccnt,		/* # of new columns being added */
	      int rcnt,		/* # of new rows being added */
	      REAL *rh,         /* vector of rcnt right-hand-sides */
	      short *ctype,	/* vector of rcnt constraint types */
	      int *rmatbeg,	/* index of new row i's non-zeros */
	      int *rmatind,	/* column # of non-zero entries */
	      REAL *rmatval)	/* non-zero entries */
{
  int i, j, nzcnt;
  int *count;
  matrec **ptr1, **ptr2, *tmat;
  matrec *mp, *src, *dst, *end, *save;

  /* Handle new columns... */
  if(ccnt < 0 || rcnt < 0)
    error("Too few rows or columns.");

  if(rmatbeg[0] != 0)
    error("Invalid 'rmatbeg' array.");
  for(i = 0; i < rcnt; i++)
    if(rmatbeg[i+1] <= rmatbeg[i])
      error("Invalid 'rmatbeg' array.");
  nzcnt = rmatbeg[rcnt];

  for(i = 0; i < nzcnt; i++)
    if(rmatind[i] < 0 || rmatind[i] >= lp->columns + ccnt)
      error("Invalid 'rmatind' array.");

  lp->non_zeros += nzcnt;
  inc_mat_space(lp, 0);

  if(ccnt > 0)
    {
      lp->columns += ccnt;
      inc_col_space(lp);
      for(i = lp->sum + 1 - ccnt; i <= lp->sum; i++)
	{
	  lp->orig_lowbo[i] = 0;
	  lp->orig_upbo[i] = lp->infinite;
	}
    }

  if(rcnt > 0)
    {
      /* Handle new rows... */
      lp->rows += rcnt;
      lp->sum += rcnt;
      inc_row_space(lp);

      if(lp->scaling_used)
        {
	  /* shift scale */
	  for(i = lp->sum; i > lp->rows; i--)
	    lp->scale[i]=lp->scale[i-rcnt];
	  for(; i > lp->rows - rcnt; i--)
	    lp->scale[i] = 1;
	  if(lp->columns_scaled)
	    for(i = 0; i < nzcnt; i++)
	      rmatval[i] *= lp->scale[lp->rows + 1 + rmatind[i]];
	}

      if(lp->names_used)
	for(i = lp->rows + 1 - rcnt; i <= lp->rows; i++)
	  sprintf(lp->row_name[i], "r_%d", i);

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

      for(i = 0; i < rcnt; i++)
	switch(ctype[i])
	  {
	  case REL_LE:
	    lp->ch_sign[lp->rows + 1 - rcnt + i] = FALSE;
	    lp->orig_rh[lp->rows + 1 - rcnt + i] = rh[i];
	    lp->orig_upbo[lp->rows + 1 - rcnt + i] = lp->infinite;
	    lp->orig_lowbo[lp->rows + 1 - rcnt + i] = 0;
	    break;

	  case REL_GE:
	    lp->ch_sign[lp->rows + 1 - rcnt + i] = TRUE;
	    lp->orig_rh[lp->rows + 1 - rcnt + i] = - rh[i];
	    lp->orig_upbo[lp->rows + 1 - rcnt + i] = lp->infinite;
	    lp->orig_lowbo[lp->rows + 1 - rcnt + i] = 0;
	    for(j = rmatbeg[i]; j < rmatbeg[i+1]; j++)
	      rmatval[j] = -rmatval[j];
	    break;

	  case REL_EQ:
	    lp->ch_sign[lp->rows + 1 - rcnt + i] = FALSE;
	    lp->orig_rh[lp->rows + 1 - rcnt + i] = rh[i];
	    lp->orig_upbo[lp->rows + 1 - rcnt + i] = 0;
	    lp->orig_lowbo[lp->rows + 1 - rcnt + i] = 0;
	    break;

	  default:
	    fprintf(stderr, "Wrong constraint type\n");
	    exit (EXIT_FAILURE);
	  }

      for(i = 1; i <= lp->rows - rcnt; i++)
	if(lp->bas[i] > lp->rows - rcnt)
	  lp->bas[i] += rcnt;

      for(i = lp->rows + 1 - rcnt; i <= lp->rows; i++)
	{
	  lp->bas[i] = i;
	  lp->basis[i] = TRUE;
	  lp->lower[i] = TRUE;
	  lp->must_be_int[i] = FALSE;
	}

      lp->row_end_valid = FALSE;

      /* Now integrate the new coefficients into the matrix */
      MALLOC(count, lp->columns);
      MALLOC(ptr1, lp->columns+1);
      MALLOC(ptr2, lp->columns+1);
      MALLOC(tmat, nzcnt);

      /* reorganize as columns... */
      for(i = 0; i < lp->columns; i++)
	count[i] = 0;
      for(i = 0; i < nzcnt; i++)
	++(count[rmatind[i]]);

      mp = tmat;
      for(i = 0; i < lp->columns; i++)
	{
	  ptr1[i] = mp;
	  ptr2[i] = mp;
	  mp += count[i];
	}
      ptr1[lp->columns] = mp;
      ptr2[lp->columns] = mp;

      for(i = 0; i < rcnt; i++)
	for(j = rmatbeg[i]; j < rmatbeg[i+1]; j++)
	  {
	    mp = (ptr2[rmatind[j]])++;
	    mp->row_nr = lp->rows + 1 - rcnt + i;
	    mp->value = rmatval[j];
	  }

      /* merge new columns into matrix */
      dst = lp->mat + lp->non_zeros;
      for(i = lp->columns; i > 0; i--)
	{
	  /* Insert new rows of this column... */
	  save = lp->mat + lp->col_end[i];
	  lp->col_end[i] = dst - lp->mat;
	  src = ptr1[i];
	  end = ptr1[i-1];
	  while (src > end)
	    *--dst = *--src;

	  /* Copy the old rows of this column... */
	  src = save;
	  end = lp->mat + lp->col_end[i-1];
	  while (src > end)
	    *--dst = *--src;
	}

      free(tmat);
      free(ptr2);
      free(ptr1);
      free(count);
    }

  lp->eta_valid = FALSE;
}

void del_constraint(lprec *lp, int del_row)
{
  int i, j;
  unsigned elmnr;
  int startcol;

  if(del_row<1 || del_row>lp->rows)
    {
      fprintf(stderr, "There is no constraint nr. %d\n", del_row);
      exit(EXIT_FAILURE);
    }

  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)
      strcpy(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];
    lp->must_be_int[i] = lp->must_be_int[i + 1];
    if(lp->scaling_used)
      lp->scale[i] = lp->scale[i + 1];
  }

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

  lp->row_end_valid = FALSE;
  lp->eta_valid     = FALSE;
  lp->basis_valid   = FALSE; 
}

void delete_row_set(lprec *lp, int *row_flags)
{
  int i, j, k;
  unsigned elmnr;
  int startcol;
  int * renum;
  int firstrow;
  int num_del;

  /* Find lowest numbered row to delete. */
  firstrow = -1;
  for (i = 0; i <= lp->rows; i++)
    if(row_flags[i])
      {
	firstrow = i;
	break;
      }

  if(firstrow < 0)
    return;

  /* Build a map to renumber the rows and columns. */
  CALLOC(renum, lp->sum + 1);
  num_del = 0;
  j = 0;
  for(i = 0; i <= lp->rows; i++)
    if(row_flags[i])
      {
	renum[i] = -1;
	++num_del;
      }
    else
      renum[i] = j++;
  for(; i <= lp->sum; i++)
    renum[i] = j++;

  elmnr=0;
  startcol=0;

  for(i = 1; i <= lp->columns; i++)
    {
      for(j=startcol; j < lp->col_end[i]; j++)
	{
	  k = renum[lp->mat[j].row_nr];
	  if(k >= 0)
	    {
	      lp->mat[elmnr].row_nr = k;
	      lp->mat[elmnr].value  = lp->mat[j].value;
	      elmnr++;
	    }
	  else
	    lp->non_zeros--;
	}
      startcol=lp->col_end[i];
      lp->col_end[i]=elmnr;
    }
  for(i = firstrow + 1; i <= lp->rows; i++)
    {
      k = renum[i];
      if (k < 0) continue;
      lp->orig_rh[k] = lp->orig_rh[i];
      lp->ch_sign[k] = lp->ch_sign[i];
      if(lp->names_used)
        strcpy(lp->row_name[k], lp->row_name[i]);
    }
  for(i = 1; i <= lp->rows; i++)
    {
      j = renum[i];
      if(j < 0)
	{
	  lp->basis[lp->bas[i]] = 0;
	  lp->lower[lp->bas[i]] = 1;
	  continue;
	}
      k = renum[lp->bas[i]];
      if (k < 0)
	{
	  lp->basis_valid = 0;
	  break;
	}
      lp->bas[j] = k;
    }

  for(i = firstrow + 1; i <= lp->sum; i++)
    {
      k = renum[i];
      if(k < 0) continue;
      lp->lower[k]=lp->lower[i];
      lp->basis[k]=lp->basis[i];
      lp->orig_upbo[k]=lp->orig_upbo[i];
      lp->orig_lowbo[k]=lp->orig_lowbo[i];
      lp->must_be_int[k]=lp->must_be_int[i];
      if(lp->scaling_used)
        lp->scale[k]=lp->scale[i];
    }

  free(renum);

  lp->rows -= num_del;
  lp->sum -= num_del;;

  lp->row_end_valid=FALSE;
  lp->eta_valid=FALSE;
}

void add_lag_con(lprec *lp, REAL *row, short con_type, REAL rhs)
{
  int i;
  REAL sign;
  if(con_type == REL_LE || con_type == REL_EQ)
    sign = 1;
  else if(con_type == REL_GE)
    sign = -1;
  else
    error("con_type not implemented\n");

  lp->nr_lagrange++;
  if(lp->nr_lagrange == 1) {
    CALLOC(lp->lag_row, lp->nr_lagrange);
    CALLOC(lp->lag_rhs, lp->nr_lagrange);
    CALLOC(lp->lambda, lp->nr_lagrange);
    CALLOC(lp->lag_con_type, lp->nr_lagrange);
  }
  else {
    REALLOC(lp->lag_row, lp->nr_lagrange);
    REALLOC(lp->lag_rhs, lp->nr_lagrange);
    REALLOC(lp->lambda, lp->nr_lagrange);
    REALLOC(lp->lag_con_type, lp->nr_lagrange);
  }
  CALLOC(lp->lag_row[lp->nr_lagrange-1], lp->columns+1);
  lp->lag_rhs[lp->nr_lagrange-1] = rhs * sign;
  for( i = 1; i <= lp->columns; i++)
    lp->lag_row[lp->nr_lagrange-1][i] = row[i] * sign;
  lp->lambda[lp->nr_lagrange-1] = 0;
  lp->lag_con_type[lp->nr_lagrange-1]=(con_type == REL_EQ);
}

void str_add_lag_con(lprec *lp, char *row, short con_type, REAL rhs)
{
  int  i;
  REAL *a_row;
  char *p, *new_p;
  CALLOC(a_row, lp->columns + 1);
  p = row;
 
  for(i = 1; i <= lp->columns; i++) {
    a_row[i] = (REAL) strtod(p, &new_p);
    if(p == new_p)
      error("Bad string in str_add_lag_con");
    else
      p = new_p; 
  }
  add_lag_con(lp, a_row, con_type, rhs);
  free(a_row);
}


void add_column(lprec *lp, REAL *column)
{
  int i, elmnr;

  /* if the column has only one entry, this should be handled as a bound, but
     this currently is not the case */

  lp->columns++;
  lp->sum++;
  inc_col_space(lp);
  inc_mat_space(lp, lp->rows + 1);

  if(lp->scaling_used) {
    for(i = 0; i <= lp->rows; i++)
      column[i] *= lp->scale[i];
    lp->scale[lp->sum] = 1;
  }

  elmnr = lp->col_end[lp->columns - 1];
  for(i = 0 ; i <= lp->rows ; i++)
    if(column[i] != 0) {
      lp->mat[elmnr].row_nr = i;
      if(lp->ch_sign[i])
	lp->mat[elmnr].value = -column[i];
      else
	lp->mat[elmnr].value = column[i];
      lp->non_zeros++;
      elmnr++;
    }
  lp->col_end[lp->columns] = elmnr;
  lp->orig_lowbo[lp->sum] = 0;
  lp->orig_upbo[lp->sum] = lp->infinite;
  lp->lower[lp->sum] = TRUE;
  lp->basis[lp->sum] = FALSE;
  lp->must_be_int[lp->sum] = FALSE;
  if(lp->names_used)
    sprintf(lp->col_name[lp->columns], "var_%d", lp->columns);
 
  lp->row_end_valid = FALSE;
}

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

void del_column(lprec *lp, int column)
{
  int i, j, from_elm, to_elm, elm_in_col;
  if(column > lp->columns || column < 1)
    error("Column out of range in del_column");
  for(i = 1; i <= lp->rows; i++) {
    if(lp->bas[i] == lp->rows + column)
      lp->basis_valid = FALSE;
    else if(lp->bas[i] > lp->rows + column)
      lp->bas[i]--;
  }
  for(i = lp->rows + column; i < lp->sum; i++) {
    if(lp->names_used)
      strcpy(lp->col_name[i - lp->rows], lp->col_name[i - lp->rows + 1]);
    lp->must_be_int[i] = lp->must_be_int[i + 1];
    lp->orig_upbo[i] = lp->orig_upbo[i + 1];
    lp->orig_lowbo[i] = lp->orig_lowbo[i + 1];
    lp->upbo[i] = lp->upbo[i + 1];
    lp->lowbo[i] = lp->lowbo[i + 1];
    lp->basis[i] = lp->basis[i + 1];
    lp->lower[i] = lp->lower[i + 1];
    if(lp->scaling_used)
      lp->scale[i] = lp->scale[i + 1];
  }
  for(i = 0; i < lp->nr_lagrange; i++)
    for(j = column; j <= lp->columns; j++)
      lp->lag_row[i][j] = lp->lag_row[i][j+1];
  to_elm = lp->col_end[column-1];
  from_elm = lp->col_end[column];
  elm_in_col = from_elm-to_elm;
  for(i = from_elm; i < lp->non_zeros; i++) {
    lp->mat[to_elm] = lp->mat[i];
    to_elm++;
  }
  for(i = column; i < lp->columns; i++)
    lp->col_end[i] = lp->col_end[i + 1] - elm_in_col;
  lp->non_zeros -= elm_in_col;
  lp->row_end_valid = FALSE;
  lp->eta_valid = FALSE;

  lp->sum--;
  lp->columns--;
}

void set_upbo(lprec *lp, int column, REAL value)
{
  if(column > lp->columns || column < 1)
    error("Column out of range");