solve.c

来自「生成直角Steiner树的程序包」· C语言 代码 · 共 2,233 行 · 第 1/4 页

static int solvelp(lprec *lp)
{
  int    i, singular_count, lost_feas_count;
  short  feasible;
  REAL	 x;

  if(lp->do_presolve)
    presolve(lp);

  lp->iter = 0;

  singular_count = 0;
  lost_feas_count = 0;

  for(;;) {
    /* Check whether we are feasible or infeasible. */
    feasible = TRUE;
    for(i = 1; i <= lp->rows; i++) {
      x = lp->rhs[i];
      if((x < 0) || (x > lp->upbo[lp->bas[i]])) {
	feasible = FALSE;
	break;
      }
    }

    if(feasible) {
      primloop(lp);
    }
    else {
      dualloop(lp);
      if(Status == SWITCH_TO_PRIMAL)
	primloop(lp);
    }

    if(Status == SINGULAR_BASIS) {
      ++singular_count;
      if(singular_count >= 5) {
	printf("%% SINGULAR BASIS!  Too many singularities - aborting.\n");
	Status = FAILURE;
	break;
      }
      printf("%% SINGULAR BASIS!  Will attempt to recover.\n");
      /* Singular pivots are simply skipped by the inversion, leaving */
      /* a row's slack var in the basis instead of the singular problem */
      /* var.  This basis could be feasible or infeasible.  Check how to */
      /* restart. */
      continue;
    }
    else if(Status == LOST_PRIMAL_FEASIBILITY) {
      ++lost_feas_count;
      if(lost_feas_count >= 5) {
	printf("%% LOST PRIMAL FEASIBILITY too many times, aborting.\n");
	Status = FAILURE;
	break;
      }
      printf("%% LOST PRIMAL FEASIBILITY!  Recovering.\n");
      continue;
    }
    else {
      /* One of the "normal" status codes -- we are done! */
      break;
    }
  }

  lp->total_iter += lp->iter;

  return(Status);
} /* solvelp */


static short is_int(lprec *lp, int i)
{
  REAL   value, error;

  value = lp->solution[i];
  error = value - (REAL)floor((double)value);

  if(error < lp->epsilon)
    return(TRUE);

  if(error > (1 - lp->epsilon))
    return(TRUE);

  return(FALSE);
} /* is_int */


static void construct_solution(lprec *lp)
{
  int    i, j, basi;
  REAL   f;

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

  lp->solution[0] = -lp->orig_rh[0];

  if(lp->scaling_used) {
    lp->solution[0] /= lp->scale[0];

    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i] * lp->scale[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += lp->rhs[i] * lp->scale[basi];
    }
    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i] * lp->scale[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	  lp->solution[lp->mat[i].row_nr] += (f / lp->scale[lp->rows+j])
	    * (lp->mat[i].value / lp->scale[lp->mat[i].row_nr]);
    }

    for(i = 0; i <= lp->rows; i++) {
      if(my_abs(lp->solution[i]) < lp->epsb)
	lp->solution[i] = 0;
      else if(lp->ch_sign[i])
	lp->solution[i] = -lp->solution[i];
    }
  }
  else { /* no scaling */
    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += lp->rhs[i];
    }

    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	  lp->solution[lp->mat[i].row_nr] += f * lp->mat[i].value;
    }

    for(i = 0; i <= lp->rows; i++) {
      if(my_abs(lp->solution[i]) < lp->epsb)
	lp->solution[i] = 0;
      else if(lp->ch_sign[i])
	lp->solution[i] = -lp->solution[i];
    }
  }
} /* construct_solution */

static void calculate_duals(lprec *lp)
{
  int i;

  /* initialize */
  lp->duals[0] = 1;
  for(i = 1; i <= lp->rows; i++)
    lp->duals[i] = 0;

  btran(lp, lp->duals);

  if(lp->scaling_used)
    for(i = 1; i <= lp->rows; i++)
      lp->duals[i] *= lp->scale[i] / lp->scale[0];

  /* the dual values are the reduced costs of the slacks */
  /* When the slack is at its upper bound, change the sign. */
  for(i = 1; i <= lp->rows; i++) {
    if(lp->basis[i])
      lp->duals[i] = 0;
    /* added a test if variable is different from 0 because sometime you get
       -0 and this is different from 0 on for example INTEL processors (ie 0
       != -0 on INTEL !) PN */
    else if((lp->ch_sign[0] == lp->ch_sign[i]) && lp->duals[i])
      lp->duals[i] = - lp->duals[i];
  }
} /* calculate_duals */

static void check_if_less(REAL x,
			  REAL y,
			  REAL value)
{
  if(x >= y) {
    printf("%% Error: new upper or lower bound is not more restrictive\n");
    printf("%% bound 1: %g, bound 2: %g, value: %g\n",
	    (double)x, (double)y, (double)value);
    /* exit(EXIT_FAILURE); */
  }
}

