if(lp->solution[i] > upbo[i] * lp->scale[i] + CHECK_EPS) {
	fprintf(stderr,
		"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) {
	fprintf(stderr,
		"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) {
	fprintf(stderr,
		"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)
    fprintf(stderr, "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)
	fprintf(stderr, "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 */
	    fprintf(stderr,
		    "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]);
	    fprintf(stderr, "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)
	fprintf(stderr, "level %d OPT     value %g\n", Level,
		(double)lp->solution[0]);
      else
	fprintf(stderr, "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))
	  fprintf(stderr,
		  "*** 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. */
  return(failure);
} /* milpsolve */


int solve(lprec *lp)
{
  int result, i;

  lp->total_iter  = 0;
  lp->max_level   = 1;
  lp->total_nodes = 0;

  if(isvalid(lp)) {
    if(lp->maximise && lp->obj_bound == lp->infinite)
      lp->best_solution[0] = -lp->infinite;
    else if(!lp->maximise && lp->obj_bound == -lp->infinite)
      lp->best_solution[0] = lp->infinite;
    else
      lp->best_solution[0] = lp->obj_bound;

    Level = 0;

    if(!lp->basis_valid) {
      for(i = 0; i <= lp->rows; i++) {
	lp->basis[i] = TRUE;
	lp->bas[i]   = i;
      }

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

      for(i = 0; i <= lp->sum; i++)
	lp->lower[i] = TRUE;

      lp->basis_valid = TRUE;
    }

    lp->eta_valid = FALSE;
    Break_bb      = FALSE;
    result        = milpsolve(lp, lp->orig_upbo, lp->orig_lowbo, lp->basis,
			      lp->lower, lp->bas, FALSE); 
    return(result);
  }

  /* if we get here, isvalid(lp) failed. I suggest we return FAILURE */
  fprintf(stderr, "Error, the current LP seems to be invalid\n");
  return(FAILURE);
} /* solve */

int lag_solve(lprec *lp, REAL start_bound, int num_iter, short verbose)
{
  int i, j, result, citer;
  short status, OrigFeas, AnyFeas, same_basis;
  REAL *OrigObj, *ModObj, *SubGrad, *BestFeasSol;
  REAL Zub, Zlb, Ztmp, pie;
  REAL rhsmod, Step, SqrsumSubGrad;
  int   *old_bas;
  short *old_lower;

  /* allocate mem */  
  MALLOC(OrigObj, lp->columns + 1);
  CALLOC(ModObj, lp->columns + 1);
  CALLOC(SubGrad, lp->nr_lagrange);
  CALLOC(BestFeasSol, lp->sum + 1);
  MALLOCCPY(old_bas, lp->bas, lp->rows + 1);
  MALLOCCPY(old_lower, lp->lower, lp->sum + 1);

  get_row(lp, 0, OrigObj);
 
  pie = 2;  

  if(lp->maximise) {
    Zub = DEF_INFINITE;
    Zlb = start_bound;
  }
  else {
    Zlb = -DEF_INFINITE;
    Zub = start_bound;
  }
  status   = RUNNING; 
  Step     = 1;
  OrigFeas = FALSE;
  AnyFeas  = FALSE;
  citer    = 0;

  for(i = 0 ; i < lp->nr_lagrange; i++)
    lp->lambda[i] = 0;

  while(status == RUNNING) {
    citer++;

    for(i = 1; i <= lp->columns; i++) {
      ModObj[i] = OrigObj[i];
      for(j = 0; j < lp->nr_lagrange; j++) {
	if(lp->maximise)
	  ModObj[i] -= lp->lambda[j] * lp->lag_row[j][i]; 
	else  
	  ModObj[i] += lp->lambda[j] * lp->lag_row[j][i];
      }
    }
    for(i = 1; i <= lp->columns; i++) {  
      set_mat(lp, 0, i, ModObj[i]);
    }
    rhsmod = 0;
    for(i = 0; i < lp->nr_lagrange; i++)
      if(lp->maximise)
	rhsmod += lp->lambda[i] * lp->lag_rhs[i];
      else
	rhsmod -= lp->lambda[i] * lp->lag_rhs[i];
 
    if(verbose) {
      fprintf(stderr, "Zub: %10g Zlb: %10g Step: %10g pie: %10g Feas %d\n",
	      (double)Zub, (double)Zlb, (double)Step, (double)pie, OrigFeas);
      for(i = 0; i < lp->nr_lagrange; i++)
	fprintf(stderr, "%3d SubGrad %10g lambda %10g\n", i,
		(double)SubGrad[i], (double)lp->lambda[i]);
    }

    if(verbose && lp->sum < 20)
      print_lp(lp);

    result = solve(lp);

    if(verbose && lp->sum < 20) { 
      print_solution(lp);
    }

    same_basis = TRUE;
    i = 1;
    while(same_basis && i < lp->rows) {
      same_basis = (old_bas[i] == lp->bas[i]);
      i++;
    }
    i = 1;
    while(same_basis && i < lp->sum) {
      same_basis=(old_lower[i] == lp->lower[i]);
      i++;
    }
    if(!same_basis) {
      memcpy(old_lower, lp->lower, (lp->sum+1) * sizeof(short));
      memcpy(old_bas, lp->bas, (lp->rows+1) * sizeof(int));
      pie *= 0.95;
    }

    if(verbose)
      fprintf(stderr, "result: %d  same basis: %d\n", result, same_basis);
      
    if(result == UNBOUNDED) {
      for(i = 1; i <= lp->columns; i++)
	fprintf(stderr, "%g ", (double)ModObj[i]);
      exit(EXIT_FAILURE);
    }

    if(result == FAILURE)
      status = FAILURE;

    if(result == INFEASIBLE)
      status = INFEASIBLE;
      
    SqrsumSubGrad = 0;
    for(i = 0; i < lp->nr_lagrange; i++) {
      SubGrad[i]= -lp->lag_rhs[i];
      for(j = 1; j <= lp->columns; j++)
	SubGrad[i] += lp->best_solution[lp->rows + j] * lp->lag_row[i][j];
      SqrsumSubGrad += SubGrad[i] * SubGrad[i];
    }

    OrigFeas = TRUE;
    for(i = 0; i < lp->nr_lagrange; i++)
      if(lp->lag_con_type[i]) {
	if(my_abs(SubGrad[i]) > lp->epsb)
	  OrigFeas = FALSE;
      }
      else if(SubGrad[i] > lp->epsb)
	OrigFeas = FALSE;

    if(OrigFeas) {
      AnyFeas = TRUE;
      Ztmp = 0;
      for(i = 1; i <= lp->columns; i++)
	Ztmp += lp->best_solution[lp->rows + i] * OrigObj[i];
      if((lp->maximise) && (Ztmp > Zlb)) {
	Zlb = Ztmp;
	for(i = 1; i <= lp->sum; i++)
	  BestFeasSol[i] = lp->best_solution[i];
	BestFeasSol[0] = Zlb;
	if(verbose)
	  fprintf(stderr, "Best feasible solution: %g\n", (double)Zlb);
      }
      else if(Ztmp < Zub) {
	Zub = Ztmp;
	for(i = 1; i <= lp->sum; i++)
	  BestFeasSol[i] = lp->best_solution[i];
	BestFeasSol[0] = Zub;
	if(verbose)
	  fprintf(stderr, "Best feasible solution: %g\n", (double)Zub);
      }
        }      

    if(lp->maximise)
      Zub = my_min(Zub, rhsmod + lp->best_solution[0]);
    else
      Zlb = my_max(Zlb, rhsmod + lp->best_solution[0]);

    if(my_abs(Zub-Zlb)<0.001) {  
      status = OPTIMAL;
    }
    Step = pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad;  
 
    for(i = 0; i < lp->nr_lagrange; i++) {
      lp->lambda[i] += Step * SubGrad[i];
      if(!lp->lag_con_type[i] && lp->lambda[i] < 0)
	lp->lambda[i] = 0;
    }
 
    if(citer == num_iter && status==RUNNING) {
      if(AnyFeas)
	status = FEAS_FOUND;
      else
	status = NO_FEAS_FOUND;
    }
  }

  for(i = 0; i <= lp->sum; i++)
    lp->best_solution[i] = BestFeasSol[i];
 
  for(i = 1; i <= lp->columns; i++)
    set_mat(lp, 0, i, OrigObj[i]);

  if(lp->maximise)
    lp->lag_bound = Zub;
  else
    lp->lag_bound = Zlb;
  free(BestFeasSol);
  free(SubGrad);
  free(OrigObj);
  free(ModObj);
  free(old_bas);
  free(old_lower);
  
  return(status);
}