static void check_solution(lprec *lp,
			   REAL *upbo,
			   REAL *lowbo)
{
  int i;

  /* check if all solution values are within the bounds, but allow some margin
     for numerical errors */

#define CHECK_EPS 1e-2

  if(lp->columns_scaled)
    for(i = lp->rows + 1; i <= lp->sum; i++) {
      if(lp->solution[i] < lowbo[i] * lp->scale[i] - CHECK_EPS) {
	printf("%% Error: variable %d (%s) has a solution (%g) smaller than its lower bound (%g)\n",
		i - lp->rows, lp->col_name[i - lp->rows],
		(double)lp->solution[i], (double)lowbo[i] * lp->scale[i]);
	/* abort(); */
      }

      if(lp->solution[i] > upbo[i] * lp->scale[i] + CHECK_EPS) {
	printf("%% Error: variable %d (%s) has a solution (%g) larger than its upper bound (%g)\n",
		i - lp->rows, lp->col_name[i - lp->rows],
		(double)lp->solution[i], (double)upbo[i] * lp->scale[i]);
	/* abort(); */
      }
    }
  else /* columns not scaled */
    for(i = lp->rows + 1; i <= lp->sum; i++) {
      if(lp->solution[i] < lowbo[i] - CHECK_EPS) {
	printf("%% Error: variable %d (%s) has a solution (%g) smaller than its lower bound (%g)\n",
		i - lp->rows, lp->col_name[i - lp->rows],
		(double)lp->solution[i], (double)lowbo[i]);
	/* abort(); */
      }

      if(lp->solution[i] > upbo[i] + CHECK_EPS) {
	printf("%% Error: variable %d (%s) has a solution (%g) larger than its upper bound (%g)\n",
		i - lp->rows, lp->col_name[i - lp->rows],
		(double)lp->solution[i], (double)upbo[i]);
	/* abort(); */
      }
    }
} /* check_solution */


static int milpsolve(lprec *lp,
		     REAL   *upbo,
		     REAL   *lowbo,
		     short  *sbasis,
		     short  *slower,
		     int    *sbas,
		     int     recursive)
{
  int i, j, failure, notint, is_worse;
  REAL theta, tmpreal;

  if(Break_bb)
    return(BREAK_BB);

  Level++;
  lp->total_nodes++;

  if(Level > lp->max_level)
    lp->max_level = Level;

  debug_print(lp, "starting solve");

  /* make fresh copies of upbo, lowbo, rh as solving changes them */
  memcpy(lp->upbo,  upbo,    (lp->sum + 1)  * sizeof(REAL));
  memcpy(lp->lowbo, lowbo,   (lp->sum + 1)  * sizeof(REAL));
  memcpy(lp->rh,    lp->orig_rh, (lp->rows + 1) * sizeof(REAL));

  /* make shure we do not do memcpy(lp->basis, lp->basis ...) ! */
  if(recursive) {
    memcpy(lp->basis, sbasis,  (lp->sum + 1)  * sizeof(short));
    memcpy(lp->lower, slower,  (lp->sum + 1)  * sizeof(short));
    memcpy(lp->bas,   sbas,    (lp->rows + 1) * sizeof(int));
  }

  if(lp->anti_degen) { /* randomly disturb bounds */
    for(i = 1; i <= lp->columns; i++) {
      tmpreal = (REAL) (rand() % 100 * 0.00001);
      if(tmpreal > lp->epsb)
	lp->lowbo[i + lp->rows] -= tmpreal;
      tmpreal = (REAL) (rand() % 100 * 0.00001);
      if(tmpreal > lp->epsb)
	lp->upbo[i + lp->rows] += tmpreal;
    }
    lp->eta_valid = FALSE;
  }

  if(!lp->eta_valid) {
    /* transform to all lower bounds to zero */
    for(i = 1; i <= lp->columns; i++)
      if((theta = lp->lowbo[lp->rows + i]) != 0) {
	if(lp->upbo[lp->rows + i] < lp->infinite)
	  lp->upbo[lp->rows + i] -= theta;
	for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	  lp->rh[lp->mat[j].row_nr] -= theta * lp->mat[j].value;
      }
    invert(lp);
    lp->eta_valid = TRUE;
  }

  failure = solvelp(lp);

  if(lp->anti_degen && (failure == OPTIMAL)) {
    /* restore to original problem, solve again starting from the basis found
       for the disturbed problem */

    /* restore original problem */
    memcpy(lp->upbo,  upbo,        (lp->sum + 1)  * sizeof(REAL));
    memcpy(lp->lowbo, lowbo,       (lp->sum + 1)  * sizeof(REAL));
    memcpy(lp->rh,    lp->orig_rh, (lp->rows + 1) * sizeof(REAL));

    /* transform to all lower bounds zero */
    for(i = 1; i <= lp->columns; i++)
      if((theta = lp->lowbo[lp->rows + i] != 0)) {
	if(lp->upbo[lp->rows + i] < lp->infinite)
	  lp->upbo[lp->rows + i] -= theta;
	for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	  lp->rh[lp->mat[j].row_nr] -= theta * lp->mat[j].value;
      }
    invert(lp);
    lp->eta_valid = TRUE;
    failure = solvelp(lp); /* and solve again */
  }

  if(failure != OPTIMAL)
    debug_print(lp, "this problem has no solution, it is %s",
		(failure == UNBOUNDED) ? "unbounded" : "infeasible");

  if(failure == INFEASIBLE && lp->verbose)
    printf("%% level %d INF\n", Level);

  if(failure == OPTIMAL) { /* there is a good solution */
    construct_solution(lp);

    /* because of reports of solution > upbo */
    /* check_solution(lp, upbo, lowbo); get too many hits ?? */

    debug_print(lp, "a solution was found");
    debug_print_solution(lp);

    /* if this solution is worse than the best sofar, this branch must die */

    /* if we can only have integer OF values, we might consider requiring to
       be at least 1 better than the best sofar, MB */

    if(lp->maximise)
      is_worse = lp->solution[0] <= lp->best_solution[0];
    else /* minimising! */
      is_worse = lp->solution[0] >= lp->best_solution[0];

    if(is_worse) {
      if(lp->verbose)
	printf("%% level %d OPT NOB value %g bound %g\n",
		Level, (double)lp->solution[0],
		(double)lp->best_solution[0]);
      debug_print(lp, "but it was worse than the best sofar, discarded");
      Level--;
      return(MILP_FAIL);
    }

    /* check if solution contains enough ints */
    if(lp->bb_rule == FIRST_NI) {
      for(notint = 0, i = lp->rows + 1;
	  i <= lp->sum && notint == 0;
	  i++) {
	if(lp->must_be_int[i] && !is_int(lp, i)) {
	  if(lowbo[i] == upbo[i]) { /* this var is already fixed */
	    printf("%% Warning: integer var %d is already fixed at %d, but has non-integer value %g\n",
		    i - lp->rows, (int)lowbo[i],
		    (double)lp->solution[i]);
	    printf("%% Perhaps the -e option should be used\n");
	  }
	  else
	    notint = i;
	}
      }
    }
    if(lp->bb_rule == RAND_NI) {
      int nr_not_int, select_not_int;
      nr_not_int = 0;

      for(i = lp->rows + 1; i <= lp->sum; i++)
	if(lp->must_be_int[i] && !is_int(lp, i))
	  nr_not_int++;

      if(nr_not_int == 0)
	notint = 0;
      else {
	select_not_int = (rand() % nr_not_int) + 1;
	i = lp->rows + 1;
	while(select_not_int > 0) {
	  if(lp->must_be_int[i] && !is_int(lp, i))
	    select_not_int--;
	  i++;
	}
	notint = i - 1;
      }
    }

    if(lp->verbose) {
      if(notint)
	printf("%% level %d OPT     value %g\n", Level,
		(double)lp->solution[0]);
      else
	printf("%% level %d OPT INT value %g\n", Level,
		(double)lp->solution[0]);
    }

    if(notint) { /* there is at least one value not yet int */
      /* set up two new problems */
      REAL   *new_upbo, *new_lowbo;
      REAL   new_bound;
      short  *new_lower,*new_basis;
      int    *new_bas;
      int     resone, restwo;

      /* allocate room for them */
      MALLOC(new_upbo,  lp->sum + 1);
      MALLOC(new_lowbo, lp->sum + 1);
      MALLOC(new_lower, lp->sum + 1);
      MALLOC(new_basis, lp->sum + 1);
      MALLOC(new_bas,   lp->rows + 1);
      memcpy(new_upbo,  upbo,      (lp->sum + 1)  * sizeof(REAL));
      memcpy(new_lowbo, lowbo,     (lp->sum + 1)  * sizeof(REAL));
      memcpy(new_lower, lp->lower, (lp->sum + 1)  * sizeof(short));
      memcpy(new_basis, lp->basis, (lp->sum + 1)  * sizeof(short));
      memcpy(new_bas,   lp->bas,   (lp->rows + 1) * sizeof(int));

      if(lp->names_used)
	debug_print(lp, "not enough ints. Selecting var %s, val: %g",
		    lp->col_name[notint - lp->rows],
		    (double)lp->solution[notint]);
      else
	debug_print(lp,
		    "not enough ints. Selecting Var [%d], val: %g",
		    notint, (double)lp->solution[notint]);
      debug_print(lp, "current bounds:\n");
      debug_print_bounds(lp, upbo, lowbo);

      if(lp->floor_first) {
	new_bound = ceil(lp->solution[notint]) - 1;

	/* this bound might conflict */
	if(new_bound < lowbo[notint]) {
	  debug_print(lp,
		      "New upper bound value %g conflicts with old lower bound %g\n",
		      (double)new_bound, (double)lowbo[notint]);
	  resone = MILP_FAIL;
	}
	else { /* bound feasible */
	  check_if_less(new_bound, upbo[notint], lp->solution[notint]);
	  new_upbo[notint] = new_bound;
	  debug_print(lp, "starting first subproblem with bounds:");
	  debug_print_bounds(lp, new_upbo, lowbo);
	  lp->eta_valid = FALSE;
	  resone = milpsolve(lp, new_upbo, lowbo, new_basis, new_lower,
			     new_bas, TRUE);
	  lp->eta_valid = FALSE;
	}
	new_bound += 1;
	if(new_bound > upbo[notint]) {
	  debug_print(lp,
		      "New lower bound value %g conflicts with old upper bound %g\n",
		      (double)new_bound, (double)upbo[notint]);
	  restwo = MILP_FAIL;
	}
	else { /* bound feasible */
	  check_if_less(lowbo[notint], new_bound, lp->solution[notint]);
	  new_lowbo[notint] = new_bound;
	  debug_print(lp, "starting second subproblem with bounds:");
	  debug_print_bounds(lp, upbo, new_lowbo);
	  lp->eta_valid = FALSE;
	  restwo = milpsolve(lp, upbo, new_lowbo, new_basis, new_lower,
			     new_bas, TRUE);
	  lp->eta_valid = FALSE;
	}
      }
      else { /* take ceiling first */
	new_bound = ceil(lp->solution[notint]);
	/* this bound might conflict */
	if(new_bound > upbo[notint]) {
	  debug_print(lp,
		      "New lower bound value %g conflicts with old upper bound %g\n",
		      (double)new_bound, (double)upbo[notint]);
	  resone = MILP_FAIL;
	}
	else { /* bound feasible */
	  check_if_less(lowbo[notint], new_bound, lp->solution[notint]);
	  new_lowbo[notint] = new_bound;
	  debug_print(lp, "starting first subproblem with bounds:");
	  debug_print_bounds(lp, upbo, new_lowbo);
	  lp->eta_valid = FALSE;
	  resone = milpsolve(lp, upbo, new_lowbo, new_basis, new_lower,
			     new_bas, TRUE);
	  lp->eta_valid = FALSE;
	}
	new_bound -= 1;
	if(new_bound < lowbo[notint]) {
	  debug_print(lp,
		      "New upper bound value %g conflicts with old lower bound %g\n",
		      (double)new_bound, (double)lowbo[notint]);
	  restwo = MILP_FAIL;
	}
	else { /* bound feasible */
	  check_if_less(new_bound, upbo[notint], lp->solution[notint]);
	  new_upbo[notint] = new_bound;
	  debug_print(lp, "starting second subproblem with bounds:");
	  debug_print_bounds(lp, new_upbo, lowbo);
	  lp->eta_valid = FALSE;
	  restwo = milpsolve(lp, new_upbo, lowbo, new_basis, new_lower,
			     new_bas, TRUE);
	  lp->eta_valid = FALSE;
	}
      }
      if(resone && restwo) /* both failed and must have been infeasible */
	failure = INFEASIBLE;
      else
	failure = OPTIMAL;

      free(new_upbo);
      free(new_lowbo);
      free(new_basis);
      free(new_lower);
      free(new_bas);
    }
    else { /* all required values are int */
      debug_print(lp, "--> valid solution found");

      if(lp->maximise)
	is_worse = lp->solution[0] < lp->best_solution[0];
      else
	is_worse = lp->solution[0] > lp->best_solution[0];

      if(!is_worse) { /* Current solution better */
	if(lp->debug || (lp->verbose && !lp->print_sol))
	  printf("%% *** new best solution: old: %g, new: %g ***\n",
		  (double)lp->best_solution[0], (double)lp->solution[0]);
	memcpy(lp->best_solution, lp->solution, (lp->sum + 1) * sizeof(REAL));
	calculate_duals(lp);

	if(lp->print_sol)
	  print_solution(lp);

	if(lp->break_at_int) {
	  if(lp->maximise && (lp->best_solution[0] > lp->break_value))
	    Break_bb = TRUE;

	  if(!lp->maximise && (lp->best_solution[0] < lp->break_value))
	    Break_bb = TRUE;
	}
      }
    }
  }

  Level--;

  /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */